Commit b1aba572 authored by Benjamin Huber's avatar Benjamin Huber

using sparse QR to solve sparse systems instead of umfpack

parent 5b0b2ffd
Pipeline #173 passed with stages
......@@ -18,7 +18,7 @@ BLAS_LIBRARIES = -lopenblas -lgfortran # Openblas, serial
LAPACK_LIBRARIES = -llapacke -llapack # Standard Lapack + Lapacke libraries
# Uncomment or add the appropriate CXSparse library
SUITESPARSE = -lcholmod -lumfpack -lspqr
SUITESPARSE = -lcholmod -lspqr
# custom include paths
......
......@@ -12,7 +12,7 @@ LOGGING += -D INFO_ # Information that is not linked to any
#=================================================================================================
BLAS_LIBRARIES = -lopenblas -lgfortran # Openblas, serial
LAPACK_LIBRARIES = -llapacke -llapack # Standard Lapack + Lapacke libraries
SUITESPARSE = -lcholmod -lumfpack -lspqr
SUITESPARSE = -lcholmod -lspqr
# custom include paths
......
......@@ -17,7 +17,7 @@ LOGGING += -D INFO_ # Information that is not linked to any
#=================================================================================================
BLAS_LIBRARIES = -lopenblas -lgfortran # Openblas, serial
LAPACK_LIBRARIES = -llapacke -llapack # Standard Lapack + Lapacke libraries
SUITESPARSE = -lcholmod -lumfpack -lspqr
SUITESPARSE = -lcholmod -lspqr
# custom include paths
......
......@@ -120,7 +120,7 @@ LAPACK_LIBRARIES = -llapacke -llapack # Standard Lapack + Lapacke librarie
# LAPACK_LIBRARIES = ../lib/lapack/liblapacke.a ../lib/lapack/liblapack.a # Custom
# Uncomment or add the appropriate CXSparse library
SUITESPARSE = -lcholmod -lumfpack -lspqr
SUITESPARSE = -lcholmod -lspqr
# Custom include paths
......
......@@ -142,5 +142,5 @@ UNIT_TEST(Tensor, solve_transposed,
x3(i) = r(j) / A(i,j);
x4(i) = r(j) / At(j,i);
MTEST(frob_norm(x3-x4) < 1e-14, "d " << frob_norm(x3-x4));
MTEST(frob_norm(x1-x3)/frob_norm(x1) < 3e-14, "sd " << frob_norm(x1-x3)/frob_norm(x1));
MTEST(frob_norm(x1-x3)/frob_norm(x1) < 5e-14, "sd " << frob_norm(x1-x3)/frob_norm(x1));
)
......@@ -28,8 +28,6 @@
#include <xerus/misc/performanceAnalysis.h>
#include <xerus/misc/basicArraySupport.h>
#include <suitesparse/umfpack.h>
#include <xerus/misc/check.h>
namespace xerus { namespace internal {
......@@ -176,9 +174,13 @@ namespace xerus { namespace internal {
const std::map<size_t, double>& _b,
size_t _bDim)
{
LOG(fatal, "not yet implemented");
const CholmodSparse A(_A, _transposeA?_xDim:_bDim, _transposeA?_bDim:_xDim, _transposeA);
const CholmodSparse b(_b, 1, _bDim, true); // avoids transpose call by giving transposed dimensions
CholmodSparse x(SuiteSparseQR<double>(0, xerus::EPSILON, A.matrix.get(), b.matrix.get(), cholmodObject.get()));
_x = x.to_map();
}
void CholmodSparse::solve_dense_rhs(double * _x,
size_t _xDim,
const std::map<size_t, double>& _A,
......@@ -186,17 +188,24 @@ namespace xerus { namespace internal {
const double* _b,
size_t _bDim)
{
REQUIRE(_xDim == _bDim, "solving sparse systems only implemented for square matrices so far");
REQUIRE(_xDim < std::numeric_limits<long>::max() && _bDim < std::numeric_limits<long>::max() && _A.size() < std::numeric_limits<long>::max(),
"sparse matrix given to qc decomposition too large for suitesparse"
);
const CholmodSparse A(_A, _transposeA?_xDim:_bDim, _transposeA?_bDim:_xDim, _transposeA);
void *symbolic, *numeric;
// double control[UMFPACK_CONTROL], info[UMFPACK_INFO];
// TODO check return values
umfpack_dl_symbolic(long(_bDim), long(_xDim), static_cast<long*>(A.matrix->p), static_cast<long*>(A.matrix->i), static_cast<double*>(A.matrix->x), &symbolic, nullptr, nullptr);
umfpack_dl_numeric(static_cast<long*>(A.matrix->p), static_cast<long*>(A.matrix->i), static_cast<double*>(A.matrix->x), symbolic, &numeric, nullptr, nullptr);
umfpack_dl_free_symbolic(&symbolic);
umfpack_dl_solve(UMFPACK_A, static_cast<long*>(A.matrix->p), static_cast<long*>(A.matrix->i), static_cast<double*>(A.matrix->x), _x, _b, numeric, nullptr, nullptr);
umfpack_dl_free_numeric(&numeric);
cholmod_dense b{
_bDim, 1, _bDim, _bDim,
static_cast<void*>(const_cast<double*>(_b)), nullptr,
CHOLMOD_REAL, CHOLMOD_DOUBLE
};
std::unique_ptr<cholmod_dense, std::function<void(cholmod_dense*)>> x(
SuiteSparseQR<double>(A.matrix.get(), &b, cholmodObject.get()),
[](cholmod_dense* _toDelete) {
cholmod_l_free_dense(&_toDelete, cholmodObject.get());
}
);
REQUIRE(x->z == nullptr, "IE");
REQUIRE(x->nrow == _xDim, "IE");
misc::copy(_x, static_cast<double*>(x->x), _xDim);
}
......@@ -207,15 +216,18 @@ namespace xerus { namespace internal {
size_t _n,
bool _fullrank)
{
REQUIRE(_m < std::numeric_limits<long>::max() && _n < std::numeric_limits<long>::max() && _A.size() < std::numeric_limits<long>::max(),
"sparse matrix given to qc decomposition too large for suitesparse"
);
CholmodSparse A(_A, _m, _n, _transposeA);
cholmod_sparse *Q, *R;
SuiteSparse_long *E;
size_t rank = SuiteSparseQR<double>(0, xerus::EPSILON, _fullrank?std::min(_m,_n):1, A.matrix.get(), &Q, &R, &E, cholmodObject.get());
long rank = SuiteSparseQR<double>(0, xerus::EPSILON, _fullrank?long(std::min(_m,_n)):1, A.matrix.get(), &Q, &R, &E, cholmodObject.get());
CholmodSparse Qs(Q);
CholmodSparse Rs(R);
REQUIRE(E == nullptr, "IE: sparse QR returned a permutation despite fixed ordering?!");
REQUIRE((_fullrank?std::min(_m,_n):rank) == Qs.matrix->ncol, "IE: strange rank deficiency after sparse qr " << (_fullrank?std::min(_m,_n):rank) << " vs " << Qs.matrix->ncol);
return std::make_tuple(Qs.to_map(), Rs.to_map(), _fullrank?std::min(_m,_n):rank);
REQUIRE((_fullrank?std::min(_m,_n):size_t(rank)) == Qs.matrix->ncol, "IE: strange rank deficiency after sparse qr " << (_fullrank?long(std::min(_m,_n)):rank) << " vs " << Qs.matrix->ncol);
return std::make_tuple(Qs.to_map(), Rs.to_map(), _fullrank?std::min(_m,_n):size_t(rank));
}
std::tuple<std::map<size_t, double>, std::map<size_t, double>, size_t> CholmodSparse::cq(
......@@ -225,19 +237,22 @@ namespace xerus { namespace internal {
size_t _n,
bool _fullrank)
{
REQUIRE(_m < std::numeric_limits<long>::max() && _n < std::numeric_limits<long>::max() && _A.size() < std::numeric_limits<long>::max(),
"sparse matrix given to qc decomposition too large for suitesparse"
);
CholmodSparse A(_A, _n, _m, !_transposeA);
cholmod_sparse *Q, *R;
SuiteSparse_long *E;
// decompose A^T = q^T*r^T
size_t rank = SuiteSparseQR<double>(0, xerus::EPSILON, _fullrank?std::min(_m,_n):1, A.matrix.get(), &Q, &R, &E, cholmodObject.get());
long rank = SuiteSparseQR<double>(0, xerus::EPSILON, _fullrank?long(std::min(_m,_n)):1, A.matrix.get(), &Q, &R, &E, cholmodObject.get());
CholmodSparse Qs(Q);
CholmodSparse Rs(R);
REQUIRE(E == nullptr, "IE: sparse QR returned a permutation despite fixed ordering?!");
REQUIRE((_fullrank?std::min(_m,_n):rank) == Qs.matrix->ncol, "IE: strange rank deficiency after sparse qr " << (_fullrank?std::min(_m,_n):rank) << " vs " << Qs.matrix->ncol);
REQUIRE((_fullrank?std::min(_m,_n):size_t(rank)) == Qs.matrix->ncol, "IE: strange rank deficiency after sparse qr " << (_fullrank?long(std::min(_m,_n)):rank) << " vs " << Qs.matrix->ncol);
//transpose q and r to get r*q=A
Qs.transpose();
Rs.transpose();
return std::make_tuple(Rs.to_map(), Qs.to_map(), _fullrank?std::min(_m,_n):rank);
return std::make_tuple(Rs.to_map(), Qs.to_map(), _fullrank?std::min(_m,_n):size_t(rank));
}
......
......@@ -85,10 +85,7 @@ namespace xerus {
// We need tmp objects for A and b, because Lapacke wants to destroys the input
// and we need to cast b if it is sparse as we cannot handle that yet
if (_b.tensorObjectReadOnly->is_sparse()) {
LOG_ONCE(info, "solve operator cannot make use of a sparse rhs yet. casting rhs to dense representation first from now on.");
}
if (_b.tensorObjectReadOnly->is_sparse() || !_a.tensorObjectReadOnly->is_sparse()) {
if (!_a.tensorObjectReadOnly->is_sparse()) {
saveSlotB.reset(new internal::IndexedTensor<Tensor>(
new Tensor(std::move(dimensionsB), Tensor::Representation::Dense, Tensor::Initialisation::None), orderB, true
));
......
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