Commit 3fff1401 authored by Michael Goette's avatar Michael Goette

added abitrary tensor networks as operator for special eigfenvalue iterations

parent df47de0f
......@@ -49,7 +49,7 @@
namespace xerus {
class Tensor;
class TensorNetwork;
/**
* @brief In this namespace the minimal wrappers for the ARPACK functions are collected.
* @details As an end user of xerus it should never be nessecary to call any of these functions, unless
......@@ -66,9 +66,9 @@ namespace xerus {
//TODO check if this can be simplified!!
///@brief: Solves Ax = lambda*x for x, this calls the Arpack Routine dsaupd
void solve_ev_dmrg_special(double* const _x, const Tensor& _l, const Tensor& _A, const Tensor& _A1, const Tensor& _r, 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);
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
void solve_ev_smallest_dmrg_special(double* const _x, const Tensor& _l, const Tensor& _A, const Tensor& _A1, const Tensor& _r, 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_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);
......
......@@ -1035,7 +1035,7 @@ namespace xerus {
* @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);
void get_smallest_eigenvalue_iterative_dmrg_special(Tensor& _X, const Tensor& _l, const Tensor& _A, const Tensor& _A1, const Tensor& _r, 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);
#endif
......
......@@ -59,6 +59,7 @@
#include <xerus/misc/basicArraySupport.h>
#include <xerus/misc/math.h>
#include <xerus/misc/internal.h>
#include "xerus/tensorNetwork.h"
#include "xerus/tensor.h"
......@@ -173,13 +174,14 @@ namespace xerus {
}
/// Solves Ax = x*lambda for x and lambda for the _k smallest eigenvalues
void solve_ev_dmrg_special(double* const _x, const Tensor& _l, const Tensor& _A, const Tensor& _A1, const Tensor& _r, 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) {
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) {
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
std::vector<size_t> dims(_op.dimensions.begin(), _op.dimensions.begin() +_op.dimensions.size()/2);
int ido = 0;
auto bmat_option = arpack::bmat::identity;
int nev = static_cast<int>(_k);
......@@ -236,12 +238,12 @@ namespace xerus {
);
if (ido == -1 or ido == 1){
Tensor tmpX({_l.dimensions[2], _A.dimensions[2], _A1.dimensions[2], _r.dimensions[2]});
Tensor tmpX1({_l.dimensions[0], _A.dimensions[1], _A1.dimensions[1], _r.dimensions[0]});
Tensor tmpX(dims);
Tensor tmpX1(dims);
auto tmpX_ptr = tmpX.override_dense_data();
misc::copy(tmpX_ptr, workd.get() + ipntr[0] - 1, _n);
Index i1,i2,i3,i4,j1,j2,j3,j4,k1,k2,k3;
tmpX1(i1,i2,i3,i4) = _l(i1, k1, j1)* (_A(k1, i2, j2, k2) * (_A1(k2,i3,j3,k3)* (_r(i4, k3, j4) * tmpX(j1,j2,j3,j4))));
tmpX1(i1&0) = _op(i1/2,j1/2) * tmpX(j1&0);
auto tmpX1_ptr = tmpX1.override_dense_data();
misc::copy(workd.get() + ipntr[1] - 1,tmpX1_ptr, _n);
......@@ -294,8 +296,8 @@ namespace xerus {
solve_ev (_x, _A, _ev, _k, _n, _resid, _maxiter, _eps, arpack::which::largest_algebraic, _info);
}
void solve_ev_smallest_dmrg_special(double* const _x, const Tensor& _l, const Tensor& _A, const Tensor& _A1, const Tensor& _r, 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_dmrg_special (_x, _l, _A, _A1, _r, _ev, _k, _n, _resid, _maxiter, _eps, arpack::which::smallest_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);
}
} // namespace arpackWrapper
......
......@@ -1805,27 +1805,24 @@ namespace xerus {
return;
}
void get_smallest_eigenvalue_iterative_dmrg_special(Tensor& _X, const Tensor& _l, const Tensor& _A, const Tensor& _A1, const Tensor& _r, double* const _ev, int _info, const size_t _miter, const double _eps) {
REQUIRE(_l.is_dense() && _A.is_dense() && _A1.is_dense() && _r.is_dense(), "for now only dense is implemented"); //TODO implement sparse
void get_smallest_eigenvalue_iterative(Tensor& _X, const TensorNetwork& _op, double* const _ev, int _info, const size_t _miter, const double _eps) {
//REQUIRE(&_X != &_A, "Not supportet yet");
REQUIRE(_l.degree() == 3, "The tensor _l needs to be of order 3");
REQUIRE(_A.degree() == 4, "The tensor _A needs to be of order 4");
REQUIRE(_A1.degree() == 4, "The tensor _A1 needs to be of order 4");
REQUIRE(_r.degree() == 3, "The tensor _r needs to be of order 3");
REQUIRE(_eps > 0 && _eps < 1, "epsilon must be betweeen 0 and 1, given " << _eps);
REQUIRE(_op.degree() % 2 == 0, "operator degree must be positive");
const size_t degN = 4;
const size_t degN = _op.degree() / 2;
int info = _info;
// Calculate multDimensions
const size_t m = _l.dimensions[0] * _A.dimensions[1] * _A1.dimensions[1] * _r.dimensions[0];
const size_t n = _l.dimensions[2] * _A.dimensions[2] * _A1.dimensions[2] * _r.dimensions[2];
const size_t m = misc::product(_op.dimensions, 0, degN);
const size_t n = misc::product(_op.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 || _X.dimensions[0] !=_l.dimensions[2] || _X.dimensions[1] !=_A.dimensions[2] || _X.dimensions[2] !=_A1.dimensions[2] || _X.dimensions[3] !=_r.dimensions[2])
if( _X.degree() != degN || !std::equal(_X.dimensions.begin(), _X.dimensions.begin() + degN, _op.dimensions.begin() + degN))
{
_X.reset({_l.dimensions[2], _A.dimensions[2], _A1.dimensions[2], _r.dimensions[2]}, Tensor::Representation::Dense, Tensor::Initialisation::None);
Tensor::DimensionTuple newDimX;
newDimX.insert(newDimX.end(), _op.dimensions.begin()+degN, _op.dimensions.end());
_X.reset(std::move(newDimX), Tensor::Representation::Dense, Tensor::Initialisation::None);
info = 0; // info must be 0 if X is not random
}
......@@ -1835,9 +1832,9 @@ namespace xerus {
misc::copy(res.get(), _X.get_unsanitized_dense_data(), n);
// Note that A is dense here
arpackWrapper::solve_ev_smallest_dmrg_special(
arpackWrapper::solve_ev_smallest_special(
rev.get(), // right ritz vectors
_l,_A,_A1,_r,
_op,
_ev, 1, n,
res.get(),
_miter,
......
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