Commit efa261c1 authored by Michael Goette's avatar Michael Goette

with macrtos but not zes working as rexpected

parent ab43eb80
Pipeline #968 failed with stages
in 3 minutes
......@@ -131,6 +131,9 @@ LAPACK_LIBRARIES = -llapacke -llapack # Standard Lapack +
# Uncomment or add the appropriate CXSparse library
SUITESPARSE = -lcholmod -lspqr
# (optional) Uncomment if not needed, iterative eigenvalue solver ARPACK-ng, see https://github.com/opencollab/arpack-ng
# ARPACK_LIBRARIES = -larpack
# Uncomment or add the appropriate boost python library
BOOST_PYTHON2 = -lboost_python-py27
BOOST_PYTHON3 = -lboost_python-py35
......
......@@ -21,9 +21,28 @@
* @file
* @brief Header file for the arpack wrapper functions.
*/
#ifdef ARPACK_LIBRARIES
#pragma once
#include <complex.h>
// fix for non standard-conform complex implementation
#undef I
#ifdef __has_include
#if __has_include("arpack.hpp")
#include "arpack.hpp"
#elif __has_include("arpack/arpack.hpp")
#include "arpack/arpack.hpp"
#else
#pragma error no arpack found
#endif
#else
#include "arpack/arpack.hpp"
#endif
#include "misc/standard.h"
#include <memory>
......@@ -36,8 +55,19 @@ namespace xerus {
*/
namespace arpackWrapper {
///@brief: Solves Ax = lambda*x for x, this calls the Arpack Routine dsaupd
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);
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);
}
}
#endif
......@@ -1023,7 +1023,20 @@ namespace xerus {
*/
double get_smallest_eigenvalue(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.
* @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 _miter maximal number of iterations of algorithm
* @param _eps tolerance for iterative algorithm
*/
void get_smallest_eigenvalue_iterative(Tensor& _X, const Tensor& _A, double* const _ev, int _info, const size_t _miter, const double _eps);
#endif
/**
* @brief calculates the entrywise product of two Tensors
*/
......
......@@ -71,7 +71,7 @@ static misc::UnitTest tensor_solve("Tensor", "solve_Ax_equals_b", [](){
TEST((x2[{1,1}]) < 1e-14);
});
static misc::UnitTest tensor_solve_smallest_ev("Tensor", "get smallest eigenvalue of Matrix", [](){
static misc::UnitTest tensor_solve_smallest_ev("Tensor", "get smallest eigenvalue of Matrix (direct)", [](){
Index i,j,k,l,m,n,o,p;
Tensor A1 = Tensor::random({4,3,4,3});
......@@ -82,9 +82,41 @@ static misc::UnitTest tensor_solve_smallest_ev("Tensor", "get smallest eigenvalu
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);
TEST(frob_norm(A2(i,j,k,l,m,n,o,p)*x2(m,n,o,p) - lambda2 * x2(i,j,k,l)) < 5e-2);
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;
Tensor A1 = Tensor::random({4,3,4,3}) + (-1)*Tensor::identity({4,3,4,3});
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
get_smallest_eigenvalue_iterative(x1,A1,ev1.get(), 0, 10000, 1e-8);
double lambda1 = ev1[0];
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)));
Tensor A2 = Tensor::random({2,2,2,2,2,2,2,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];
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)));
Tensor A3 = Tensor::random({2,2,2,2,2,2,2,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];
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)));
});
//#endif
static misc::UnitTest solve_vs_lsqr("Tensor", "solve vs least squares", [](){
const size_t N = 500;
......
......@@ -22,6 +22,7 @@
* @brief Implementation of the blas and lapack wrapper functions.
*/
#ifdef ARPACK_LIBRARIES
#include <complex.h>
// fix for non standard-conform complex implementation
......@@ -61,22 +62,18 @@
namespace xerus {
namespace arpackWrapper {
/// Solves Ax = x*lambda for x and lambda
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) {
/// Solves Ax = x*lambda for x and lambda for the _k smallest eigenvalues
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) {
REQUIRE(_n <= static_cast<size_t>(std::numeric_limits<int>::max()), "Dimension to large for ARPACK");
REQUIRE(_info >= 0, "Info == 0, random; Info > 0, take residual; info is " << _info);
REQUIRE(_k < _n, "For some reason the number of Eigenvalues must be smaller than the dimension. see https://github.com/opencollab/arpack-ng/blob/master/SRC/dsaupd.f and error code -3" );
//REQUIRE(is_symmetric(_A, _n), "A must be symmetric");
// iteration variables dsaupd
int ido = 0;
auto bmat_option = arpack::bmat::identity;
auto ritz_option = arpack::which::smallest_algebraic;
int nev = static_cast<int>(_k);
int ncv = nev * 3; // TODO check for best value here
int ncv = nev * 3 >= static_cast<int>(_n) ? static_cast<int>(_n) : nev * 3; // TODO check for best value here
std::unique_ptr<double[]> v(new double[_n * ncv]);
std::unique_ptr<int[]> iparam(new int[11]);
std::unique_ptr<int[]> ipntr(new int[11]);
......@@ -84,13 +81,13 @@ namespace xerus {
int lworkl = ncv*(ncv + 8);
std::unique_ptr<double[]> workl(new double[lworkl]);
int info = 0;
// ev extraction variables dseupd
bool rvec = true;
auto howmny_option = arpack::howmny::ritz_vectors;
std::unique_ptr<int[]> select(new int[ncv]);
double sigma = 0.0;
int info1 = _info;
// intialization of iparam, parameters for arpack
iparam[0] = 1;
......@@ -98,12 +95,12 @@ namespace xerus {
iparam[6] = 1;
for (size_t i = 0; i < _maxiter; ++i){
XERUS_LOG(info, "Iter = " << i);
//LOG(info, "iter = " << i);
arpack::saupd( //NOTE for more details see https://github.com/opencollab/arpack-ng/blob/master/SRC/dsaupd.f
ido, // Reverse Communication flag
bmat_option, // Standard or Generalized EV problem
static_cast<int>(_n), // Dimension of the EV problem
ritz_option, /* Specify which of the Ritz values of OP to compute.
_ritz_option, /* Specify which of the Ritz values of OP to compute.
'LA' - compute the NEV largest (algebraic) eigenvalues.
'SA' - compute the NEV smallest (algebraic) eigenvalues.
'LM' - compute the NEV largest (in magnitude) eigenvalues.
......@@ -122,32 +119,20 @@ namespace xerus {
workd.get(), // Distributed array to be used in the basic Arnoldi iteration for reverse communication.
workl.get(), // rivate (replicated) array on each PE or array allocated on the front end.
lworkl, // LWORKL must be at least NCV**2 + 8*NCV .
info /* If INFO .EQ. 0, a randomly initial residual vector is used.
info1 /* If INFO .EQ. 0, a randomly initial residual vector is used.
If INFO .NE. 0, RESID contains the initial residual vector,
possibly from a previous run.
Error flag on output. */
);
XERUS_LOG(info, "saupd done, ido = " << ido);
if (ido == -1 or ido == 1){
XERUS_LOG(info, "ipntr = " << ipntr[0]);
XERUS_LOG(info, "ipntr = " << ipntr[1]);
XERUS_LOG(info, "A = " << _A[0] << " " << _A[1] << " " << _A[2] << " " << _A[3] << " " << _A[4] << " " << _A[5] << " " << _A[6] << " " << _A[7] << " " << _A[8] << " " << _A[9] << " " << _A[10] << " " << _A[11] << " " << _A[12] << " " << _A[13] << " " << _A[14] << " " << _A[15] << " ");
auto x = workd.get() + ipntr[0] - 1;
auto y = workd.get() + ipntr[1] - 1;
XERUS_LOG(info, "x = " << x[0] << " " << x[1] << " " << x[2] << " " << x[3] );
XERUS_LOG(info, "y = " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] );
blasWrapper::matrix_vector_product(workd.get() + ipntr[0] - 1,_n,1.0, _A, _n, false, workd.get() + ipntr[1] - 1);
XERUS_LOG(info, "y = " << y[0] << " " << y[1] << " " << y[2] << " " << y[3] );
blasWrapper::matrix_vector_product(workd.get() + ipntr[1] - 1,_n,1.0, _A, _n, false, workd.get() + ipntr[0] - 1);
}
else {
break;
}
}
XERUS_LOG(info, "saupd done");
REQUIRE(info1 >= 0, "ARPACK exited with error, see https://github.com/opencollab/arpack-ng/blob/master/SRC/dsaupd.f, error code is " << info1);
arpack::seupd( // NOTE for more details see https://github.com/opencollab/arpack-ng/blob/master/SRC/dseupd.f
rvec, // Specifies whether Ritz vectors corresponding to the Ritz value approximations to the eigenproblem A*z = lambda*B*z are computed.
howmny_option, /* Specifies how many Ritz vectors are wanted and the form of Z
......@@ -165,7 +150,7 @@ namespace xerus {
// Same as above
bmat_option,
static_cast<int>(_n),
ritz_option,
_ritz_option,
nev,
static_cast<double>(_eps),
_resid,
......@@ -177,19 +162,25 @@ namespace xerus {
workd.get(),
workl.get(),
lworkl,
info
info1
);
XERUS_LOG(info, "seupd done");
XERUS_LOG(info, "ev = " << _ev[0]);
XERUS_LOG(info, "iparam = " << iparam[0]);
XERUS_LOG(info, "ipntr = " << ipntr[0]);
return;
}
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) {
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) {
solve_ev (_x, _A, _ev, _k, _n, _resid, _maxiter, _eps, arpack::which::largest_algebraic, _info);
}
} // namespace arpackWrapper
} // namespace xerus
#endif
......@@ -118,9 +118,6 @@ namespace xerus {
REQUIRE(_m <= static_cast<size_t>(std::numeric_limits<int>::max()), "Dimension to large for BLAS/Lapack");
REQUIRE(_n <= static_cast<size_t>(std::numeric_limits<int>::max()), "Dimension to large for BLAS/Lapack");
XERUS_LOG(info, "y = " << _y[0] << " " << _y[1] << " " << _y[2] << " " << _y[3] );
XERUS_LOG(info, "x = " << _x[0] << " " << _x[1] << " " << _x[2] << " " << _x[3] );
XERUS_LOG(info, "A = " << _A[0] << " " << _A[1] << " " << _A[2] << " " << _A[3] << " " << _A[4] << " " << _A[5] << " " << _A[6] << " " << _A[7] << " " << _A[8] << " " << _A[9] << " " << _A[10] << " " << _A[11] << " " << _A[12] << " " << _A[13] << " " << _A[14] << " " << _A[15] << " ");
XERUS_PA_START;
if(!_transposed) {
......@@ -128,8 +125,6 @@ namespace xerus {
} else {
cblas_dgemv(CblasRowMajor, CblasTrans, static_cast<int>(_n), static_cast<int>(_m), _alpha, _A, static_cast<int>(_m) , _y, 1, 0.0, _x, 1);
}
XERUS_LOG(info, "y = " << _y[0] << " " << _y[1] << " " << _y[2] << " " << _y[3] );
XERUS_PA_END("Dense BLAS", "Matrix Vector Product", misc::to_string(_m)+"x"+misc::to_string(_n)+" * "+misc::to_string(_n));
}
......
......@@ -30,6 +30,7 @@
#include <xerus/misc/math.h>
#include <xerus/misc/internal.h>
#include <xerus/arpackWrapper.h>
#include <xerus/blasLapackWrapper.h>
#include <xerus/cholmod_wrapper.h>
#include <xerus/sparseTimesFullContraction.h>
......@@ -1758,7 +1759,53 @@ namespace xerus {
return ev;
}
#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) {
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;
// Calculate multDimensions
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, _A.dimensions.begin() + degN))
{
Tensor::DimensionTuple newDimX;
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
}
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);
// Note that A is dense here
arpackWrapper::solve_ev_smallest(
rev.get(), // right ritz vectors
_A.get_unsanitized_dense_data(),
_ev, 1, n,
res.get(),
_miter,
_eps, info
);
// eigenvector of smallest eigenvalue
auto tmpX = _X.override_dense_data();
for (size_t i = 0; i < n; ++i)
tmpX[i] = rev[i];
return;
}
#endif
Tensor entrywise_product(const Tensor &_A, const Tensor &_B) {
REQUIRE(_A.dimensions == _B.dimensions, "Entrywise product ill-defined for non-equal dimensions.");
if(_A.is_dense() && _B.is_dense()) {
......
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