Commit 180b35d7 authored by Sebastian Wolf's avatar Sebastian Wolf

Merge branch 'development' of git:xerus/xerus into development

parents f56ae533 3ad08c42
Pipeline #1114 passed with stages
in 22 minutes and 57 seconds
......@@ -90,7 +90,7 @@ static misc::UnitTest tensor_svd_rnd512("Tensor", "SVD_Random_512x512", [](){
(res1(i,j,k,o), res2(o,p), res3(p,l,m,n)) = SVD(A(l,i,m,k,j,n));
res4(k,i,m,l,j,n) = res1(i,j,l,o)*res2(o,p)*res3(p,k,m,n);
TEST(approx_equal(res4, A, 1e-14));
MTEST(approx_equal(res4, A, 1e-14), frob_norm(res4-A));
MTEST(frob_norm(res1(i^3, m)*res1(i^3, n) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " U not orthogonal");
MTEST(frob_norm(res3(m, i^3)*res3(n, i^3) - Tensor::identity(res2.dimensions)(m, n)) < 1e-12, " Vt not orthogonal");
......
......@@ -65,11 +65,14 @@ extern "C"
namespace xerus {
namespace blasWrapper {
// the following routines use a work array as temporary storage
// to avoid the overhead of repeated reallocation for many small calls, every thread pre-allocates a small scratch-space
// const size_t DEFAULT_WORKSPACE_SIZE = 1024*1024;
// thread_local value_t defaultWorkspace[DEFAULT_WORKSPACE_SIZE]; // NOTE recheck compatibility with eigen (dolfin) when reinserting this!
/// @brief stores in @a _out the transpose of the @a _m x @a _n matrix @a _in
void low_level_transpose(double * _out, double * _in, size_t _m, size_t _n) {
for (size_t i=0; i<_m; ++i) {
for (size_t j=0; j<_n; ++j) {
_out[j*_m + i] = _in[i*_n+j];
}
}
}
//----------------------------------------------- LEVEL I BLAS ----------------------------------------------------------
......@@ -206,26 +209,42 @@ namespace xerus {
}
lapack_int dgesdd_get_workarray_size(lapack_int m, lapack_int n) {
lapack_int info = 0;
char job = 'S';
double work = 0;
lapack_int lwork = -1;
lapack_int min = std::min(m,n);
LAPACK_dgesdd( &job, &m, &n, nullptr, &m, nullptr, nullptr, &m, nullptr, &min, &work, &lwork, nullptr, &info );
REQUIRE(info == 0, "work array size query of dgesdd returned " << info);
return lapack_int(work);
}
void dgesdd_work(lapack_int m, lapack_int n, double* a, double* s, double* u, double* vt, double* work, lapack_int lwork, lapack_int* iwork ) {
REQUIRE(lwork > 0, "");
lapack_int info = 0;
char job = 'S';
lapack_int min = std::min(m,n);
// if A = U*S*V^T, then A^T = V^T*S*U^T, so we can simply call lapack without transposing our matrices
// by simply changing the order of U and Vt
LAPACK_dgesdd( &job, &n, &m, a, &n, s, vt, &n, u, &min, work, &lwork, iwork, &info );
REQUIRE(info == 0, "dgesdd failed with info " << info);
}
void svd_destructive( double* const _U, double* const _S, double* const _Vt, double* const _A, const size_t _m, const size_t _n) {
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;
std::unique_ptr<double[]> tmpA(new double[_m*_n]);
misc::copy(tmpA.get(), _A, _m*_n);
lapack_int m = lapack_int(_m);
lapack_int n = lapack_int(_n);
std::unique_ptr<lapack_int[]> iwork(new lapack_int[std::max(1ul,8u*std::min(_m,_n))]);
lapack_int lwork = dgesdd_get_workarray_size(m, n);
std::unique_ptr<double[]> work(new double[size_t(lwork)]);
int lapackAnswer = LAPACKE_dgesdd(LAPACK_ROW_MAJOR, 'S', static_cast<int>(_m), static_cast<int>(_n), _A, static_cast<int>(_n), _S, _U, static_cast<int>(std::min(_m, _n)), _Vt, static_cast<int>(_n));
CHECK(lapackAnswer == 0, warning, "Lapack failed to compute SVD. Answer is: " << lapackAnswer);
CHECK(lapackAnswer == 0, warning, "Call was: LAPACKE_dgesdd(LAPACK_ROW_MAJOR, 'S', " << static_cast<int>(_m) << ", " << static_cast<int>(_n) << ", " << _A << ", " << static_cast<int>(_n) <<", "
<< _S <<", " << _U << ", " << static_cast<int>(std::min(_m, _n)) << ", " << _Vt << ", " << static_cast<int>(_n) << ");");
if(lapackAnswer != 0) {
LOG(warning, "SVD failed ");
// for(size_t i=0; i < _m; ++i) {
// for(size_t j=0; j < _n; ++j) {
// LOG(warning, tmpA[i*_n+j]);
// }
// }
}
dgesdd_work( m, n, _A, _S, _U, _Vt, work.get(), lwork, iwork.get());
XERUS_PA_END("Dense LAPACK", "Singular Value Decomposition", misc::to_string(_m)+"x"+misc::to_string(_n));
}
......@@ -679,7 +698,7 @@ namespace xerus {
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));
XERUS_PA_END("Dense LAPACK", "Solve (DGEEV)", misc::to_string(_n)+"x"+misc::to_string(_n));
return;
}
......
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