Commit 5bb317e6 authored by Michael Goette's avatar Michael Goette

preliminary results on EV solver integration, wrapper is ready

parent a1387a6a
Pipeline #942 passed with stages
in 10 minutes and 1 second
......@@ -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
......
......@@ -650,6 +650,42 @@ 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,41 @@ 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);
}
// // Note that A is dense here
// blasWrapper::solve_ev(
// _X.override_dense_data(),
// _A.get_unsanitized_dense_data(), m, n,
// _B.get_unsanitized_dense_data(), p);
//
// // Propagate the constant factor
// _X.factor = _B.factor / _A.factor;
return 0.0;
}
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