From 984d010b7bcb6a03f0319e79b8a768587be85422 Mon Sep 17 00:00:00 2001 From: Ralf Hannemann-Tamas Date: Mon, 8 Feb 2021 22:00:31 +0000 Subject: add specialization of check_sparse_solving() for SuperLU solver, in order to test adjoint and transpose solves --- test/sparse_solver.h | 131 +++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 131 insertions(+) (limited to 'test/sparse_solver.h') diff --git a/test/sparse_solver.h b/test/sparse_solver.h index f45d7ef80..58927944b 100644 --- a/test/sparse_solver.h +++ b/test/sparse_solver.h @@ -9,6 +9,7 @@ #include "sparse.h" #include +#include #include template @@ -144,6 +145,136 @@ void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, } } +// specialization of generic check_sparse_solving for SuperLU in order to also test adjoint and transpose solves +template +void check_sparse_solving(Eigen::SparseLU >& solver, const typename Eigen::SparseMatrix& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db) +{ + typedef typename Eigen::SparseMatrix Mat; + typedef typename Mat::StorageIndex StorageIndex; + typedef typename Eigen::SparseLU > Solver; + + // reference solutions computed by dense QR solver + DenseRhs refX1 = dA.householderQr().solve(db); // solution of A x = db + DenseRhs refX2 = dA.transpose().householderQr().solve(db); // solution of A^T * x = db (use transposed matrix A^T) + DenseRhs refX3 = dA.adjoint().householderQr().solve(db); // solution of A^* * x = db (use adjoint matrix A^*) + + + { + Rhs x1(A.cols(), b.cols()); + Rhs x2(A.cols(), b.cols()); + Rhs x3(A.cols(), b.cols()); + Rhs oldb = b; + + solver.compute(A); + if (solver.info() != Success) + { + std::cerr << "ERROR | sparse solver testing, factorization failed (" << typeid(Solver).name() << ")\n"; + VERIFY(solver.info() == Success); + } + x1 = solver.solve(b); + if (solver.info() != Success) + { + std::cerr << "WARNING | sparse solver testing: solving failed (" << typeid(Solver).name() << ")\n"; + return; + } + VERIFY(oldb.isApprox(b,0.0) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(x1.isApprox(refX1,test_precision())); + + // test solve with transposed + x2 = solver.transpose().solve(b); + VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(x2.isApprox(refX2,test_precision())); + + + // test solve with adjoint + //solver.template _solve_impl_transposed(b, x3); + x3 = solver.adjoint().solve(b); + VERIFY(oldb.isApprox(b,0.0) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(x3.isApprox(refX3,test_precision())); + + x1.setZero(); + solve_with_guess(solver, b, x1, x1); + VERIFY(solver.info() == Success && "solving failed when using analyzePattern/factorize API"); + VERIFY(oldb.isApprox(b,0.0) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(x1.isApprox(refX1,test_precision())); + + x1.setZero(); + x2.setZero(); + x3.setZero(); + // test the analyze/factorize API + solver.analyzePattern(A); + solver.factorize(A); + VERIFY(solver.info() == Success && "factorization failed when using analyzePattern/factorize API"); + x1 = solver.solve(b); + x2 = solver.transpose().solve(b); + x3 = solver.adjoint().solve(b); + + VERIFY(solver.info() == Success && "solving failed when using analyzePattern/factorize API"); + VERIFY(oldb.isApprox(b,0.0) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(x1.isApprox(refX1,test_precision())); + VERIFY(x2.isApprox(refX2,test_precision())); + VERIFY(x3.isApprox(refX3,test_precision())); + + x1.setZero(); + // test with Map + MappedSparseMatrix Am(A.rows(), A.cols(), A.nonZeros(), const_cast(A.outerIndexPtr()), const_cast(A.innerIndexPtr()), const_cast(A.valuePtr())); + solver.compute(Am); + VERIFY(solver.info() == Success && "factorization failed when using Map"); + DenseRhs dx(refX1); + dx.setZero(); + Map xm(dx.data(), dx.rows(), dx.cols()); + Map bm(db.data(), db.rows(), db.cols()); + xm = solver.solve(bm); + VERIFY(solver.info() == Success && "solving failed when using Map"); + VERIFY(oldb.isApprox(bm,0.0) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(xm.isApprox(refX1,test_precision())); + } + + // if not too large, do some extra check: + if(A.rows()<2000) + { + // test initialization ctor + { + Rhs x(b.rows(), b.cols()); + Solver solver2(A); + VERIFY(solver2.info() == Success); + x = solver2.solve(b); + VERIFY(x.isApprox(refX1,test_precision())); + } + + // test dense Block as the result and rhs: + { + DenseRhs x(refX1.rows(), refX1.cols()); + DenseRhs oldb(db); + x.setZero(); + x.block(0,0,x.rows(),x.cols()) = solver.solve(db.block(0,0,db.rows(),db.cols())); + VERIFY(oldb.isApprox(db,0.0) && "sparse solver testing: the rhs should not be modified!"); + VERIFY(x.isApprox(refX1,test_precision())); + } + + // test uncompressed inputs + { + Mat A2 = A; + A2.reserve((ArrayXf::Random(A.outerSize())+2).template cast().eval()); + solver.compute(A2); + Rhs x = solver.solve(b); + VERIFY(x.isApprox(refX1,test_precision())); + } + + // test expression as input + { + solver.compute(0.5*(A+A)); + Rhs x = solver.solve(b); + VERIFY(x.isApprox(refX1,test_precision())); + + Solver solver2(0.5*(A+A)); + Rhs x2 = solver2.solve(b); + VERIFY(x2.isApprox(refX1,test_precision())); + } + } +} + + template void check_sparse_solving_real_cases(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const typename Solver::MatrixType& fullA, const Rhs& refX) { -- cgit v1.2.3