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 14template<typename Solver, typename Rhs, typename Guess,typename Result> 15void 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 28template<typename Solver, typename Rhs, typename Guess,typename Result> 29void 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 36template<typename Solver, typename Rhs, typename Guess,typename Result> 37void solve_with_guess(SparseSolverBase<Solver>& solver, const SparseMatrixBase<Rhs>& b, const Guess& , Result& x) { 38 x = solver.derived().solve(b); 39} 40 41template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs> 42void 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 143template<typename Solver, typename Rhs> 144void 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} 176template<typename Solver, typename DenseMat> 177void 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} 192template<typename Solver, typename DenseMat> 193void 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 210template<typename Solver, typename DenseMat> 211int 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 239template<typename Scalar> 240inline 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} 249std::string sym_to_string(int sym) 250{ 251 if(sym==Symmetric) return "Symmetric "; 252 if(sym==SPD) return "SPD "; 253 return ""; 254} 255template<typename Derived> 256std::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} 262template<typename Derived> 263std::string solver_stats(const SparseSolverBase<Derived> &/*solver*/) 264{ 265 return ""; 266} 267#endif 268 269template<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 354template<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 371template<typename Solver, typename DenseMat> 372Index 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 389struct prune_column { 390 Index m_col; 391 prune_column(Index col) : m_col(col) {} 392 template<class Scalar> 393 bool operator()(Index, Index col, const Scalar&) const { 394 return col != m_col; 395 } 396}; 397 398 399template<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 477template<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 500template<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 516template<typename Solver, typename DenseMat> 517void 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 532template<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