1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011 Gael Guennebaud <g.gael@free.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10 #include "sparse.h"
11 #include <Eigen/SparseCore>
12 #include <sstream>
13
14 template<typename Solver, typename Rhs, typename Guess,typename Result>
solve_with_guess(IterativeSolverBase<Solver> & solver,const MatrixBase<Rhs> & b,const Guess & g,Result & x)15 void solve_with_guess(IterativeSolverBase<Solver>& solver, const MatrixBase<Rhs>& b, const Guess& g, Result &x) {
16 if(internal::random<bool>())
17 {
18 // With a temporary through evaluator<SolveWithGuess>
19 x = solver.derived().solveWithGuess(b,g) + Result::Zero(x.rows(), x.cols());
20 }
21 else
22 {
23 // direct evaluation within x through Assignment<Result,SolveWithGuess>
24 x = solver.derived().solveWithGuess(b.derived(),g);
25 }
26 }
27
28 template<typename Solver, typename Rhs, typename Guess,typename Result>
solve_with_guess(SparseSolverBase<Solver> & solver,const MatrixBase<Rhs> & b,const Guess &,Result & x)29 void solve_with_guess(SparseSolverBase<Solver>& solver, const MatrixBase<Rhs>& b, const Guess& , Result& x) {
30 if(internal::random<bool>())
31 x = solver.derived().solve(b) + Result::Zero(x.rows(), x.cols());
32 else
33 x = solver.derived().solve(b);
34 }
35
36 template<typename Solver, typename Rhs, typename Guess,typename Result>
solve_with_guess(SparseSolverBase<Solver> & solver,const SparseMatrixBase<Rhs> & b,const Guess &,Result & x)37 void solve_with_guess(SparseSolverBase<Solver>& solver, const SparseMatrixBase<Rhs>& b, const Guess& , Result& x) {
38 x = solver.derived().solve(b);
39 }
40
41 template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs>
check_sparse_solving(Solver & solver,const typename Solver::MatrixType & A,const Rhs & b,const DenseMat & dA,const DenseRhs & db)42 void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db)
43 {
44 typedef typename Solver::MatrixType Mat;
45 typedef typename Mat::Scalar Scalar;
46 typedef typename Mat::StorageIndex StorageIndex;
47
48 DenseRhs refX = dA.householderQr().solve(db);
49 {
50 Rhs x(A.cols(), b.cols());
51 Rhs oldb = b;
52
53 solver.compute(A);
54 if (solver.info() != Success)
55 {
56 std::cerr << "ERROR | sparse solver testing, factorization failed (" << typeid(Solver).name() << ")\n";
57 VERIFY(solver.info() == Success);
58 }
59 x = solver.solve(b);
60 if (solver.info() != Success)
61 {
62 std::cerr << "WARNING | sparse solver testing: solving failed (" << typeid(Solver).name() << ")\n";
63 return;
64 }
65 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!");
66 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
67
68 x.setZero();
69 solve_with_guess(solver, b, x, x);
70 VERIFY(solver.info() == Success && "solving failed when using analyzePattern/factorize API");
71 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!");
72 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
73
74 x.setZero();
75 // test the analyze/factorize API
76 solver.analyzePattern(A);
77 solver.factorize(A);
78 VERIFY(solver.info() == Success && "factorization failed when using analyzePattern/factorize API");
79 x = solver.solve(b);
80 VERIFY(solver.info() == Success && "solving failed when using analyzePattern/factorize API");
81 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!");
82 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
83
84 x.setZero();
85 // test with Map
86 MappedSparseMatrix<Scalar,Mat::Options,StorageIndex> Am(A.rows(), A.cols(), A.nonZeros(), const_cast<StorageIndex*>(A.outerIndexPtr()), const_cast<StorageIndex*>(A.innerIndexPtr()), const_cast<Scalar*>(A.valuePtr()));
87 solver.compute(Am);
88 VERIFY(solver.info() == Success && "factorization failed when using Map");
89 DenseRhs dx(refX);
90 dx.setZero();
91 Map<DenseRhs> xm(dx.data(), dx.rows(), dx.cols());
92 Map<const DenseRhs> bm(db.data(), db.rows(), db.cols());
93 xm = solver.solve(bm);
94 VERIFY(solver.info() == Success && "solving failed when using Map");
95 VERIFY(oldb.isApprox(bm) && "sparse solver testing: the rhs should not be modified!");
96 VERIFY(xm.isApprox(refX,test_precision<Scalar>()));
97 }
98
99 // if not too large, do some extra check:
100 if(A.rows()<2000)
101 {
102 // test initialization ctor
103 {
104 Rhs x(b.rows(), b.cols());
105 Solver solver2(A);
106 VERIFY(solver2.info() == Success);
107 x = solver2.solve(b);
108 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
109 }
110
111 // test dense Block as the result and rhs:
112 {
113 DenseRhs x(refX.rows(), refX.cols());
114 DenseRhs oldb(db);
115 x.setZero();
116 x.block(0,0,x.rows(),x.cols()) = solver.solve(db.block(0,0,db.rows(),db.cols()));
117 VERIFY(oldb.isApprox(db) && "sparse solver testing: the rhs should not be modified!");
118 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
119 }
120
121 // test uncompressed inputs
122 {
123 Mat A2 = A;
124 A2.reserve((ArrayXf::Random(A.outerSize())+2).template cast<typename Mat::StorageIndex>().eval());
125 solver.compute(A2);
126 Rhs x = solver.solve(b);
127 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
128 }
129
130 // test expression as input
131 {
132 solver.compute(0.5*(A+A));
133 Rhs x = solver.solve(b);
134 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
135
136 Solver solver2(0.5*(A+A));
137 Rhs x2 = solver2.solve(b);
138 VERIFY(x2.isApprox(refX,test_precision<Scalar>()));
139 }
140 }
141 }
142
143 template<typename Solver, typename Rhs>
check_sparse_solving_real_cases(Solver & solver,const typename Solver::MatrixType & A,const Rhs & b,const typename Solver::MatrixType & fullA,const Rhs & refX)144 void check_sparse_solving_real_cases(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const typename Solver::MatrixType& fullA, const Rhs& refX)
145 {
146 typedef typename Solver::MatrixType Mat;
147 typedef typename Mat::Scalar Scalar;
148 typedef typename Mat::RealScalar RealScalar;
149
150 Rhs x(A.cols(), b.cols());
151
152 solver.compute(A);
153 if (solver.info() != Success)
154 {
155 std::cerr << "ERROR | sparse solver testing, factorization failed (" << typeid(Solver).name() << ")\n";
156 VERIFY(solver.info() == Success);
157 }
158 x = solver.solve(b);
159
160 if (solver.info() != Success)
161 {
162 std::cerr << "WARNING | sparse solver testing, solving failed (" << typeid(Solver).name() << ")\n";
163 return;
164 }
165
166 RealScalar res_error = (fullA*x-b).norm()/b.norm();
167 VERIFY( (res_error <= test_precision<Scalar>() ) && "sparse solver failed without noticing it");
168
169
170 if(refX.size() != 0 && (refX - x).norm()/refX.norm() > test_precision<Scalar>())
171 {
172 std::cerr << "WARNING | found solution is different from the provided reference one\n";
173 }
174
175 }
176 template<typename Solver, typename DenseMat>
check_sparse_determinant(Solver & solver,const typename Solver::MatrixType & A,const DenseMat & dA)177 void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
178 {
179 typedef typename Solver::MatrixType Mat;
180 typedef typename Mat::Scalar Scalar;
181
182 solver.compute(A);
183 if (solver.info() != Success)
184 {
185 std::cerr << "WARNING | sparse solver testing: factorization failed (check_sparse_determinant)\n";
186 return;
187 }
188
189 Scalar refDet = dA.determinant();
190 VERIFY_IS_APPROX(refDet,solver.determinant());
191 }
192 template<typename Solver, typename DenseMat>
check_sparse_abs_determinant(Solver & solver,const typename Solver::MatrixType & A,const DenseMat & dA)193 void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
194 {
195 using std::abs;
196 typedef typename Solver::MatrixType Mat;
197 typedef typename Mat::Scalar Scalar;
198
199 solver.compute(A);
200 if (solver.info() != Success)
201 {
202 std::cerr << "WARNING | sparse solver testing: factorization failed (check_sparse_abs_determinant)\n";
203 return;
204 }
205
206 Scalar refDet = abs(dA.determinant());
207 VERIFY_IS_APPROX(refDet,solver.absDeterminant());
208 }
209
210 template<typename Solver, typename DenseMat>
211 int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300)
212 {
213 typedef typename Solver::MatrixType Mat;
214 typedef typename Mat::Scalar Scalar;
215 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
216
217 int size = internal::random<int>(1,maxSize);
218 double density = (std::max)(8./(size*size), 0.01);
219
220 Mat M(size, size);
221 DenseMatrix dM(size, size);
222
223 initSparse<Scalar>(density, dM, M, ForceNonZeroDiag);
224
225 A = M * M.adjoint();
226 dA = dM * dM.adjoint();
227
228 halfA.resize(size,size);
229 if(Solver::UpLo==(Lower|Upper))
230 halfA = A;
231 else
232 halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M);
233
234 return size;
235 }
236
237
238 #ifdef TEST_REAL_CASES
239 template<typename Scalar>
get_matrixfolder()240 inline std::string get_matrixfolder()
241 {
242 std::string mat_folder = TEST_REAL_CASES;
243 if( internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value )
244 mat_folder = mat_folder + static_cast<std::string>("/complex/");
245 else
246 mat_folder = mat_folder + static_cast<std::string>("/real/");
247 return mat_folder;
248 }
sym_to_string(int sym)249 std::string sym_to_string(int sym)
250 {
251 if(sym==Symmetric) return "Symmetric ";
252 if(sym==SPD) return "SPD ";
253 return "";
254 }
255 template<typename Derived>
solver_stats(const IterativeSolverBase<Derived> & solver)256 std::string solver_stats(const IterativeSolverBase<Derived> &solver)
257 {
258 std::stringstream ss;
259 ss << solver.iterations() << " iters, error: " << solver.error();
260 return ss.str();
261 }
262 template<typename Derived>
solver_stats(const SparseSolverBase<Derived> &)263 std::string solver_stats(const SparseSolverBase<Derived> &/*solver*/)
264 {
265 return "";
266 }
267 #endif
268
269 template<typename Solver> void check_sparse_spd_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000)
270 {
271 typedef typename Solver::MatrixType Mat;
272 typedef typename Mat::Scalar Scalar;
273 typedef typename Mat::StorageIndex StorageIndex;
274 typedef SparseMatrix<Scalar,ColMajor, StorageIndex> SpMat;
275 typedef SparseVector<Scalar, 0, StorageIndex> SpVec;
276 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
277 typedef Matrix<Scalar,Dynamic,1> DenseVector;
278
279 // generate the problem
280 Mat A, halfA;
281 DenseMatrix dA;
282 for (int i = 0; i < g_repeat; i++) {
283 int size = generate_sparse_spd_problem(solver, A, halfA, dA, maxSize);
284
285 // generate the right hand sides
286 int rhsCols = internal::random<int>(1,16);
287 double density = (std::max)(8./(size*rhsCols), 0.1);
288 SpMat B(size,rhsCols);
289 DenseVector b = DenseVector::Random(size);
290 DenseMatrix dB(size,rhsCols);
291 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
292 SpVec c = B.col(0);
293 DenseVector dc = dB.col(0);
294
295 CALL_SUBTEST( check_sparse_solving(solver, A, b, dA, b) );
296 CALL_SUBTEST( check_sparse_solving(solver, halfA, b, dA, b) );
297 CALL_SUBTEST( check_sparse_solving(solver, A, dB, dA, dB) );
298 CALL_SUBTEST( check_sparse_solving(solver, halfA, dB, dA, dB) );
299 CALL_SUBTEST( check_sparse_solving(solver, A, B, dA, dB) );
300 CALL_SUBTEST( check_sparse_solving(solver, halfA, B, dA, dB) );
301 CALL_SUBTEST( check_sparse_solving(solver, A, c, dA, dc) );
302 CALL_SUBTEST( check_sparse_solving(solver, halfA, c, dA, dc) );
303
304 // check only once
305 if(i==0)
306 {
307 b = DenseVector::Zero(size);
308 check_sparse_solving(solver, A, b, dA, b);
309 }
310 }
311
312 // First, get the folder
313 #ifdef TEST_REAL_CASES
314 // Test real problems with double precision only
315 if (internal::is_same<typename NumTraits<Scalar>::Real, double>::value)
316 {
317 std::string mat_folder = get_matrixfolder<Scalar>();
318 MatrixMarketIterator<Scalar> it(mat_folder);
319 for (; it; ++it)
320 {
321 if (it.sym() == SPD){
322 A = it.matrix();
323 if(A.diagonal().size() <= maxRealWorldSize)
324 {
325 DenseVector b = it.rhs();
326 DenseVector refX = it.refX();
327 PermutationMatrix<Dynamic, Dynamic, StorageIndex> pnull;
328 halfA.resize(A.rows(), A.cols());
329 if(Solver::UpLo == (Lower|Upper))
330 halfA = A;
331 else
332 halfA.template selfadjointView<Solver::UpLo>() = A.template triangularView<Eigen::Lower>().twistedBy(pnull);
333
334 std::cout << "INFO | Testing " << sym_to_string(it.sym()) << "sparse problem " << it.matname()
335 << " (" << A.rows() << "x" << A.cols() << ") using " << typeid(Solver).name() << "..." << std::endl;
336 CALL_SUBTEST( check_sparse_solving_real_cases(solver, A, b, A, refX) );
337 std::string stats = solver_stats(solver);
338 if(stats.size()>0)
339 std::cout << "INFO | " << stats << std::endl;
340 CALL_SUBTEST( check_sparse_solving_real_cases(solver, halfA, b, A, refX) );
341 }
342 else
343 {
344 std::cout << "INFO | Skip sparse problem \"" << it.matname() << "\" (too large)" << std::endl;
345 }
346 }
347 }
348 }
349 #else
350 EIGEN_UNUSED_VARIABLE(maxRealWorldSize);
351 #endif
352 }
353
check_sparse_spd_determinant(Solver & solver)354 template<typename Solver> void check_sparse_spd_determinant(Solver& solver)
355 {
356 typedef typename Solver::MatrixType Mat;
357 typedef typename Mat::Scalar Scalar;
358 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
359
360 // generate the problem
361 Mat A, halfA;
362 DenseMatrix dA;
363 generate_sparse_spd_problem(solver, A, halfA, dA, 30);
364
365 for (int i = 0; i < g_repeat; i++) {
366 check_sparse_determinant(solver, A, dA);
367 check_sparse_determinant(solver, halfA, dA );
368 }
369 }
370
371 template<typename Solver, typename DenseMat>
372 Index generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300, int options = ForceNonZeroDiag)
373 {
374 typedef typename Solver::MatrixType Mat;
375 typedef typename Mat::Scalar Scalar;
376
377 Index size = internal::random<int>(1,maxSize);
378 double density = (std::max)(8./(size*size), 0.01);
379
380 A.resize(size,size);
381 dA.resize(size,size);
382
383 initSparse<Scalar>(density, dA, A, options);
384
385 return size;
386 }
387
388
389 struct prune_column {
390 Index m_col;
prune_columnprune_column391 prune_column(Index col) : m_col(col) {}
392 template<class Scalar>
operatorprune_column393 bool operator()(Index, Index col, const Scalar&) const {
394 return col != m_col;
395 }
396 };
397
398
399 template<typename Solver> void check_sparse_square_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000, bool checkDeficient = false)
400 {
401 typedef typename Solver::MatrixType Mat;
402 typedef typename Mat::Scalar Scalar;
403 typedef SparseMatrix<Scalar,ColMajor, typename Mat::StorageIndex> SpMat;
404 typedef SparseVector<Scalar, 0, typename Mat::StorageIndex> SpVec;
405 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
406 typedef Matrix<Scalar,Dynamic,1> DenseVector;
407
408 int rhsCols = internal::random<int>(1,16);
409
410 Mat A;
411 DenseMatrix dA;
412 for (int i = 0; i < g_repeat; i++) {
413 Index size = generate_sparse_square_problem(solver, A, dA, maxSize);
414
415 A.makeCompressed();
416 DenseVector b = DenseVector::Random(size);
417 DenseMatrix dB(size,rhsCols);
418 SpMat B(size,rhsCols);
419 double density = (std::max)(8./(size*rhsCols), 0.1);
420 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
421 B.makeCompressed();
422 SpVec c = B.col(0);
423 DenseVector dc = dB.col(0);
424 CALL_SUBTEST(check_sparse_solving(solver, A, b, dA, b));
425 CALL_SUBTEST(check_sparse_solving(solver, A, dB, dA, dB));
426 CALL_SUBTEST(check_sparse_solving(solver, A, B, dA, dB));
427 CALL_SUBTEST(check_sparse_solving(solver, A, c, dA, dc));
428
429 // check only once
430 if(i==0)
431 {
432 b = DenseVector::Zero(size);
433 check_sparse_solving(solver, A, b, dA, b);
434 }
435 // regression test for Bug 792 (structurally rank deficient matrices):
436 if(checkDeficient && size>1) {
437 Index col = internal::random<int>(0,int(size-1));
438 A.prune(prune_column(col));
439 solver.compute(A);
440 VERIFY_IS_EQUAL(solver.info(), NumericalIssue);
441 }
442 }
443
444 // First, get the folder
445 #ifdef TEST_REAL_CASES
446 // Test real problems with double precision only
447 if (internal::is_same<typename NumTraits<Scalar>::Real, double>::value)
448 {
449 std::string mat_folder = get_matrixfolder<Scalar>();
450 MatrixMarketIterator<Scalar> it(mat_folder);
451 for (; it; ++it)
452 {
453 A = it.matrix();
454 if(A.diagonal().size() <= maxRealWorldSize)
455 {
456 DenseVector b = it.rhs();
457 DenseVector refX = it.refX();
458 std::cout << "INFO | Testing " << sym_to_string(it.sym()) << "sparse problem " << it.matname()
459 << " (" << A.rows() << "x" << A.cols() << ") using " << typeid(Solver).name() << "..." << std::endl;
460 CALL_SUBTEST(check_sparse_solving_real_cases(solver, A, b, A, refX));
461 std::string stats = solver_stats(solver);
462 if(stats.size()>0)
463 std::cout << "INFO | " << stats << std::endl;
464 }
465 else
466 {
467 std::cout << "INFO | SKIP sparse problem \"" << it.matname() << "\" (too large)" << std::endl;
468 }
469 }
470 }
471 #else
472 EIGEN_UNUSED_VARIABLE(maxRealWorldSize);
473 #endif
474
475 }
476
check_sparse_square_determinant(Solver & solver)477 template<typename Solver> void check_sparse_square_determinant(Solver& solver)
478 {
479 typedef typename Solver::MatrixType Mat;
480 typedef typename Mat::Scalar Scalar;
481 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
482
483 for (int i = 0; i < g_repeat; i++) {
484 // generate the problem
485 Mat A;
486 DenseMatrix dA;
487
488 int size = internal::random<int>(1,30);
489 dA.setRandom(size,size);
490
491 dA = (dA.array().abs()<0.3).select(0,dA);
492 dA.diagonal() = (dA.diagonal().array()==0).select(1,dA.diagonal());
493 A = dA.sparseView();
494 A.makeCompressed();
495
496 check_sparse_determinant(solver, A, dA);
497 }
498 }
499
check_sparse_square_abs_determinant(Solver & solver)500 template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver)
501 {
502 typedef typename Solver::MatrixType Mat;
503 typedef typename Mat::Scalar Scalar;
504 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
505
506 for (int i = 0; i < g_repeat; i++) {
507 // generate the problem
508 Mat A;
509 DenseMatrix dA;
510 generate_sparse_square_problem(solver, A, dA, 30);
511 A.makeCompressed();
512 check_sparse_abs_determinant(solver, A, dA);
513 }
514 }
515
516 template<typename Solver, typename DenseMat>
517 void generate_sparse_leastsquare_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300, int options = ForceNonZeroDiag)
518 {
519 typedef typename Solver::MatrixType Mat;
520 typedef typename Mat::Scalar Scalar;
521
522 int rows = internal::random<int>(1,maxSize);
523 int cols = internal::random<int>(1,rows);
524 double density = (std::max)(8./(rows*cols), 0.01);
525
526 A.resize(rows,cols);
527 dA.resize(rows,cols);
528
529 initSparse<Scalar>(density, dA, A, options);
530 }
531
check_sparse_leastsquare_solving(Solver & solver)532 template<typename Solver> void check_sparse_leastsquare_solving(Solver& solver)
533 {
534 typedef typename Solver::MatrixType Mat;
535 typedef typename Mat::Scalar Scalar;
536 typedef SparseMatrix<Scalar,ColMajor, typename Mat::StorageIndex> SpMat;
537 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
538 typedef Matrix<Scalar,Dynamic,1> DenseVector;
539
540 int rhsCols = internal::random<int>(1,16);
541
542 Mat A;
543 DenseMatrix dA;
544 for (int i = 0; i < g_repeat; i++) {
545 generate_sparse_leastsquare_problem(solver, A, dA);
546
547 A.makeCompressed();
548 DenseVector b = DenseVector::Random(A.rows());
549 DenseMatrix dB(A.rows(),rhsCols);
550 SpMat B(A.rows(),rhsCols);
551 double density = (std::max)(8./(A.rows()*rhsCols), 0.1);
552 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
553 B.makeCompressed();
554 check_sparse_solving(solver, A, b, dA, b);
555 check_sparse_solving(solver, A, dB, dA, dB);
556 check_sparse_solving(solver, A, B, dA, dB);
557
558 // check only once
559 if(i==0)
560 {
561 b = DenseVector::Zero(A.rows());
562 check_sparse_solving(solver, A, b, dA, b);
563 }
564 }
565 }
566