Commit 564f04cf authored by Philipp  Trunschke's avatar Philipp Trunschke

Merge branch 'development' of git.hemio.de:xerus/xerus into development

parents be157b50 df47de0f
Pipeline #1052 passed with stages
in 7 minutes and 40 seconds
......@@ -151,7 +151,7 @@ build/libxerus_misc.so: $(MINIMAL_DEPS) $(MISC_SOURCES)
build/libxerus.so: $(MINIMAL_DEPS) $(XERUS_SOURCES) build/libxerus_misc.so
mkdir -p $(dir $@)
$(CXX) -shared -fPIC -Wl,-soname,libxerus.so $(FLAGS) -I include $(XERUS_SOURCES) -L ./build/ -Wl,--as-needed -lxerus_misc $(SUITESPARSE) $(LAPACK_LIBRARIES) $(BLAS_LIBRARIES) -o build/libxerus.so
$(CXX) -shared -fPIC -Wl,-soname,libxerus.so $(FLAGS) -I include $(XERUS_SOURCES) -L ./build/ -Wl,--as-needed -lxerus_misc $(SUITESPARSE) $(LAPACK_LIBRARIES) $(ARPACK_LIBRARIES) $(BLAS_LIBRARIES) -o build/libxerus.so
python2: build/python2/xerus.so
......@@ -221,7 +221,7 @@ endif
$(TEST_NAME): $(MINIMAL_DEPS) $(UNIT_TEST_OBJECTS) $(TEST_OBJECTS) build/libxerus.a build/libxerus_misc.a
$(CXX) -D XERUS_UNITTEST $(FLAGS) $(UNIT_TEST_OBJECTS) $(TEST_OBJECTS) build/libxerus.a build/libxerus_misc.a $(SUITESPARSE) $(LAPACK_LIBRARIES) $(BLAS_LIBRARIES) $(CALLSTACK_LIBS) -o $(TEST_NAME)
$(CXX) -D XERUS_UNITTEST $(FLAGS) $(UNIT_TEST_OBJECTS) $(TEST_OBJECTS) build/libxerus.a build/libxerus_misc.a $(SUITESPARSE) $(LAPACK_LIBRARIES) $(ARPACK_LIBRARIES) $(BLAS_LIBRARIES) $(CALLSTACK_LIBS) -o $(TEST_NAME)
build/print_boost_version: src/print_boost_version.cpp
@$(CXX) -o $@ $<
......@@ -258,7 +258,7 @@ clean:
benchmark: $(MINIMAL_DEPS) $(LOCAL_HEADERS) benchmark.cxx $(LIB_NAME_STATIC)
$(CXX) $(FLAGS) benchmark.cxx $(LIB_NAME_STATIC) $(SUITESPARSE) $(LAPACK_LIBRARIES) $(BLAS_LIBRARIES) $(CALLSTACK_LIBS) -lboost_filesystem -lboost_system -o Benchmark
$(CXX) $(FLAGS) benchmark.cxx $(LIB_NAME_STATIC) $(SUITESPARSE) $(LAPACK_LIBRARIES) $(ARPACK_LIBRARIES) $(BLAS_LIBRARIES) $(CALLSTACK_LIBS) -lboost_filesystem -lboost_system -o Benchmark
# Build rule for normal misc objects
build/.miscObjects/%.o: %.cpp $(MINIMAL_DEPS)
......@@ -299,7 +299,7 @@ endif
# Build and execution rules for tutorials
build/.tutorialObjects/%: %.cpp $(MINIMAL_DEPS) build/libxerus.a build/libxerus_misc.a
mkdir -p $(dir $@)
$(CXX) -I include $< build/libxerus.a build/libxerus_misc.a $(SUITESPARSE) $(LAPACK_LIBRARIES) $(BLAS_LIBRARIES) $(CALLSTACK_LIBS) $(FLAGS) -MMD -o $@
$(CXX) -I include $< build/libxerus.a build/libxerus_misc.a $(SUITESPARSE) $(LAPACK_LIBRARIES) $(ARPACK_LIBRARIES) $(BLAS_LIBRARIES) $(CALLSTACK_LIBS) $(FLAGS) -MMD -o $@
# Build rule for the preCompileHeader
......
......@@ -134,6 +134,10 @@ 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
# OTHER += -DARPACK_LIBRARIES
# Uncomment or add the appropriate boost python library
BOOST_PYTHON2 = -lboost_python-py27
BOOST_PYTHON3 = -lboost_python-py35
......
......@@ -36,6 +36,9 @@
// All the xerus headers
#include "xerus/blasLapackWrapper.h"
#ifdef ARPACK_LIBRARIES
#include "xerus/arpackWrapper.h"
#endif
#include "xerus/index.h"
#include "xerus/indexedTensorReadOnly.h"
#include "xerus/indexedTensorWritable.h"
......
// Xerus - A General Purpose Tensor Library
// Copyright (C) 2014-2018 Benjamin Huber and Sebastian Wolf.
//
// Xerus is free software: you can redistribute it and/or modify
// it under the terms of the GNU Affero General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// Xerus is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Affero General Public License for more details.
//
// You should have received a copy of the GNU Affero General Public License
// along with Xerus. If not, see <http://www.gnu.org/licenses/>.
//
// For further information on Xerus visit https://libXerus.org
// or contact us at contact@libXerus.org.
/**
* @file
* @brief Header file for the arpack wrapper functions.
*/
#pragma once
#ifdef ARPACK_LIBRARIES
#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>
namespace xerus {
class Tensor;
/**
* @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
* a seriously low level implementation of a critical part of an algorithm is required.
*/
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, 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);
//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);
///@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);
}
}
#endif
......@@ -1023,7 +1023,22 @@ 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);
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);
#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;
......
This diff is collapsed.
......@@ -118,14 +118,13 @@ 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_PA_START;
if(!_transposed) {
cblas_dgemv(CblasRowMajor, CblasNoTrans, static_cast<int>(_m), static_cast<int>(_n), _alpha, _A, static_cast<int>(_n), _y, 1, 0.0, _x, 1);
} 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_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,99 @@ 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;
}
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
//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);
const size_t degN = 4;
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];
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])
{
_X.reset({_l.dimensions[2], _A.dimensions[2], _A1.dimensions[2], _r.dimensions[2]}, 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_dmrg_special(
rev.get(), // right ritz vectors
_l,_A,_A1,_r,
_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()) {
......
......@@ -616,7 +616,7 @@ namespace xerus {
}
}
while (exceeds_maximal_ranks()) {
while (exceeds_maximal_ranks() && !_keepRank) {
// Move left from given CorePosition
for (size_t n = _position; n > 0; --n) {
transfer_core(n+1, n, !_keepRank);
......@@ -1147,7 +1147,6 @@ namespace xerus {
}
}
}
// Use Tensor fallback
if (_other.tensorObjectReadOnly->nodes.size() > 1) {
LOG_ONCE(warning, "Assigning a general tensor network to TTOperator not yet implemented. casting to fullTensor first");
......
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