11 #include <Eigen/SparseCore>
12 #include <Eigen/SparseLU>
15 template<
typename Solver,
typename Rhs,
typename Guess,
typename Result>
17 if(internal::random<bool>())
20 x =
solver.derived().solveWithGuess(
b,
g) + Result::Zero(
x.rows(),
x.cols());
25 x =
solver.derived().solveWithGuess(
b.derived(),
g);
29 template<
typename Solver,
typename Rhs,
typename Guess,
typename Result>
31 if(internal::random<bool>())
32 x =
solver.derived().solve(
b) + Result::Zero(
x.rows(),
x.cols());
37 template<
typename Solver,
typename Rhs,
typename Guess,
typename Result>
42 template<
typename Solver,
typename Rhs,
typename DenseMat,
typename DenseRhs>
47 typedef typename Mat::StorageIndex StorageIndex;
49 DenseRhs refX = dA.householderQr().solve(db);
57 std::cerr <<
"ERROR | sparse solver testing, factorization failed (" <<
typeid(Solver).
name() <<
")\n";
63 std::cerr <<
"WARNING: sparse solver testing: solving failed (" <<
typeid(Solver).
name() <<
")\n";
70 VERIFY(oldb.isApprox(
b) &&
"sparse solver testing: the rhs should not be modified!");
71 VERIFY(
x.isApprox(refX,test_precision<Scalar>()));
76 VERIFY(oldb.isApprox(
b) &&
"sparse solver testing: the rhs should not be modified!");
77 VERIFY(
x.isApprox(refX,test_precision<Scalar>()));
83 VERIFY(
solver.info() ==
Success &&
"factorization failed when using analyzePattern/factorize API");
85 VERIFY(
solver.info() ==
Success &&
"solving failed when using analyzePattern/factorize API");
86 VERIFY(oldb.isApprox(
b) &&
"sparse solver testing: the rhs should not be modified!");
87 VERIFY(
x.isApprox(refX,test_precision<Scalar>()));
100 VERIFY(oldb.isApprox(bm) &&
"sparse solver testing: the rhs should not be modified!");
101 VERIFY(xm.isApprox(refX,test_precision<Scalar>()));
109 Rhs x(
b.rows(),
b.cols());
112 x = solver2.solve(
b);
113 VERIFY(
x.isApprox(refX,test_precision<Scalar>()));
118 DenseRhs
x(refX.rows(), refX.cols());
121 x.block(0,0,
x.rows(),
x.cols()) =
solver.solve(db.block(0,0,db.rows(),db.cols()));
122 VERIFY(oldb.isApprox(db) &&
"sparse solver testing: the rhs should not be modified!");
123 VERIFY(
x.isApprox(refX,test_precision<Scalar>()));
129 A2.reserve((ArrayXf::Random(
A.outerSize())+2).template cast<typename Mat::StorageIndex>().eval());
132 VERIFY(
x.isApprox(refX,test_precision<Scalar>()));
139 VERIFY(
x.isApprox(refX,test_precision<Scalar>()));
141 Solver solver2(0.5*(
A+
A));
142 Rhs x2 = solver2.solve(
b);
143 VERIFY(
x2.isApprox(refX,test_precision<Scalar>()));
149 template<
typename Scalar,
typename Rhs,
typename DenseMat,
typename DenseRhs>
153 typedef typename Mat::StorageIndex StorageIndex;
157 DenseRhs refX1 = dA.householderQr().solve(db);
158 DenseRhs refX2 = dA.
transpose().householderQr().solve(db);
159 DenseRhs refX3 = dA.adjoint().householderQr().solve(db);
171 std::cerr <<
"ERROR | sparse solver testing, factorization failed (" <<
typeid(Solver).
name() <<
")\n";
177 std::cerr <<
"WARNING | sparse solver testing: solving failed (" <<
typeid(Solver).
name() <<
")\n";
180 VERIFY(oldb.isApprox(
b,0.0) &&
"sparse solver testing: the rhs should not be modified!");
181 VERIFY(
x1.isApprox(refX1,test_precision<Scalar>()));
185 VERIFY(oldb.isApprox(
b) &&
"sparse solver testing: the rhs should not be modified!");
186 VERIFY(
x2.isApprox(refX2,test_precision<Scalar>()));
192 VERIFY(oldb.isApprox(
b,0.0) &&
"sparse solver testing: the rhs should not be modified!");
193 VERIFY(
x3.isApprox(refX3,test_precision<Scalar>()));
197 VERIFY(
solver.info() ==
Success &&
"solving failed when using analyzePattern/factorize API");
198 VERIFY(oldb.isApprox(
b,0.0) &&
"sparse solver testing: the rhs should not be modified!");
199 VERIFY(
x1.isApprox(refX1,test_precision<Scalar>()));
207 VERIFY(
solver.info() ==
Success &&
"factorization failed when using analyzePattern/factorize API");
212 VERIFY(
solver.info() ==
Success &&
"solving failed when using analyzePattern/factorize API");
213 VERIFY(oldb.isApprox(
b,0.0) &&
"sparse solver testing: the rhs should not be modified!");
214 VERIFY(
x1.isApprox(refX1,test_precision<Scalar>()));
215 VERIFY(
x2.isApprox(refX2,test_precision<Scalar>()));
216 VERIFY(
x3.isApprox(refX3,test_precision<Scalar>()));
229 VERIFY(oldb.isApprox(bm,0.0) &&
"sparse solver testing: the rhs should not be modified!");
230 VERIFY(xm.isApprox(refX1,test_precision<Scalar>()));
238 Rhs x(
b.rows(),
b.cols());
241 x = solver2.solve(
b);
242 VERIFY(
x.isApprox(refX1,test_precision<Scalar>()));
247 DenseRhs
x(refX1.rows(), refX1.cols());
250 x.block(0,0,
x.rows(),
x.cols()) =
solver.solve(db.block(0,0,db.rows(),db.cols()));
251 VERIFY(oldb.isApprox(db,0.0) &&
"sparse solver testing: the rhs should not be modified!");
252 VERIFY(
x.isApprox(refX1,test_precision<Scalar>()));
258 A2.reserve((ArrayXf::Random(
A.outerSize())+2).template cast<typename Mat::StorageIndex>().eval());
261 VERIFY(
x.isApprox(refX1,test_precision<Scalar>()));
268 VERIFY(
x.isApprox(refX1,test_precision<Scalar>()));
270 Solver solver2(0.5*(
A+
A));
271 Rhs x2 = solver2.solve(
b);
272 VERIFY(
x2.isApprox(refX1,test_precision<Scalar>()));
278 template<
typename Solver,
typename Rhs>
285 Rhs x(
A.cols(),
b.cols());
290 std::cerr <<
"ERROR | sparse solver testing, factorization failed (" <<
typeid(Solver).
name() <<
")\n";
297 std::cerr <<
"WARNING | sparse solver testing, solving failed (" <<
typeid(Solver).
name() <<
")\n";
302 VERIFY( (res_error <= test_precision<Scalar>() ) &&
"sparse solver failed without noticing it");
305 if(refX.size() != 0 && (refX -
x).norm()/refX.norm() > test_precision<Scalar>())
307 std::cerr <<
"WARNING | found solution is different from the provided reference one\n";
311 template<
typename Solver,
typename DenseMat>
320 std::cerr <<
"WARNING | sparse solver testing: factorization failed (check_sparse_determinant)\n";
324 Scalar refDet = dA.determinant();
327 template<
typename Solver,
typename DenseMat>
337 std::cerr <<
"WARNING | sparse solver testing: factorization failed (check_sparse_abs_determinant)\n";
345 template<
typename Solver,
typename DenseMat>
352 int size = internal::random<int>(1,maxSize);
361 dA = dM * dM.adjoint();
367 halfA.template selfadjointView<Solver::UpLo>().rankUpdate(
M);
373 #ifdef TEST_REAL_CASES
374 template<
typename Scalar>
375 inline std::string get_matrixfolder()
377 std::string mat_folder = TEST_REAL_CASES;
378 if( internal::is_same<
Scalar, std::complex<float> >::
value || internal::is_same<
Scalar, std::complex<double> >::
value )
379 mat_folder = mat_folder +
static_cast<std::string
>(
"/complex/");
381 mat_folder = mat_folder +
static_cast<std::string
>(
"/real/");
384 std::string sym_to_string(
int sym)
387 if(sym==
SPD)
return "SPD ";
390 template<
typename Derived>
393 std::stringstream
ss;
397 template<
typename Derived>
408 typedef typename Mat::StorageIndex StorageIndex;
421 int rhsCols = internal::random<int>(1,16);
442 b = DenseVector::Zero(
size);
448 #ifdef TEST_REAL_CASES
452 std::string mat_folder = get_matrixfolder<Scalar>();
458 if(
A.diagonal().size() <= maxRealWorldSize)
467 halfA.template selfadjointView<Solver::UpLo>() =
A.template triangularView<Eigen::Lower>().twistedBy(pnull);
469 std::cout <<
"INFO | Testing " << sym_to_string(it.
sym()) <<
"sparse problem " << it.
matname()
470 <<
" (" <<
A.rows() <<
"x" <<
A.cols() <<
") using " <<
typeid(Solver).
name() <<
"..." << std::endl;
474 std::cout <<
"INFO | " <<
stats << std::endl;
479 std::cout <<
"INFO | Skip sparse problem \"" << it.
matname() <<
"\" (too large)" << std::endl;
506 template<
typename Solver,
typename DenseMat>
512 Index size = internal::random<int>(1,maxSize);
527 template<
class Scalar>
543 int rhsCols = internal::random<int>(1,16);
570 if(checkDeficient &&
size>1) {
579 #ifdef TEST_REAL_CASES
583 std::string mat_folder = get_matrixfolder<Scalar>();
588 if(
A.diagonal().size() <= maxRealWorldSize)
592 std::cout <<
"INFO | Testing " << sym_to_string(it.
sym()) <<
"sparse problem " << it.
matname()
593 <<
" (" <<
A.rows() <<
"x" <<
A.cols() <<
") using " <<
typeid(Solver).
name() <<
"..." << std::endl;
597 std::cout <<
"INFO | " <<
stats << std::endl;
601 std::cout <<
"INFO | SKIP sparse problem \"" << it.
matname() <<
"\" (too large)" << std::endl;
622 int size = internal::random<int>(1,30);
625 dA = (dA.array().abs()<0.3).select(0,dA);
626 dA.diagonal() = (dA.diagonal().array()==0).select(1,dA.diagonal());
650 template<
typename Solver,
typename DenseMat>
656 int rows = internal::random<int>(1,maxSize);
657 int cols = internal::random<int>(1,
rows);
674 int rhsCols = internal::random<int>(1,16);
695 b = DenseVector::Zero(
A.rows());