Commit 79860180 authored by Michael Goette's avatar Michael Goette

changed interface for eigenvalue solvers for smallest eigenpair and added some...

changed interface for eigenvalue solvers for smallest eigenpair and added some tests for iterative eigenvalue solver
parent c02abce3
Pipeline #1073 failed with stages
in 13 minutes and 19 seconds
......@@ -1019,25 +1019,35 @@ namespace xerus {
/**
* @brief Solves the equation A*x = lambda*x for x and lambda and A symmetric. It calls the LAPACK routine DGEEV. It calculates the eigenvalue with smallest real part.
* @param _X Output Tensor for the result
* @param _X Output Tensor for the result, this is the eigenvector in the appropriate tensor format
* @param _A input Operator A symmetric with respect to matrification
* @return the smallest eigenvalue
*/
double get_smallest_eigenvalue(Tensor &_X, const Tensor &_A);
value_t get_smallest_eigenpair(Tensor &_X, const Tensor &_A);
#ifdef ARPACK_LIBRARIES
/**
* @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 algerbaic values.
* @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 symmetric with respect to matrification
* @param _ev array holding smallest eigenvalue (ev[0]) after calculation
* @param _info if 0 algorithm is randomly initialized, if > 0 it takes _X as starting point
* @param _A input Operator A in Tensor Format, symmetric with respect to matrification
* @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);
/**
* @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 _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
*/
void get_smallest_eigenvalue_iterative(Tensor& _X, const Tensor& _A, double* const _ev, int _info, const size_t _miter, const double _eps);
void get_smallest_eigenvalue_iterative(Tensor& _X, const TensorNetwork& _op, double* const _ev, int _info, const size_t _miter, const double _eps);
value_t get_smallest_eigenpair_iterative(Tensor& _X, const TensorNetwork& _A, bool _initialize=true, const size_t _miter=1000, const double _eps=EPSILON);
#endif
......
......@@ -76,45 +76,51 @@ static misc::UnitTest tensor_solve_smallest_ev("Tensor", "get smallest eigenvalu
Tensor A1 = Tensor::random({4,3,4,3});
Tensor x1({4,3});
double lambda1 = get_smallest_eigenvalue(x1,A1);
value_t lambda1 = get_smallest_eigenpair(x1,A1);
TEST(frob_norm(A1(i,j,k,l)*x1(k,l) - lambda1 * x1(i,j)) < 1e-13);
Tensor A2 = Tensor::random({4,3,4,5,4,3,4,5});
Tensor x2({4,3,4,5});
double lambda2 = get_smallest_eigenvalue(x2,A2);
value_t lambda2 = get_smallest_eigenpair(x2,A2);
TEST(frob_norm(A2(i,j,k,l,m,n,o,p)*x2(m,n,o,p) - lambda2 * x2(i,j,k,l)) < 5e-12);
});
#ifdef ARPACK_LIBRARIES
static misc::UnitTest tensor_solve_smallest_ev_iterative("Tensor", "get smallest 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});
std::unique_ptr<double[]> ev1(new double[1]); // real eigenvalues
std::unique_ptr<double[]> ev2(new double[1]); // real eigenvalues
std::unique_ptr<double[]> ev3(new double[1]); // real eigenvalues
Tensor x11({4,3});
get_smallest_eigenvalue_iterative(x1,A1,ev1.get(), 0, 10000, 1e-8);
double lambda1 = ev1[0];
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);
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)));
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});
get_smallest_eigenvalue_iterative(x2,A2,ev2.get(), 0, 10000, 1e-8);
double lambda2 = ev2[0];
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);
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});
get_smallest_eigenvalue_iterative(x3,A3,ev3.get(), 1, 10000, 1e-8);
double lambda3 = ev3[0];
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);
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)));
});
#endif
......
......@@ -1715,7 +1715,7 @@ namespace xerus {
_X.factor = _B.factor / _A.factor;
}
double get_smallest_eigenvalue(Tensor& _X, const Tensor& _A) {
double get_smallest_eigenpair(Tensor& _X, const 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");
......@@ -1766,14 +1766,16 @@ namespace xerus {
}
#ifdef ARPACK_LIBRARIES
void get_smallest_eigenvalue_iterative(Tensor& _X, const Tensor& _A, double* const _ev, int _info, const size_t _miter, const double _eps) {
value_t get_smallest_eigenpair_iterative(Tensor& _X, const Tensor& _A, 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");
REQUIRE(_eps > 0 && _eps < 1, "epsilon must be betweeen 0 and 1, given " << _eps);
const size_t degN = _A.degree() / 2;
int info = _info;
int info = _initialize ? 0 : 1;
std::unique_ptr<value_t[]> ev(new value_t[1]); // resulting eigenvalue
// Calculate multDimensions
const size_t m = misc::product(_A.dimensions, 0, degN);
const size_t n = misc::product(_A.dimensions, degN, 2*degN);
......@@ -1798,7 +1800,7 @@ namespace xerus {
arpackWrapper::solve_ev_smallest(
rev.get(), // right ritz vectors
_A.get_unsanitized_dense_data(),
_ev, 1, n,
ev.get(), 1, n,
res.get(),
_miter,
_eps, info
......@@ -1808,26 +1810,27 @@ namespace xerus {
auto tmpX = _X.override_dense_data();
for (size_t i = 0; i < n; ++i)
tmpX[i] = rev[i];
return;
return ev[0];
}
void get_smallest_eigenvalue_iterative(Tensor& _X, const TensorNetwork& _op, double* const _ev, int _info, const size_t _miter, const double _eps) {
value_t get_smallest_eigenpair_iterative(Tensor& _X, const TensorNetwork& _A, 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(_op.degree() % 2 == 0, "operator degree must be positive");
REQUIRE(_A.degree() % 2 == 0, "operator degree must be positive");
const size_t degN = _op.degree() / 2;
int info = _info;
const size_t degN = _A.degree() / 2;
int info = _initialize ? 0 : 1;
std::unique_ptr<value_t[]> ev(new value_t[1]); // resulting eigenvalue
// Calculate multDimensions
const size_t m = misc::product(_op.dimensions, 0, degN);
const size_t n = misc::product(_op.dimensions, degN, 2*degN);
const size_t m = misc::product(_A.dimensions, 0, degN);
const size_t n = misc::product(_A.dimensions, degN, 2*degN);
XERUS_REQUIRE(m == n, "the dimensions of A do not agree, m != n, m x n = " << m << "x" << n);
// Make sure X has right dimensions
if( _X.degree() != degN || !std::equal(_X.dimensions.begin(), _X.dimensions.begin() + degN, _op.dimensions.begin() + degN))
if( _X.degree() != degN || !std::equal(_X.dimensions.begin(), _X.dimensions.begin() + degN, _A.dimensions.begin() + degN))
{
Tensor::DimensionTuple newDimX;
newDimX.insert(newDimX.end(), _op.dimensions.begin()+degN, _op.dimensions.end());
newDimX.insert(newDimX.end(), _A.dimensions.begin()+degN, _A.dimensions.end());
_X.reset(std::move(newDimX), Tensor::Representation::Dense, Tensor::Initialisation::None);
info = 0; // info must be 0 if X is not random
}
......@@ -1840,8 +1843,8 @@ namespace xerus {
// Note that A is dense here
arpackWrapper::solve_ev_smallest_special(
rev.get(), // right ritz vectors
_op,
_ev, 1, n,
_A,
ev.get(), 1, n,
res.get(),
_miter,
_eps, info
......@@ -1851,7 +1854,7 @@ namespace xerus {
auto tmpX = _X.override_dense_data();
for (size_t i = 0; i < n; ++i)
tmpX[i] = rev[i];
return;
return ev[0];
}
#endif
......
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