|
Periodic Eigendecomposition
|
00001 00011 #ifndef PED_H 00012 #define PED_H 00013 00014 #include <Eigen/Dense> 00015 #include <Eigen/Sparse> 00016 #include <Eigen/Eigenvalues> 00017 #include <Eigen/SparseLU> 00018 #include <vector> 00019 #include <utility> 00020 #include <tuple> 00021 00022 00023 using std::vector; 00024 using std::pair; using std::make_pair; 00025 using Eigen::MatrixXd; 00026 using Eigen::VectorXd; 00027 using Eigen::Matrix2d; using Eigen::Vector2d; 00028 using Eigen::Matrix2cd; using Eigen::Vector2cd; 00029 using Eigen::MatrixXcd; using Eigen::VectorXcd; 00030 using Eigen::Map; using Eigen::Ref; 00031 using Eigen::EigenSolver; 00032 using Eigen::Triplet; 00033 using Eigen::SparseMatrix; 00034 using Eigen::SparseLU; 00035 using Eigen::COLAMDOrdering; 00036 00037 /*============================================================ * 00038 * Class : periodic Eigendecomposition * 00039 *============================================================ */ 00040 class PED{ 00041 00042 public: 00043 MatrixXd 00044 EigVals(MatrixXd &J, const int MaxN = 100, 00045 const double tol = 1e-16 , bool Print = true); 00046 pair<MatrixXd, MatrixXd> 00047 EigVecs(MatrixXd &J, const int MaxN = 100, 00048 const double tol = 1e-16, bool Print = true); 00049 std::tuple<MatrixXd, vector<int>, MatrixXd> 00050 eigenvalues(MatrixXd &J, const int MaxN = 100, 00051 const double tol = 1e-16, bool Print = true); 00052 pair<MatrixXd, vector<int> > 00053 PerSchur(MatrixXd &J, const int MaxN = 100, 00054 const double tol = 1e-16, bool Print = true); 00055 MatrixXd 00056 HessTrian(MatrixXd &G); 00057 vector<int> 00058 PeriodicQR(MatrixXd &J, MatrixXd &Q, const int L, const int U, 00059 const int MaxN, const double tol, bool Print); 00060 //protected: 00061 void Givens(Ref<MatrixXd> A, Ref<MatrixXd> B, Ref<MatrixXd> C, 00062 const int k); 00063 void Givens(Ref<MatrixXd> A, Ref<MatrixXd> B, Ref<MatrixXd> C, 00064 const Vector2d &v, const int k); 00065 void HouseHolder(Ref<MatrixXd> A, Ref<MatrixXd> B, Ref<MatrixXd> C, const int k, 00066 bool subDiag = false); 00067 void GivensOneIter(MatrixXd &J, MatrixXd &Q, const Vector2d &v, 00068 const int L, const int U); 00069 void GivensOneRound(MatrixXd &J, MatrixXd &Q, const Vector2d &v, 00070 const int k); 00071 vector<int> 00072 checkSubdiagZero(const Ref<const MatrixXd> &J0, const int L, 00073 const int U,const double tol); 00074 vector<int> 00075 padEnds(const vector<int> &v, const int &left, const int &right); 00076 double 00077 deltaMat2(const Matrix2d &A); 00078 Vector2d 00079 vecMat2(const Matrix2d &A); 00080 pair<Vector2d, Matrix2d> 00081 complexEigsMat2(const Matrix2d &A); 00082 int 00083 sgn(const double &num); 00084 pair<double, int> 00085 product1dDiag(const MatrixXd &J, const int k); 00086 pair<Matrix2d, double> 00087 product2dDiag(const MatrixXd &J, const int k); 00088 vector<int> 00089 realIndex(const vector<int> &complexIndex, const int N); 00090 vector<Triplet<double> > 00091 triDenseMat(const Ref<const MatrixXd> &A, const size_t M = 0, 00092 const size_t N = 0); 00093 vector<Triplet<double> > 00094 triDenseMatKron(const size_t I, const Ref<const MatrixXd> &A, 00095 const size_t M = 0, const size_t N = 0); 00096 vector<Triplet<double> > 00097 triDiagMat(const size_t n, const double x, 00098 const size_t M = 0, const size_t N = 0 ); 00099 vector<Triplet<double> > 00100 triDiagMatKron(const size_t I, const Ref<const MatrixXd> &A, 00101 const size_t M = 0, const size_t N = 0 ); 00102 pair<SparseMatrix<double>, VectorXd> 00103 PerSylvester(const MatrixXd &J, const int &P, 00104 const bool &isReal, const bool &Print); 00105 MatrixXd 00106 oneEigVec(const MatrixXd &J, const int &P, 00107 const bool &isReal, const bool &Print); 00108 void 00109 fixPhase(MatrixXd &EigVecs, const VectorXd &realComplexIndex); 00110 void 00111 reverseOrder(MatrixXd &J); 00112 void 00113 reverseOrderSize(MatrixXd &J); 00114 }; 00115 00116 #endif /* PED_H */ 00117
1.7.6.1
