eigen.cpp
Go to the documentation of this file.
#ifndef _MSC_VER
# include <itpp/config.h>
#else
# include <itpp/config_msvc.h>
#endif
#if defined(HAVE_LAPACK)
# include <itpp/base/algebra/lapack.h>
#endif
#include <itpp/base/algebra/eigen.h>
#include <itpp/base/converters.h>
namespace itpp
{
#if defined(HAVE_LAPACK)
bool eig_sym(const mat &A, vec &d, mat &V)
{
it_assert_debug(A.rows() == A.cols(), "eig_sym: Matrix is not symmetric");
// Test for symmetric?
char jobz = 'V', uplo = 'U';
int n, lda, lwork, info;
n = lda = A.rows();
lwork = std::max(1, 3 * n - 1); // This may be choosen better!
d.set_size(n, false);
vec work(lwork);
V = A; // The routine overwrites input matrix with eigenvectors
dsyev_(&jobz, &uplo, &n, V._data(), &lda, d._data(), work._data(), &lwork, &info);
return (info == 0);
}
bool eig_sym(const mat &A, vec &d)
{
it_assert_debug(A.rows() == A.cols(), "eig_sym: Matrix is not symmetric");
// Test for symmetric?
char jobz = 'N', uplo = 'U';
int n, lda, lwork, info;
n = lda = A.rows();
lwork = std::max(1, 3 * n - 1); // This may be choosen better!
d.set_size(n, false);
vec work(lwork);
mat B(A); // The routine overwrites input matrix
dsyev_(&jobz, &uplo, &n, B._data(), &lda, d._data(), work._data(), &lwork, &info);
return (info == 0);
}
bool eig_sym(const cmat &A, vec &d, cmat &V)
{
it_assert_debug(A.rows() == A.cols(), "eig_sym: Matrix is not hermitian");
// Test for symmetric?
char jobz = 'V', uplo = 'U';
int n, lda, lwork, info;
n = lda = A.rows();
lwork = std::max(1, 2 * n - 1); // This may be choosen better!
d.set_size(n, false);
cvec work(lwork);
vec rwork(std::max(1, 3*n - 2)); // This may be choosen better!
V = A; // The routine overwrites input matrix with eigenvectors
zheev_(&jobz, &uplo, &n, V._data(), &lda, d._data(), work._data(), &lwork, rwork._data(), &info);
return (info == 0);
}
bool eig_sym(const cmat &A, vec &d)
{
it_assert_debug(A.rows() == A.cols(), "eig_sym: Matrix is not hermitian");
// Test for symmetric?
char jobz = 'N', uplo = 'U';
int n, lda, lwork, info;
n = lda = A.rows();
lwork = std::max(1, 2 * n - 1); // This may be choosen better!
d.set_size(n, false);
cvec work(lwork);
vec rwork(std::max(1, 3*n - 2)); // This may be choosen better!
cmat B(A); // The routine overwrites input matrix
zheev_(&jobz, &uplo, &n, B._data(), &lda, d._data(), work._data(), &lwork, rwork._data(), &info);
return (info == 0);
}
// Non-symmetric matrix
bool eig(const mat &A, cvec &d, cmat &V)
{
it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square");
char jobvl = 'N', jobvr = 'V';
int n, lda, ldvl, ldvr, lwork, info;
n = lda = A.rows();
ldvl = 1;
ldvr = n;
lwork = std::max(1, 4 * n); // This may be choosen better!
vec work(lwork);
vec rwork(std::max(1, 2*n)); // This may be choosen better
vec wr(n), wi(n);
mat vl, vr(n, n);
mat B(A); // The routine overwrites input matrix
dgeev_(&jobvl, &jobvr, &n, B._data(), &lda, wr._data(), wi._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, &info);
d = to_cvec(wr, wi);
// 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 mat &A, cvec &d)
{
it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square");
char jobvl = 'N', jobvr = 'N';
int n, lda, ldvl, ldvr, lwork, info;
n = lda = A.rows();
ldvl = 1;
ldvr = 1;
lwork = std::max(1, 4 * n); // This may be choosen better!
vec work(lwork);
vec rwork(std::max(1, 2*n)); // This may be choosen better
vec wr(n), wi(n);
mat vl, vr;
mat B(A); // The routine overwrites input matrix
dgeev_(&jobvl, &jobvr, &n, B._data(), &lda, wr._data(), wi._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, &info);
d = to_cvec(wr, wi);
return (info == 0);
}
bool eig(const cmat &A, cvec &d, cmat &V)
{
it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square");
char jobvl = 'N', jobvr = 'V';
int n, lda, ldvl, ldvr, lwork, info;
n = lda = A.rows();
ldvl = 1;
ldvr = n;
lwork = std::max(1, 2 * n); // This may be choosen better!
d.set_size(n, false);
V.set_size(n, n, false);
cvec work(lwork);
vec rwork(std::max(1, 2*n)); // This may be choosen better!
cmat vl;
cmat B(A); // The routine overwrites input matrix
zgeev_(&jobvl, &jobvr, &n, B._data(), &lda, d._data(), vl._data(), &ldvl, V._data(), &ldvr, work._data(), &lwork, rwork._data(), &info);
return (info == 0);
}
bool eig(const cmat &A, cvec &d)
{
it_assert_debug(A.rows() == A.cols(), "eig: Matrix is not square");
char jobvl = 'N', jobvr = 'N';
int n, lda, ldvl, ldvr, lwork, info;
n = lda = A.rows();
ldvl = 1;
ldvr = 1;
lwork = std::max(1, 2 * n); // This may be choosen better!
d.set_size(n, false);
cvec work(lwork);
vec rwork(std::max(1, 2*n)); // This may be choosen better!
cmat vl, vr;
cmat B(A); // The routine overwrites input matrix
zgeev_(&jobvl, &jobvr, &n, B._data(), &lda, d._data(), vl._data(), &ldvl, vr._data(), &ldvr, work._data(), &lwork, rwork._data(), &info);
return (info == 0);
}
#else
bool eig_sym(const mat &A, vec &d, mat &V)
{
it_error("LAPACK library is needed to use eig_sym() function");
return false;
}
bool eig_sym(const mat &A, vec &d)
{
it_error("LAPACK library is needed to use eig_sym() function");
return false;
}
bool eig_sym(const cmat &A, vec &d, cmat &V)
{
it_error("LAPACK library is needed to use eig_sym() function");
return false;
}
bool eig_sym(const cmat &A, vec &d)
{
it_error("LAPACK library is needed to use eig_sym() function");
return false;
}
bool eig(const mat &A, cvec &d, cmat &V)
{
it_error("LAPACK library is needed to use eig() function");
return false;
}
bool eig(const mat &A, cvec &d)
{
it_error("LAPACK library is needed to use eig() function");
return false;
}
bool eig(const cmat &A, cvec &d, cmat &V)
{
it_error("LAPACK library is needed to use eig() function");
return false;
}
bool eig(const cmat &A, cvec &d)
{
it_error("LAPACK library is needed to use eig() function");
return false;
}
#endif // HAVE_LAPACK
vec eig_sym(const mat &A)
{
vec d;
eig_sym(A, d);
return d;
}
vec eig_sym(const cmat &A)
{
vec d;
eig_sym(A, d);
return d;
}
cvec eig(const mat &A)
{
cvec d;
eig(A, d);
return d;
}
cvec eig(const cmat &A)
{
cvec d;
eig(A, d);
return d;
}
} //namespace itpp
eigen.h
Definitions of eigenvalue decomposition functions.
itpp
itpp namespace
Definition: itmex.h:36
itpp::conj
cvec conj(const cvec &x)
Conjugate of complex value.
Definition: elem_math.h:226
it_assert_debug
#define it_assert_debug(t, s)
Abort if t is not true and NDEBUG is not defined.
Definition: itassert.h:107
itpp::eig
bool eig(const mat &A, cvec &d, cmat &V)
Calculates the eigenvalues and eigenvectors of a real non-symmetric matrix.
Definition: eigen.cpp:277
itpp::to_cvec
cvec to_cvec(const Vec< T > &v)
Converts a Vec<T> to cvec.
Definition: converters.h:107
itpp::max
T max(const Vec< T > &v)
Maximum value of vector.
Definition: min_max.h:45
it_error
#define it_error(s)
Abort unconditionally.
Definition: itassert.h:126
itpp::eig_sym
bool eig_sym(const mat &A, vec &d, mat &V)
Calculates the eigenvalues and eigenvectors of a symmetric real matrix.
Definition: eigen.cpp:252
converters.h
Definitions of converters between different vector and matrix types.

Generated on Mon Mar 23 2020 06:18:03 for IT++ by Doxygen 1.8.17