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 <Eigen/SparseLU>
13 #include <sstream>
14
15 template<typename Solver, typename Rhs, typename Guess,typename Result>
solve_with_guess(IterativeSolverBase<Solver> & solver,const MatrixBase<Rhs> & b,const Guess & g,Result & x)16 void solve_with_guess(IterativeSolverBase<Solver>& solver, const MatrixBase<Rhs>& b, const Guess& g, Result &x) {
17 if(internal::random<bool>())
18 {
19 // With a temporary through evaluator<SolveWithGuess>
20 x = solver.derived().solveWithGuess(b,g) + Result::Zero(x.rows(), x.cols());
21 }
22 else
23 {
24 // direct evaluation within x through Assignment<Result,SolveWithGuess>
25 x = solver.derived().solveWithGuess(b.derived(),g);
26 }
27 }
28
29 template<typename Solver, typename Rhs, typename Guess,typename Result>
solve_with_guess(SparseSolverBase<Solver> & solver,const MatrixBase<Rhs> & b,const Guess &,Result & x)30 void solve_with_guess(SparseSolverBase<Solver>& solver, const MatrixBase<Rhs>& b, const Guess& , Result& x) {
31 if(internal::random<bool>())
32 x = solver.derived().solve(b) + Result::Zero(x.rows(), x.cols());
33 else
34 x = solver.derived().solve(b);
35 }
36
37 template<typename Solver, typename Rhs, typename Guess,typename Result>
solve_with_guess(SparseSolverBase<Solver> & solver,const SparseMatrixBase<Rhs> & b,const Guess &,Result & x)38 void solve_with_guess(SparseSolverBase<Solver>& solver, const SparseMatrixBase<Rhs>& b, const Guess& , Result& x) {
39 x = solver.derived().solve(b);
40 }
41
42 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)43 void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db)
44 {
45 typedef typename Solver::MatrixType Mat;
46 typedef typename Mat::Scalar Scalar;
47 typedef typename Mat::StorageIndex StorageIndex;
48
49 DenseRhs refX = dA.householderQr().solve(db);
50 {
51 Rhs x(A.cols(), b.cols());
52 Rhs oldb = b;
53
54 solver.compute(A);
55 if (solver.info() != Success)
56 {
57 std::cerr << "ERROR | sparse solver testing, factorization failed (" << typeid(Solver).name() << ")\n";
58 VERIFY(solver.info() == Success);
59 }
60 x = solver.solve(b);
61 if (solver.info() != Success)
62 {
63 std::cerr << "WARNING: sparse solver testing: solving failed (" << typeid(Solver).name() << ")\n";
64 // dump call stack:
65 g_test_level++;
66 VERIFY(solver.info() == Success);
67 g_test_level--;
68 return;
69 }
70 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!");
71 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
72
73 x.setZero();
74 solve_with_guess(solver, b, x, x);
75 VERIFY(solver.info() == Success && "solving failed when using solve_with_guess API");
76 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!");
77 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
78
79 x.setZero();
80 // test the analyze/factorize API
81 solver.analyzePattern(A);
82 solver.factorize(A);
83 VERIFY(solver.info() == Success && "factorization failed when using analyzePattern/factorize API");
84 x = solver.solve(b);
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>()));
88
89 x.setZero();
90 // test with Map
91 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()));
92 solver.compute(Am);
93 VERIFY(solver.info() == Success && "factorization failed when using Map");
94 DenseRhs dx(refX);
95 dx.setZero();
96 Map<DenseRhs> xm(dx.data(), dx.rows(), dx.cols());
97 Map<const DenseRhs> bm(db.data(), db.rows(), db.cols());
98 xm = solver.solve(bm);
99 VERIFY(solver.info() == Success && "solving failed when using Map");
100 VERIFY(oldb.isApprox(bm) && "sparse solver testing: the rhs should not be modified!");
101 VERIFY(xm.isApprox(refX,test_precision<Scalar>()));
102 }
103
104 // if not too large, do some extra check:
105 if(A.rows()<2000)
106 {
107 // test initialization ctor
108 {
109 Rhs x(b.rows(), b.cols());
110 Solver solver2(A);
111 VERIFY(solver2.info() == Success);
112 x = solver2.solve(b);
113 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
114 }
115
116 // test dense Block as the result and rhs:
117 {
118 DenseRhs x(refX.rows(), refX.cols());
119 DenseRhs oldb(db);
120 x.setZero();
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>()));
124 }
125
126 // test uncompressed inputs
127 {
128 Mat A2 = A;
129 A2.reserve((ArrayXf::Random(A.outerSize())+2).template cast<typename Mat::StorageIndex>().eval());
130 solver.compute(A2);
131 Rhs x = solver.solve(b);
132 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
133 }
134
135 // test expression as input
136 {
137 solver.compute(0.5*(A+A));
138 Rhs x = solver.solve(b);
139 VERIFY(x.isApprox(refX,test_precision<Scalar>()));
140
141 Solver solver2(0.5*(A+A));
142 Rhs x2 = solver2.solve(b);
143 VERIFY(x2.isApprox(refX,test_precision<Scalar>()));
144 }
145 }
146 }
147
148 // specialization of generic check_sparse_solving for SuperLU in order to also test adjoint and transpose solves
149 template<typename Scalar, typename Rhs, typename DenseMat, typename DenseRhs>
check_sparse_solving(Eigen::SparseLU<Eigen::SparseMatrix<Scalar>> & solver,const typename Eigen::SparseMatrix<Scalar> & A,const Rhs & b,const DenseMat & dA,const DenseRhs & db)150 void check_sparse_solving(Eigen::SparseLU<Eigen::SparseMatrix<Scalar> >& solver, const typename Eigen::SparseMatrix<Scalar>& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db)
151 {
152 typedef typename Eigen::SparseMatrix<Scalar> Mat;
153 typedef typename Mat::StorageIndex StorageIndex;
154 typedef typename Eigen::SparseLU<Eigen::SparseMatrix<Scalar> > Solver;
155
156 // reference solutions computed by dense QR solver
157 DenseRhs refX1 = dA.householderQr().solve(db); // solution of A x = db
158 DenseRhs refX2 = dA.transpose().householderQr().solve(db); // solution of A^T * x = db (use transposed matrix A^T)
159 DenseRhs refX3 = dA.adjoint().householderQr().solve(db); // solution of A^* * x = db (use adjoint matrix A^*)
160
161
162 {
163 Rhs x1(A.cols(), b.cols());
164 Rhs x2(A.cols(), b.cols());
165 Rhs x3(A.cols(), b.cols());
166 Rhs oldb = b;
167
168 solver.compute(A);
169 if (solver.info() != Success)
170 {
171 std::cerr << "ERROR | sparse solver testing, factorization failed (" << typeid(Solver).name() << ")\n";
172 VERIFY(solver.info() == Success);
173 }
174 x1 = solver.solve(b);
175 if (solver.info() != Success)
176 {
177 std::cerr << "WARNING | sparse solver testing: solving failed (" << typeid(Solver).name() << ")\n";
178 return;
179 }
180 VERIFY(oldb.isApprox(b,0.0) && "sparse solver testing: the rhs should not be modified!");
181 VERIFY(x1.isApprox(refX1,test_precision<Scalar>()));
182
183 // test solve with transposed
184 x2 = solver.transpose().solve(b);
185 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!");
186 VERIFY(x2.isApprox(refX2,test_precision<Scalar>()));
187
188
189 // test solve with adjoint
190 //solver.template _solve_impl_transposed<true>(b, x3);
191 x3 = solver.adjoint().solve(b);
192 VERIFY(oldb.isApprox(b,0.0) && "sparse solver testing: the rhs should not be modified!");
193 VERIFY(x3.isApprox(refX3,test_precision<Scalar>()));
194
195 x1.setZero();
196 solve_with_guess(solver, b, x1, x1);
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>()));
200
201 x1.setZero();
202 x2.setZero();
203 x3.setZero();
204 // test the analyze/factorize API
205 solver.analyzePattern(A);
206 solver.factorize(A);
207 VERIFY(solver.info() == Success && "factorization failed when using analyzePattern/factorize API");
208 x1 = solver.solve(b);
209 x2 = solver.transpose().solve(b);
210 x3 = solver.adjoint().solve(b);
211
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>()));
217
218 x1.setZero();
219 // test with Map
220 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()));
221 solver.compute(Am);
222 VERIFY(solver.info() == Success && "factorization failed when using Map");
223 DenseRhs dx(refX1);
224 dx.setZero();
225 Map<DenseRhs> xm(dx.data(), dx.rows(), dx.cols());
226 Map<const DenseRhs> bm(db.data(), db.rows(), db.cols());
227 xm = solver.solve(bm);
228 VERIFY(solver.info() == Success && "solving failed when using Map");
229 VERIFY(oldb.isApprox(bm,0.0) && "sparse solver testing: the rhs should not be modified!");
230 VERIFY(xm.isApprox(refX1,test_precision<Scalar>()));
231 }
232
233 // if not too large, do some extra check:
234 if(A.rows()<2000)
235 {
236 // test initialization ctor
237 {
238 Rhs x(b.rows(), b.cols());
239 Solver solver2(A);
240 VERIFY(solver2.info() == Success);
241 x = solver2.solve(b);
242 VERIFY(x.isApprox(refX1,test_precision<Scalar>()));
243 }
244
245 // test dense Block as the result and rhs:
246 {
247 DenseRhs x(refX1.rows(), refX1.cols());
248 DenseRhs oldb(db);
249 x.setZero();
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>()));
253 }
254
255 // test uncompressed inputs
256 {
257 Mat A2 = A;
258 A2.reserve((ArrayXf::Random(A.outerSize())+2).template cast<typename Mat::StorageIndex>().eval());
259 solver.compute(A2);
260 Rhs x = solver.solve(b);
261 VERIFY(x.isApprox(refX1,test_precision<Scalar>()));
262 }
263
264 // test expression as input
265 {
266 solver.compute(0.5*(A+A));
267 Rhs x = solver.solve(b);
268 VERIFY(x.isApprox(refX1,test_precision<Scalar>()));
269
270 Solver solver2(0.5*(A+A));
271 Rhs x2 = solver2.solve(b);
272 VERIFY(x2.isApprox(refX1,test_precision<Scalar>()));
273 }
274 }
275 }
276
277
278 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)279 void check_sparse_solving_real_cases(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const typename Solver::MatrixType& fullA, const Rhs& refX)
280 {
281 typedef typename Solver::MatrixType Mat;
282 typedef typename Mat::Scalar Scalar;
283 typedef typename Mat::RealScalar RealScalar;
284
285 Rhs x(A.cols(), b.cols());
286
287 solver.compute(A);
288 if (solver.info() != Success)
289 {
290 std::cerr << "ERROR | sparse solver testing, factorization failed (" << typeid(Solver).name() << ")\n";
291 VERIFY(solver.info() == Success);
292 }
293 x = solver.solve(b);
294
295 if (solver.info() != Success)
296 {
297 std::cerr << "WARNING | sparse solver testing, solving failed (" << typeid(Solver).name() << ")\n";
298 return;
299 }
300
301 RealScalar res_error = (fullA*x-b).norm()/b.norm();
302 VERIFY( (res_error <= test_precision<Scalar>() ) && "sparse solver failed without noticing it");
303
304
305 if(refX.size() != 0 && (refX - x).norm()/refX.norm() > test_precision<Scalar>())
306 {
307 std::cerr << "WARNING | found solution is different from the provided reference one\n";
308 }
309
310 }
311 template<typename Solver, typename DenseMat>
check_sparse_determinant(Solver & solver,const typename Solver::MatrixType & A,const DenseMat & dA)312 void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
313 {
314 typedef typename Solver::MatrixType Mat;
315 typedef typename Mat::Scalar Scalar;
316
317 solver.compute(A);
318 if (solver.info() != Success)
319 {
320 std::cerr << "WARNING | sparse solver testing: factorization failed (check_sparse_determinant)\n";
321 return;
322 }
323
324 Scalar refDet = dA.determinant();
325 VERIFY_IS_APPROX(refDet,solver.determinant());
326 }
327 template<typename Solver, typename DenseMat>
check_sparse_abs_determinant(Solver & solver,const typename Solver::MatrixType & A,const DenseMat & dA)328 void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
329 {
330 using std::abs;
331 typedef typename Solver::MatrixType Mat;
332 typedef typename Mat::Scalar Scalar;
333
334 solver.compute(A);
335 if (solver.info() != Success)
336 {
337 std::cerr << "WARNING | sparse solver testing: factorization failed (check_sparse_abs_determinant)\n";
338 return;
339 }
340
341 Scalar refDet = abs(dA.determinant());
342 VERIFY_IS_APPROX(refDet,solver.absDeterminant());
343 }
344
345 template<typename Solver, typename DenseMat>
346 int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300)
347 {
348 typedef typename Solver::MatrixType Mat;
349 typedef typename Mat::Scalar Scalar;
350 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
351
352 int size = internal::random<int>(1,maxSize);
353 double density = (std::max)(8./(size*size), 0.01);
354
355 Mat M(size, size);
356 DenseMatrix dM(size, size);
357
358 initSparse<Scalar>(density, dM, M, ForceNonZeroDiag);
359
360 A = M * M.adjoint();
361 dA = dM * dM.adjoint();
362
363 halfA.resize(size,size);
364 if(Solver::UpLo==(Lower|Upper))
365 halfA = A;
366 else
367 halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M);
368
369 return size;
370 }
371
372
373 #ifdef TEST_REAL_CASES
374 template<typename Scalar>
get_matrixfolder()375 inline std::string get_matrixfolder()
376 {
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/");
380 else
381 mat_folder = mat_folder + static_cast<std::string>("/real/");
382 return mat_folder;
383 }
sym_to_string(int sym)384 std::string sym_to_string(int sym)
385 {
386 if(sym==Symmetric) return "Symmetric ";
387 if(sym==SPD) return "SPD ";
388 return "";
389 }
390 template<typename Derived>
solver_stats(const IterativeSolverBase<Derived> & solver)391 std::string solver_stats(const IterativeSolverBase<Derived> &solver)
392 {
393 std::stringstream ss;
394 ss << solver.iterations() << " iters, error: " << solver.error();
395 return ss.str();
396 }
397 template<typename Derived>
solver_stats(const SparseSolverBase<Derived> &)398 std::string solver_stats(const SparseSolverBase<Derived> &/*solver*/)
399 {
400 return "";
401 }
402 #endif
403
404 template<typename Solver> void check_sparse_spd_solving(Solver& solver, int maxSize = (std::min)(300,EIGEN_TEST_MAX_SIZE), int maxRealWorldSize = 100000)
405 {
406 typedef typename Solver::MatrixType Mat;
407 typedef typename Mat::Scalar Scalar;
408 typedef typename Mat::StorageIndex StorageIndex;
409 typedef SparseMatrix<Scalar,ColMajor, StorageIndex> SpMat;
410 typedef SparseVector<Scalar, 0, StorageIndex> SpVec;
411 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
412 typedef Matrix<Scalar,Dynamic,1> DenseVector;
413
414 // generate the problem
415 Mat A, halfA;
416 DenseMatrix dA;
417 for (int i = 0; i < g_repeat; i++) {
418 int size = generate_sparse_spd_problem(solver, A, halfA, dA, maxSize);
419
420 // generate the right hand sides
421 int rhsCols = internal::random<int>(1,16);
422 double density = (std::max)(8./(size*rhsCols), 0.1);
423 SpMat B(size,rhsCols);
424 DenseVector b = DenseVector::Random(size);
425 DenseMatrix dB(size,rhsCols);
426 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
427 SpVec c = B.col(0);
428 DenseVector dc = dB.col(0);
429
430 CALL_SUBTEST( check_sparse_solving(solver, A, b, dA, b) );
431 CALL_SUBTEST( check_sparse_solving(solver, halfA, b, dA, b) );
432 CALL_SUBTEST( check_sparse_solving(solver, A, dB, dA, dB) );
433 CALL_SUBTEST( check_sparse_solving(solver, halfA, dB, dA, dB) );
434 CALL_SUBTEST( check_sparse_solving(solver, A, B, dA, dB) );
435 CALL_SUBTEST( check_sparse_solving(solver, halfA, B, dA, dB) );
436 CALL_SUBTEST( check_sparse_solving(solver, A, c, dA, dc) );
437 CALL_SUBTEST( check_sparse_solving(solver, halfA, c, dA, dc) );
438
439 // check only once
440 if(i==0)
441 {
442 b = DenseVector::Zero(size);
443 check_sparse_solving(solver, A, b, dA, b);
444 }
445 }
446
447 // First, get the folder
448 #ifdef TEST_REAL_CASES
449 // Test real problems with double precision only
450 if (internal::is_same<typename NumTraits<Scalar>::Real, double>::value)
451 {
452 std::string mat_folder = get_matrixfolder<Scalar>();
453 MatrixMarketIterator<Scalar> it(mat_folder);
454 for (; it; ++it)
455 {
456 if (it.sym() == SPD){
457 A = it.matrix();
458 if(A.diagonal().size() <= maxRealWorldSize)
459 {
460 DenseVector b = it.rhs();
461 DenseVector refX = it.refX();
462 PermutationMatrix<Dynamic, Dynamic, StorageIndex> pnull;
463 halfA.resize(A.rows(), A.cols());
464 if(Solver::UpLo == (Lower|Upper))
465 halfA = A;
466 else
467 halfA.template selfadjointView<Solver::UpLo>() = A.template triangularView<Eigen::Lower>().twistedBy(pnull);
468
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;
471 CALL_SUBTEST( check_sparse_solving_real_cases(solver, A, b, A, refX) );
472 std::string stats = solver_stats(solver);
473 if(stats.size()>0)
474 std::cout << "INFO | " << stats << std::endl;
475 CALL_SUBTEST( check_sparse_solving_real_cases(solver, halfA, b, A, refX) );
476 }
477 else
478 {
479 std::cout << "INFO | Skip sparse problem \"" << it.matname() << "\" (too large)" << std::endl;
480 }
481 }
482 }
483 }
484 #else
485 EIGEN_UNUSED_VARIABLE(maxRealWorldSize);
486 #endif
487 }
488
check_sparse_spd_determinant(Solver & solver)489 template<typename Solver> void check_sparse_spd_determinant(Solver& solver)
490 {
491 typedef typename Solver::MatrixType Mat;
492 typedef typename Mat::Scalar Scalar;
493 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
494
495 // generate the problem
496 Mat A, halfA;
497 DenseMatrix dA;
498 generate_sparse_spd_problem(solver, A, halfA, dA, 30);
499
500 for (int i = 0; i < g_repeat; i++) {
501 check_sparse_determinant(solver, A, dA);
502 check_sparse_determinant(solver, halfA, dA );
503 }
504 }
505
506 template<typename Solver, typename DenseMat>
507 Index generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300, int options = ForceNonZeroDiag)
508 {
509 typedef typename Solver::MatrixType Mat;
510 typedef typename Mat::Scalar Scalar;
511
512 Index size = internal::random<int>(1,maxSize);
513 double density = (std::max)(8./(size*size), 0.01);
514
515 A.resize(size,size);
516 dA.resize(size,size);
517
518 initSparse<Scalar>(density, dA, A, options);
519
520 return size;
521 }
522
523
524 struct prune_column {
525 Index m_col;
prune_columnprune_column526 prune_column(Index col) : m_col(col) {}
527 template<class Scalar>
operatorprune_column528 bool operator()(Index, Index col, const Scalar&) const {
529 return col != m_col;
530 }
531 };
532
533
534 template<typename Solver> void check_sparse_square_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000, bool checkDeficient = false)
535 {
536 typedef typename Solver::MatrixType Mat;
537 typedef typename Mat::Scalar Scalar;
538 typedef SparseMatrix<Scalar,ColMajor, typename Mat::StorageIndex> SpMat;
539 typedef SparseVector<Scalar, 0, typename Mat::StorageIndex> SpVec;
540 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
541 typedef Matrix<Scalar,Dynamic,1> DenseVector;
542
543 int rhsCols = internal::random<int>(1,16);
544
545 Mat A;
546 DenseMatrix dA;
547 for (int i = 0; i < g_repeat; i++) {
548 Index size = generate_sparse_square_problem(solver, A, dA, maxSize);
549
550 A.makeCompressed();
551 DenseVector b = DenseVector::Random(size);
552 DenseMatrix dB(size,rhsCols);
553 SpMat B(size,rhsCols);
554 double density = (std::max)(8./(size*rhsCols), 0.1);
555 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
556 B.makeCompressed();
557 SpVec c = B.col(0);
558 DenseVector dc = dB.col(0);
559 CALL_SUBTEST(check_sparse_solving(solver, A, b, dA, b));
560 CALL_SUBTEST(check_sparse_solving(solver, A, dB, dA, dB));
561 CALL_SUBTEST(check_sparse_solving(solver, A, B, dA, dB));
562 CALL_SUBTEST(check_sparse_solving(solver, A, c, dA, dc));
563
564 // check only once
565 if(i==0)
566 {
567 CALL_SUBTEST(b = DenseVector::Zero(size); check_sparse_solving(solver, A, b, dA, b));
568 }
569 // regression test for Bug 792 (structurally rank deficient matrices):
570 if(checkDeficient && size>1) {
571 Index col = internal::random<int>(0,int(size-1));
572 A.prune(prune_column(col));
573 solver.compute(A);
574 VERIFY_IS_EQUAL(solver.info(), NumericalIssue);
575 }
576 }
577
578 // First, get the folder
579 #ifdef TEST_REAL_CASES
580 // Test real problems with double precision only
581 if (internal::is_same<typename NumTraits<Scalar>::Real, double>::value)
582 {
583 std::string mat_folder = get_matrixfolder<Scalar>();
584 MatrixMarketIterator<Scalar> it(mat_folder);
585 for (; it; ++it)
586 {
587 A = it.matrix();
588 if(A.diagonal().size() <= maxRealWorldSize)
589 {
590 DenseVector b = it.rhs();
591 DenseVector refX = it.refX();
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;
594 CALL_SUBTEST(check_sparse_solving_real_cases(solver, A, b, A, refX));
595 std::string stats = solver_stats(solver);
596 if(stats.size()>0)
597 std::cout << "INFO | " << stats << std::endl;
598 }
599 else
600 {
601 std::cout << "INFO | SKIP sparse problem \"" << it.matname() << "\" (too large)" << std::endl;
602 }
603 }
604 }
605 #else
606 EIGEN_UNUSED_VARIABLE(maxRealWorldSize);
607 #endif
608
609 }
610
check_sparse_square_determinant(Solver & solver)611 template<typename Solver> void check_sparse_square_determinant(Solver& solver)
612 {
613 typedef typename Solver::MatrixType Mat;
614 typedef typename Mat::Scalar Scalar;
615 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
616
617 for (int i = 0; i < g_repeat; i++) {
618 // generate the problem
619 Mat A;
620 DenseMatrix dA;
621
622 int size = internal::random<int>(1,30);
623 dA.setRandom(size,size);
624
625 dA = (dA.array().abs()<0.3).select(0,dA);
626 dA.diagonal() = (dA.diagonal().array()==0).select(1,dA.diagonal());
627 A = dA.sparseView();
628 A.makeCompressed();
629
630 check_sparse_determinant(solver, A, dA);
631 }
632 }
633
check_sparse_square_abs_determinant(Solver & solver)634 template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver)
635 {
636 typedef typename Solver::MatrixType Mat;
637 typedef typename Mat::Scalar Scalar;
638 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
639
640 for (int i = 0; i < g_repeat; i++) {
641 // generate the problem
642 Mat A;
643 DenseMatrix dA;
644 generate_sparse_square_problem(solver, A, dA, 30);
645 A.makeCompressed();
646 check_sparse_abs_determinant(solver, A, dA);
647 }
648 }
649
650 template<typename Solver, typename DenseMat>
651 void generate_sparse_leastsquare_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300, int options = ForceNonZeroDiag)
652 {
653 typedef typename Solver::MatrixType Mat;
654 typedef typename Mat::Scalar Scalar;
655
656 int rows = internal::random<int>(1,maxSize);
657 int cols = internal::random<int>(1,rows);
658 double density = (std::max)(8./(rows*cols), 0.01);
659
660 A.resize(rows,cols);
661 dA.resize(rows,cols);
662
663 initSparse<Scalar>(density, dA, A, options);
664 }
665
check_sparse_leastsquare_solving(Solver & solver)666 template<typename Solver> void check_sparse_leastsquare_solving(Solver& solver)
667 {
668 typedef typename Solver::MatrixType Mat;
669 typedef typename Mat::Scalar Scalar;
670 typedef SparseMatrix<Scalar,ColMajor, typename Mat::StorageIndex> SpMat;
671 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
672 typedef Matrix<Scalar,Dynamic,1> DenseVector;
673
674 int rhsCols = internal::random<int>(1,16);
675
676 Mat A;
677 DenseMatrix dA;
678 for (int i = 0; i < g_repeat; i++) {
679 generate_sparse_leastsquare_problem(solver, A, dA);
680
681 A.makeCompressed();
682 DenseVector b = DenseVector::Random(A.rows());
683 DenseMatrix dB(A.rows(),rhsCols);
684 SpMat B(A.rows(),rhsCols);
685 double density = (std::max)(8./(A.rows()*rhsCols), 0.1);
686 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
687 B.makeCompressed();
688 check_sparse_solving(solver, A, b, dA, b);
689 check_sparse_solving(solver, A, dB, dA, dB);
690 check_sparse_solving(solver, A, B, dA, dB);
691
692 // check only once
693 if(i==0)
694 {
695 b = DenseVector::Zero(A.rows());
696 check_sparse_solving(solver, A, b, dA, b);
697 }
698 }
699 }
700