Commit 640d5071 authored by Michael Goette's avatar Michael Goette

added largest eigenvalue for arpack support

parent 9e3b7e46
Pipeline #1138 passed with stages
in 23 minutes and 14 seconds
......@@ -60,15 +60,18 @@ namespace xerus {
void solve_ev(double* const _x, const double* const _A, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, arpack::which const _ritz_option, int _info);
///@brief: Solves Ax = lambda*x for x, for the smallest _k eigenvalues
void solve_ev_smallest(double* const _x, const double* const _A, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info);
///@brief: Solves Ax = lambda*x for x, for the biggest _k eigenvalues
void solve_ev_biggest(double* const _x, const double* const _A, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info);
///@brief: Solves Ax = lambda*x for x, for the largest _k eigenvalues
void solve_ev_largest(double* const _x, const double* const _A, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info);
//TODO check if this can be simplified!!
///@brief: Solves Ax = lambda*x for x, this calls the Arpack Routine dsaupd
void solve_ev_special(double* const _x, const TensorNetwork& _op, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, arpack::which const _ritz_option, int _info);
///@brief: Solves Ax = lambda*x for x, for the smallest _k eigenvalues
///@brief: Solves Ax = lambda*x for x, for the smallest _k eigenvalues, takes Tensor Network as Operator
void solve_ev_smallest_special(double* const _x, const TensorNetwork& _op, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info);
//@brief: Solves Ax = lambda*x for x, for the largest _k eigenvalues, takes Tensor Network as Operator
void solve_ev_largest_special(double* const _x, const TensorNetwork& _op, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info);
......
......@@ -1024,7 +1024,7 @@ namespace xerus {
* @param _A input Operator A symmetric with respect to matrification
* @return the smallest eigenvalue
*/
value_t get_smallest_eigenpair(Tensor &_X, const Tensor &_A);
value_t get_smallest_eigenpair(Tensor &_X, Tensor &_A);
#ifdef ARPACK_LIBRARIES
......@@ -1032,23 +1032,25 @@ namespace xerus {
* @brief Solves the equation A*x = lambda*x for x and lambda and A symmetric. It calls the ARPACK routine dsaupd. It calculates the smallest algebraic values.
* @param _X Output Tensor for the result, the eigenvector of the smallest eigenvalue
* @param _A input Operator A in Tensor Format, symmetric with respect to matrification
* @param _smallest if true algorithm searches for smallest eigenvalue else algorithm searches for largest eigenvalue
* @param _initialize if true algorithm is randomly initialized, if false it takes _X as starting point for the lanczos iteration
* @param _miter maximal number of iterations of algorithm
* @param _eps tolerance for iterative algorithm
* @return smallest eigenvalue
*/
value_t get_smallest_eigenpair_iterative(Tensor& _X, const Tensor& _A, bool _initialize=true, const size_t _miter=1000, const double _eps=EPSILON);
value_t get_eigenpair_iterative(Tensor& _X, Tensor& _A, bool _smallest=false, bool _initialize=true, const size_t _miter=1000, const double _eps=EPSILON);
/**
* @brief Solves the equation A*x = lambda*x for x and lambda and A symmetric. It calls the ARPACK routine dsaupd. It calculates the smallest algebraic values.
* @param _X Output Tensor for the result, the eigenvector of the smallest eigenvalue
* @param _A input Operator A in Tensor Network Format, symmetric with respect to matrification
* @param _smallest if true algorithm searches for smallest eigenvalue else algorithm searches for largest eigenvalue
* @param _initialize if true algorithm is randomly initialized, if false it takes _X as starting point for the lanczos iteration
* @param _miter maximal number of iterations of algorithm
* @param _eps tolerance for iterative algorithm
* @return smallest eigenvalue
*/
value_t get_smallest_eigenpair_iterative(Tensor& _X, const TensorNetwork& _A, bool _initialize=true, const size_t _miter=1000, const double _eps=EPSILON);
value_t get_eigenpair_iterative(Tensor& _X, const TensorNetwork& _A, bool _smallest=false, bool _initialize=true, const size_t _miter=1000, const double _eps=EPSILON);
#endif
......
......@@ -95,10 +95,10 @@ static misc::UnitTest tensor_solve_smallest_ev_iterative("Tensor", "get_smallest
Tensor x1({4,3});
Tensor x11({4,3});
value_t lambda1 = get_smallest_eigenpair_iterative(x1,A1, true, 10000, 1e-8);
value_t lambda11 = get_smallest_eigenpair_iterative(x11,A11, true, 10000, 1e-8);
value_t lambda1 = get_eigenpair_iterative(x1,A1, true, true, 10000, 1e-8);
value_t lambda11 = get_eigenpair_iterative(x11,A11, true, true, 10000, 1e-8);
MTEST(frob_norm(A1(i,j,k,l)*x1(k,l) - lambda1* x1(i,j)) < 1e-8, frob_norm(A1(i,j,k,l)*x1(k,l) - lambda1 * x1(i,j)));
MTEST(frob_norm(A11(i,j,k,l)*x1(k,l) - lambda11* x11(i,j)) < 1e-8, frob_norm(A11(i,j,k,l)*x11(k,l) - lambda11 * x11(i,j)));
MTEST(frob_norm(A11(i,j,k,l)*x11(k,l) - lambda11* x11(i,j)) < 1e-8, frob_norm(A11(i,j,k,l)*x11(k,l) - lambda11 * x11(i,j)));
Tensor A2 = Tensor::random({2,2,2,2,2,2,2,2});
A22(i/2,j/2) = A2(i/2, k/2) * A2(j/2, k/2);
......@@ -106,8 +106,8 @@ static misc::UnitTest tensor_solve_smallest_ev_iterative("Tensor", "get_smallest
Tensor x2({2,2,2,2});
Tensor x22({2,2,2,2});
value_t lambda2 = get_smallest_eigenpair_iterative(x2,A2, true, 10000, 1e-8);
value_t lambda22 = get_smallest_eigenpair_iterative(x22,A22, true, 10000, 1e-8);
value_t lambda2 = get_eigenpair_iterative(x2,A2, true, true, 10000, 1e-8);
value_t lambda22 = get_eigenpair_iterative(x22,A22, true, true, 10000, 1e-8);
MTEST(frob_norm(A2(i,j,k,l,m,n,o,p)*x2(m,n,o,p) - lambda2 * x2(i,j,k,l)) < 1e-8, frob_norm(A2(i,j,k,l,m,n,o,p)*x2(m,n,o,p) - lambda2 * x2(i,j,k,l)));
MTEST(frob_norm(A22(i,j,k,l,m,n,o,p)*x22(m,n,o,p) - lambda22 * x22(i,j,k,l)) < 1e-8, frob_norm(A22(i,j,k,l,m,n,o,p)*x22(m,n,o,p) - lambda22 * x22(i,j,k,l)));
......@@ -117,8 +117,45 @@ static misc::UnitTest tensor_solve_smallest_ev_iterative("Tensor", "get_smallest
Tensor x3 = Tensor::random({2,2,2,2});
Tensor x33 = Tensor::random({2,2,2,2});
value_t lambda3 = get_smallest_eigenpair_iterative(x3,A3, false, 10000, 1e-8);
value_t lambda33 = get_smallest_eigenpair_iterative(x33,A33, false, 10000, 1e-8);
value_t lambda3 = get_eigenpair_iterative(x3,A3, true, false, 10000, 1e-8);
value_t lambda33 = get_eigenpair_iterative(x33,A33, true, false, 10000, 1e-8);
MTEST(frob_norm(A3(i,j,k,l,m,n,o,p)*x3(m,n,o,p) - lambda3 * x3(i,j,k,l)) < 1e-8, frob_norm(A3(i,j,k,l,m,n,o,p)*x3(m,n,o,p) - lambda3 * x3(i,j,k,l)));
MTEST(frob_norm(A33(i,j,k,l,m,n,o,p)*x33(m,n,o,p) - lambda33 * x33(i,j,k,l)) < 1e-8, frob_norm(A33(i,j,k,l,m,n,o,p)*x33(m,n,o,p) - lambda33 * x33(i,j,k,l)));
});
static misc::UnitTest tensor_solve_largest_ev_iterative("Tensor", "get_largest_eigenvalue_of_Matrix_iterative", [](){
Index i,j,k,l,m,n,o,p;
TensorNetwork A11,A22,A33;
Tensor A1 = Tensor::random({4,3,4,3}) + (-1)*Tensor::identity({4,3,4,3});
A11(i/2,j/2) = A1(i/2, k/2) * A1(j/2, k/2);
A1(i/2,j/2) = A1(i/2, k/2) * A1(j/2, k/2);
Tensor x1({4,3});
Tensor x11({4,3});
value_t lambda1 = get_eigenpair_iterative(x1,A1, false, true, 10000, 1e-10);
value_t lambda11 = get_eigenpair_iterative(x11,A11, false, true, 10000, 1e-10);
MTEST(frob_norm(A1(i,j,k,l)*x1(k,l) - lambda1* x1(i,j)) < 1e-8, frob_norm(A1(i,j,k,l)*x1(k,l) - lambda1 * x1(i,j)));
MTEST(frob_norm(A11(i,j,k,l)*x11(k,l) - lambda11* x11(i,j)) < 1e-8, frob_norm(A11(i,j,k,l)*x11(k,l) - lambda11 * x11(i,j)));
Tensor A2 = Tensor::random({2,2,2,2,2,2,2,2});
A22(i/2,j/2) = A2(i/2, k/2) * A2(j/2, k/2);
A2(i/2,j/2) = A2(i/2, k/2) * A2(j/2, k/2);
Tensor x2({2,2,2,2});
Tensor x22({2,2,2,2});
value_t lambda2 = get_eigenpair_iterative(x2,A2, false, true, 10000, 1e-10);
value_t lambda22 = get_eigenpair_iterative(x22,A22, false, true, 10000, 1e-10);
MTEST(frob_norm(A2(i,j,k,l,m,n,o,p)*x2(m,n,o,p) - lambda2 * x2(i,j,k,l)) < 1e-8, frob_norm(A2(i,j,k,l,m,n,o,p)*x2(m,n,o,p) - lambda2 * x2(i,j,k,l)));
MTEST(frob_norm(A22(i,j,k,l,m,n,o,p)*x22(m,n,o,p) - lambda22 * x22(i,j,k,l)) < 1e-8, frob_norm(A22(i,j,k,l,m,n,o,p)*x22(m,n,o,p) - lambda22 * x22(i,j,k,l)));
Tensor A3 = Tensor::random({2,2,2,2,2,2,2,2});
A33(i/2,j/2) = A3(i/2, k/2) * A3(j/2, k/2);
A3(i/2,j/2) = A3(i/2, k/2) * A3(j/2, k/2);
Tensor x3 = Tensor::random({2,2,2,2});
Tensor x33 = Tensor::random({2,2,2,2});
value_t lambda3 = get_eigenpair_iterative(x3,A3, false, false, 10000, 1e-10);
value_t lambda33 = get_eigenpair_iterative(x33,A33, false, false, 10000, 1e-10);
MTEST(frob_norm(A3(i,j,k,l,m,n,o,p)*x3(m,n,o,p) - lambda3 * x3(i,j,k,l)) < 1e-8, frob_norm(A3(i,j,k,l,m,n,o,p)*x3(m,n,o,p) - lambda3 * x3(i,j,k,l)));
MTEST(frob_norm(A33(i,j,k,l,m,n,o,p)*x33(m,n,o,p) - lambda33 * x33(i,j,k,l)) < 1e-8, frob_norm(A33(i,j,k,l,m,n,o,p)*x33(m,n,o,p) - lambda33 * x33(i,j,k,l)));
});
......
......@@ -292,13 +292,16 @@ namespace xerus {
solve_ev (_x, _A, _ev, _k, _n, _resid, _maxiter, _eps, arpack::which::smallest_algebraic, _info);
}
void solve_ev_biggest(double* const _x, const double* const _A, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info) {
void solve_ev_largest(double* const _x, const double* const _A, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info) {
solve_ev (_x, _A, _ev, _k, _n, _resid, _maxiter, _eps, arpack::which::largest_algebraic, _info);
}
void solve_ev_smallest_special(double* const _x, const TensorNetwork& _op, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info) {
solve_ev_special (_x, _op, _ev, _k, _n, _resid, _maxiter, _eps, arpack::which::smallest_algebraic, _info);
}
void solve_ev_largest_special(double* const _x, const TensorNetwork& _op, double* const _ev, const size_t _k, const size_t _n, double* const _resid, const size_t _maxiter, const double _eps, int _info) {
solve_ev_special (_x, _op, _ev, _k, _n, _resid, _maxiter, _eps, arpack::which::largest_algebraic, _info);
}
} // namespace arpackWrapper
......
......@@ -1712,7 +1712,7 @@ namespace xerus {
_X.factor = _B.factor / _A.factor;
}
double get_smallest_eigenpair(Tensor& _X, const Tensor& _A) {
double get_smallest_eigenpair(Tensor& _X, Tensor& _A) {
REQUIRE(_A.is_dense(), "for now only dense is implemented"); //TODO implement sparse
REQUIRE(&_X != &_A, "Not supportet yet");
REQUIRE(_A.degree() % 2 == 0, "The tensor A needs to be an operator, i.e. has even degree");
......@@ -1742,7 +1742,7 @@ namespace xerus {
rev.get(), // right eigenvectors
re.get(), // real eigenvalues
im.get(), // imaginary eigenvalues
_A.get_unsanitized_dense_data(), n);
_A.get_dense_data(), n);
//get eigenvalue with minimal real part
double ev = re[0];
......@@ -1763,7 +1763,7 @@ namespace xerus {
}
#ifdef ARPACK_LIBRARIES
value_t get_smallest_eigenpair_iterative(Tensor& _X, const Tensor& _A, bool _initialize, const size_t _miter, const double _eps) {
value_t get_eigenpair_iterative(Tensor& _X, Tensor& _A, bool _smallest, bool _initialize, const size_t _miter, const double _eps) {
REQUIRE(_A.is_dense(), "for now only dense is implemented"); //TODO implement sparse
REQUIRE(&_X != &_A, "Not supportet yet");
REQUIRE(_A.degree() % 2 == 0, "The tensor A needs to be an operator, i.e. has even degree");
......@@ -1791,17 +1791,28 @@ namespace xerus {
std::unique_ptr<double[]> rev(new double[n]); // right eigenvectors
std::unique_ptr<double[]> res(new double[n]); // residual
if (info > 0)
misc::copy(res.get(), _X.get_unsanitized_dense_data(), n);
misc::copy(res.get(), _X.get_dense_data(), n);
// Note that A is dense here
arpackWrapper::solve_ev_smallest(
rev.get(), // right ritz vectors
_A.get_unsanitized_dense_data(),
ev.get(), 1, n,
res.get(),
_miter,
_eps, info
);
if (_smallest){
arpackWrapper::solve_ev_smallest(
rev.get(), // right ritz vectors
_A.get_dense_data(),
ev.get(), 1, n,
res.get(),
_miter,
_eps, info
);
} else {
arpackWrapper::solve_ev_largest(
rev.get(), // right ritz vectors
_A.get_dense_data(),
ev.get(), 1, n,
res.get(),
_miter,
_eps, info
);
}
// eigenvector of smallest eigenvalue
auto tmpX = _X.override_dense_data();
......@@ -1810,7 +1821,7 @@ namespace xerus {
return ev[0];
}
value_t get_smallest_eigenpair_iterative(Tensor& _X, const TensorNetwork& _A, bool _initialize, const size_t _miter, const double _eps) {
value_t get_eigenpair_iterative(Tensor& _X, const TensorNetwork& _A, bool _smallest, bool _initialize, const size_t _miter, const double _eps) {
//REQUIRE(&_X != &_A, "Not supportet yet");
REQUIRE(_eps > 0 && _eps < 1, "epsilon must be betweeen 0 and 1, given " << _eps);
REQUIRE(_A.degree() % 2 == 0, "operator degree must be positive");
......@@ -1838,16 +1849,27 @@ namespace xerus {
misc::copy(res.get(), _X.get_unsanitized_dense_data(), n);
// Note that A is dense here
arpackWrapper::solve_ev_smallest_special(
rev.get(), // right ritz vectors
_A,
ev.get(), 1, n,
res.get(),
_miter,
_eps, info
);
if (_smallest){
arpackWrapper::solve_ev_smallest_special(
rev.get(), // right ritz vectors
_A,
ev.get(), 1, n,
res.get(),
_miter,
_eps, info
);
} else {
arpackWrapper::solve_ev_largest_special(
rev.get(), // right ritz vectors
_A,
ev.get(), 1, n,
res.get(),
_miter,
_eps, info
);
}
// eigenvector of smallest eigenvalue
// eigenvector of eigenvalue
auto tmpX = _X.override_dense_data();
for (size_t i = 0; i < n; ++i)
tmpX[i] = rev[i];
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment