Commit 9df6368f authored by Philipp  Trunschke's avatar Philipp Trunschke

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

parents edc0d5d5 b8fa7ccb
Pipeline #947 passed with stages
in 9 minutes and 55 seconds
......@@ -131,7 +131,8 @@ namespace xerus {
///@brief: Solves Ax = b for x
void solve(double* _x, const double* _A, size_t _m, size_t _n, const double* _b, size_t _p);
///@brief: Solves Ax = lambda*x for x, this calls the Lapack Routine DGEEV
void solve_ev(double* const _x, double* const _re, double* const _im, const double* const _A, const size_t _n);
///@brief: Solves min ||Ax - b||_2 for x
void solve_least_squares( double* const _x, const double* const _A, const size_t _m, const size_t _n, const double* const _b, const size_t _p);
......
......@@ -1014,6 +1014,15 @@ namespace xerus {
* @param _extraDegree number of modes that @a _x and @a _B sharefor which the solution should be computed independently.
*/
void solve(Tensor &_X, const Tensor &_A, const Tensor &_B, size_t _extraDegree = 0);
/**
* @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 _A input Operator A symmetric with respect to matrification
* @return the smallest eigenvalue
*/
double get_smallest_eigenvalue(Tensor &_X, const Tensor &_A);
/**
* @brief calculates the entrywise product of two Tensors
......
......@@ -32,7 +32,7 @@ static misc::UnitTest tensor_solve("Tensor", "solve_Ax_equals_b", [](){
A1[{1,0,1}] = 1;
A1[{2,1,0}] = 1;
A1[{3,1,1}] = 1;
Tensor b1({4});
b1[0] = 73;
b1[1] = -73;
......@@ -71,6 +71,20 @@ 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", [](){
Index i,j,k,l,m,n,o,p;
Tensor A1 = Tensor::random({4,3,4,3});
Tensor x1({4,3});
double lambda1 = get_smallest_eigenvalue(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);
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);
});
static misc::UnitTest solve_vs_lsqr("Tensor", "solve vs least squares", [](){
const size_t N = 500;
......
......@@ -650,6 +650,40 @@ namespace xerus {
XERUS_PA_END("Dense LAPACK", "Solve (LDL)", misc::to_string(_n)+"x"+misc::to_string(_n)+"x"+misc::to_string(_nrhs));
}
/// Solves Ax = x*lambda for x and lambda
void solve_ev(double* const _x, double* const _re, double* const _im, const double* const _A, const size_t _n) {
REQUIRE(_n <= static_cast<size_t>(std::numeric_limits<int>::max()), "Dimension to large for BLAS/Lapack");
const std::unique_ptr<double[]> tmpA(new double[_n*_n]);
misc::copy(tmpA.get(), _A, _n * _n);
LOG(debug, "solving with...");
//so far only non symmetric -> dgeev
LOG(debug, "DGEEV");
XERUS_PA_START;
std::unique_ptr<double[]> leftev(new double[1]);
IF_CHECK( int lapackAnswer = ) LAPACKE_dgeev(
LAPACK_ROW_MAJOR,
'N', // No left eigenvalues are computed
'V', // Right eigenvalues are computed
static_cast<int>(_n), // Dimensions of A (nxn)
tmpA.get(), // input: A, output: L and U
static_cast<int>(_n), // LDA
_re, // real part of the eigenvalues
_im, // imaginary part of the eigenvalues
leftev.get(), // output: left eigenvectors, here dummy
static_cast<int>(_n), // LDVL
_x, // right eigenvectors
static_cast<int>(_n) // LDVR TODO check size of _x
);
CHECK(lapackAnswer == 0, error, "Unable to solve Ax = lambda*x (DGEEV solver). Lapacke says: " << lapackAnswer);
XERUS_PA_END("Dense LAPACK", "Solve (DGEEV)", misc::to_string(_n)+"x"+misc::to_string(_n)+"x"+misc::to_string(_nrhs));
return;
}
void solve_least_squares( double* const _x, const double* const _A, const size_t _m, const size_t _n, const double* const _b, const size_t _p){
const std::unique_ptr<double[]> tmpA(new double[_m*_n]);
......
......@@ -1708,6 +1708,55 @@ namespace xerus {
_X.factor = _B.factor / _A.factor;
}
double get_smallest_eigenvalue(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");
const size_t degN = _A.degree() / 2;
// Calculate multDimensions
const size_t m = misc::product(_A.dimensions, 0, degN);
const size_t n = misc::product(_A.dimensions, degN, 2*degN);
REQUIRE(m == n, "Inconsistent dimensions.");
// 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);
}
std::unique_ptr<double[]> re(new double[n]); // real eigenvalues
std::unique_ptr<double[]> im(new double[n]); // imaginary eigenvalues
std::unique_ptr<double[]> rev(new double[n * n]); // right eigenvectors
// Note that A is dense here
blasWrapper::solve_ev(
rev.get(), // right eigenvectors
re.get(), // real eigenvalues
im.get(), // imaginary eigenvalues
_A.get_unsanitized_dense_data(), n);
//get eigenvalue with minimal real part
double ev = re[0];
size_t idx = 0;
for (size_t i = 1; i < n; ++i) {
if (re[i] < ev){
ev = re[i];
idx = i;
}
}
//get eigenvector oof smallest eigenvalue
auto tmpX = _X.override_dense_data();
for (size_t i = 0; i <= n; ++i)
tmpX[i] = rev[idx + i * n];
return ev;
}
Tensor entrywise_product(const Tensor &_A, const Tensor &_B) {
......
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