--- lapack.h.old 2009-08-06 11:26:52.000000000 +0800 +++ lapack.h 2009-08-06 16:14:37.000000000 +0800 @@ -90,6 +90,9 @@ #define dgees_ dgees #define zgees_ zgees +#define dggev_ dggev +#define zggevx_ zggevx + #endif // #ifdef HAVE_MKL // Exists in ATLAS @@ -161,6 +164,17 @@ int *ldvl, std::complex *vr, int *ldvr, std::complex *work, int *lwork, double *rwork, int *info); + /* Eigenvalues and eigenvectors of a generalized matrices A and B*/ + void dggev_(char *jobvl, char *jobvr, int *n, double *a, int *lda, double *b, + int *ldb, double *alphar, double *alphai, double *beta, + double *vl, int *ldvl, double *vr, int *ldvr, + double *work, int *lwork, int *info); + void zggev_(char *jobvl, char *jobvr, int *n, std::complex *a, + int *lda, std::complex *b, int *ldb, std::complex *alpha, + std::complex *beta, std::complex *vl, + int *ldvl, std::complex *vr, int *ldvr, + std::complex *work, int *lwork, double *rwork, int *info); + // In ATLAS /* Cholesky factorization */ void dpotrf_(char *uplo, int *n, double *a, int *lda, int *info); --- eigen.h.old 2009-08-12 15:34:50.000000000 +0800 +++ eigen.h 2009-08-12 15:36:38.000000000 +0800 @@ -140,7 +140,7 @@ \f$\mathbf{v}_i, \: i=0, \ldots, n-1\f$ of the real \f$n \times n\f$ matrix \f$\mathbf{A}\f$ satisfies \f[ - \mathbf{A} \mathbf{v}_i = d_i \mathbf{v}_i\: i=0, \ldots, n-1. + \mathbf{A} \mathbf{v}_i = d_i \mathbf{B} \mathbf{v}_i\: i=0, \ldots, n-1. \f] The eigenvectors are the columns of the matrix V. True is returned if the calculation was successful. Otherwise false. @@ -151,6 +151,40 @@ /*! \ingroup matrixdecomp + \brief Calculates the generalized eigenvalues and eigenvectors of a complex non-symmetric matrix + + The Eigenvalues \f$\mathbf{d}(d_0, d_1, \ldots, d_{n-1})\f$ and the eigenvectors + \f$\mathbf{v}_i, \: i=0, \ldots, n-1\f$ of the real \f$n \times n\f$ + matrix \f$\mathbf{A}\f$ satisfies + \f[ + \mathbf{A} \mathbf{v}_i = d_i \mathbf{B} \mathbf{v}_i\: i=0, \ldots, n-1. + \f] + The eigenvectors are the columns of the matrix V. + True is returned if the calculation was successful. Otherwise false. + + Uses the LAPACK routine ZGGEV. +*/ +bool eig(const mat &A, const mat &B, cvec &d, cmat &V); + +/*! + \ingroup matrixdecomp + \brief Calculates the generalized eigenvalues and eigenvectors of a real non-symmetric matrix + + The Eigenvalues \f$\mathbf{d}(d_0, d_1, \ldots, d_{n-1})\f$ and the eigenvectors + \f$\mathbf{v}_i, \: i=0, \ldots, n-1\f$ of the real \f$n \times n\f$ + matrix \f$\mathbf{A}\f$ satisfies + \f[ + \mathbf{A} \mathbf{v}_i = d_i \mathbf{v}_i\: i=0, \ldots, n-1. + \f] + The eigenvectors are the columns of the matrix V. + True is returned if the calculation was successful. Otherwise false. + + Uses the LAPACK routine DGEEV. +*/ +bool eig(const cmat &A, const cmat &B, cvec &d, cmat &V); + +/*! + \ingroup matrixdecomp \brief Calculates the eigenvalues of a real non-symmetric matrix The Eigenvalues \f$\mathbf{d}(d_0, d_1, \ldots, d_{n-1})\f$ and the eigenvectors --- eigen.cpp.old 2009-08-06 11:39:13.000000000 +0800 +++ eigen.cpp 2009-08-12 14:29:14.000000000 +0800 @@ -26,7 +26,7 @@ * * ------------------------------------------------------------------------- */ - +#include #ifndef _MSC_VER # include #else @@ -174,6 +174,97 @@ } // Non-symmetric matrix +bool eig(const mat &A, const mat &B, cvec &d, cmat &V) +{ + it_assert_debug(A.rows() == A.cols(), "eig: Matrix A is not square"); + it_assert_debug(B.rows() == B.cols(), "eig: Matrix B is not square"); + + char jobvl = 'N', jobvr = 'V'; + int n, lda, ldb, ldvl, ldvr, lwork, info; + n = lda = A.rows(); + ldb = B.rows(); + ldvl = 1; + ldvr = n; + lwork = std::max(1, 8 * n); // This may be choosen better! + + vec work(lwork); + vec rwork(std::max(1, 2*n)); // This may be choosen better + vec alphar(n), alphai(n), beta(n); + mat vl, vr(n, n); + + mat C(A); // The routine overwrites input matrix + mat D(B); // The routine overwrites input matrix + + dggev_(&jobvl, &jobvr, &n, C._data(), &lda, D._data(), &ldb, alphar._data(), alphai._data(), beta._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, &info); + + d = to_cvec(alphar, alphai); + d = elem_div(d,to_cvec(beta)); + + // Fix V + V.set_size(n, n, false); + for (int j = 0; j < n; j++) { + // if d(j) and d(j+1) are complex conjugate pairs, treat special + if ((j < n - 1) && d(j) == std::conj(d(j + 1))) { + V.set_col(j, to_cvec(vr.get_col(j), vr.get_col(j + 1))); + V.set_col(j + 1, to_cvec(vr.get_col(j), -vr.get_col(j + 1))); + j++; + } + else { + V.set_col(j, to_cvec(vr.get_col(j))); + } + } + + return (info == 0); +} + +// Non-symmetric matrix +bool eig(const cmat &A, const cmat &B, cvec &d, cmat &V) +{ + it_assert_debug(A.rows() == A.cols(), "eig: Matrix A is not square"); + it_assert_debug(B.rows() == B.cols(), "eig: Matrix B is not square"); + + char jobvl = 'N', jobvr = 'V'; + int n, lda, ldb, ldvl, ldvr, lwork, info; + n = lda = A.rows(); + ldb = B.rows(); + ldvl = 1; + ldvr = n; + lwork = std::max(1, 2*n); // This may be choosen better! + + cvec work(lwork); + vec rwork(8*n); // This may be choosen better + cvec alpha(n), beta(n); + cmat vl(1,1), vr(n, n); + + cmat C(A); // The routine overwrites input matrix + cmat D(B); // The routine overwrites input matrix + + zggev_(&jobvl, &jobvr, &n, C._data(), &lda, D._data(), &ldb, alpha._data(), beta._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, rwork._data(), &info); + + d = elem_div(alpha,beta); + + // Fix V + V.set_size(n, n, false); + V=vr; + /* for (int j = 0; j < n; j++) { + + // if d(j) and d(j+1) are complex conjugate pairs, treat special + if ((j < n - 1) && d(j) == std::conj(d(j + 1))) { + V.set_col(j, to_cvec(vr.get_col(j), vr.get_col(j + 1))); + V.set_col(j + 1, to_cvec(vr.get_col(j), -vr.get_col(j + 1))); + j++; + } + else { + V.set_col(j, to_cvec(vr.get_col(j))); + } + + V.set_col(j,vr.get_col(j)); + } +*/ + return (info == 0); +} + +// Non-symmetric matrix bool eig(const mat &A, cvec &d) { it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square");