sparse_solver.h
Go to the documentation of this file.
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>
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>
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>
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>
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>
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>
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>
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>
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;
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>
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 }
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>
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>
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;
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();
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 
489 template<typename Solver> void check_sparse_spd_determinant(Solver& solver)
490 {
491  typedef typename Solver::MatrixType Mat;
492  typedef typename Mat::Scalar Scalar;
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 {
526  prune_column(Index col) : m_col(col) {}
527  template<class Scalar>
528  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;
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 
611 template<typename Solver> void check_sparse_square_determinant(Solver& solver)
612 {
613  typedef typename Solver::MatrixType Mat;
614  typedef typename Mat::Scalar Scalar;
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 
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;
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 
666 template<typename Solver> void check_sparse_leastsquare_solving(Solver& solver)
667 {
668  typedef typename Solver::MatrixType Mat;
669  typedef typename Mat::Scalar Scalar;
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++) {
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 }
Matrix< SCALARB, Dynamic, Dynamic, opt_B > B
Definition: bench_gemm.cpp:49
Matrix< Scalar, Dynamic, Dynamic > DenseMatrix
Index cols() const
Definition: SparseMatrix.h:140
SCALAR Scalar
Definition: bench_gemm.cpp:46
internal::traits< Derived >::Scalar Scalar
#define max(a, b)
Definition: datatypes.h:20
Matrix< RealScalar, Dynamic, Dynamic > M
Definition: bench_gemm.cpp:51
void check_sparse_solving(Solver &solver, const typename Solver::MatrixType &A, const Rhs &b, const DenseMat &dA, const DenseRhs &db)
Definition: sparse_solver.h:43
Scalar * b
Definition: benchVecAdd.cpp:17
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:168
Index rows() const
Definition: SparseMatrix.h:138
#define min(a, b)
Definition: datatypes.h:19
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
A base class for sparse solvers.
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
MatrixXf MatrixType
void check_sparse_determinant(Solver &solver, const typename Solver::MatrixType &A, const DenseMat &dA)
BiCGSTAB< SparseMatrix< double > > solver
Pose3 x2(Rot3::Ypr(0.0, 0.0, 0.0), l2)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
void check_sparse_spd_determinant(Solver &solver)
Matrix< SCALARA, Dynamic, Dynamic, opt_A > A
Definition: bench_gemm.cpp:48
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
bool stats
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:17
void g(const string &key, int i)
Definition: testBTree.cpp:41
bool operator()(Index, Index col, const Scalar &) const
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
Matrix< Scalar, Dynamic, 1 > DenseVector
#define VERIFY_IS_APPROX(a, b)
void check_sparse_leastsquare_solving(Solver &solver)
Base class of any sparse matrices or sparse expressions.
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
#define VERIFY_IS_EQUAL(a, b)
Definition: main.h:386
prune_column(Index col)
void check_sparse_solving_real_cases(Solver &solver, const typename Solver::MatrixType &A, const Rhs &b, const typename Solver::MatrixType &fullA, const Rhs &refX)
a sparse vector class
Definition: SparseUtil.h:54
static int g_repeat
Definition: main.h:169
idx_t idx_t idx_t idx_t idx_t idx_t idx_t real_t real_t idx_t * options
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Iterator to browse matrices from a specified folder.
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Pose3 x3(Rot3::Ypr(M_PI/4.0, 0.0, 0.0), l2)
static std::stringstream ss
Definition: testBTree.cpp:31
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:159
Eigen::SparseMatrix< double > SpMat
int generate_sparse_spd_problem(Solver &, typename Solver::MatrixType &A, typename Solver::MatrixType &halfA, DenseMat &dA, int maxSize=300)
void check_sparse_abs_determinant(Solver &solver, const typename Solver::MatrixType &A, const DenseMat &dA)
#define CALL_SUBTEST(FUNC)
Definition: main.h:399
#define VERIFY(a)
Definition: main.h:380
#define EIGEN_TEST_MAX_SIZE
Pose3 x1
Definition: testPose3.cpp:663
Index outerSize() const
Definition: SparseMatrix.h:145
m col(1)
Matrix< Scalar, Dynamic, Dynamic > Mat
Definition: gemm_common.h:15
Index generate_sparse_square_problem(Solver &, typename Solver::MatrixType &A, DenseMat &dA, int maxSize=300, int options=ForceNonZeroDiag)
void generate_sparse_leastsquare_problem(Solver &, typename Solver::MatrixType &A, DenseMat &dA, int maxSize=300, int options=ForceNonZeroDiag)
void solve_with_guess(IterativeSolverBase< Solver > &solver, const MatrixBase< Rhs > &b, const Guess &g, Result &x)
Definition: sparse_solver.h:16
static int g_test_level
Definition: main.h:168
internal::nested_eval< T, 1 >::type eval(const T &xpr)
void check_sparse_square_determinant(Solver &solver)
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
#define abs(x)
Definition: datatypes.h:17
Base class for linear iterative solvers.
NumTraits< Scalar >::Real RealScalar
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Scalar * valuePtr() const
Definition: SparseMatrix.h:150
void check_sparse_square_abs_determinant(Solver &solver)
void check_sparse_square_solving(Solver &solver, int maxSize=300, int maxRealWorldSize=100000, bool checkDeficient=false)
#define EIGEN_UNUSED_VARIABLE(var)
Definition: Macros.h:1076
void check_sparse_spd_solving(Solver &solver, int maxSize=(std::min)(300, EIGEN_TEST_MAX_SIZE), int maxRealWorldSize=100000)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:54