sparse_solver.h (24396B)
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; 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> 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; 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 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; 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; 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 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 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 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 }