sparse_solver.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd
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 13template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs> 14void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db) 15{ 16 typedef typename Solver::MatrixType Mat; 17 typedef typename Mat::Scalar Scalar; 18 19 DenseRhs refX = dA.lu().solve(db); 20 21 Rhs x(b.rows(), b.cols()); 22 Rhs oldb = b; 23 24 solver.compute(A); 25 if (solver.info() != Success) 26 { 27 std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n"; 28 exit(0); 29 return; 30 } 31 x = solver.solve(b); 32 if (solver.info() != Success) 33 { 34 std::cerr << "sparse solver testing: solving failed\n"; 35 return; 36 } 37 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); 38 39 VERIFY(x.isApprox(refX,test_precision<Scalar>())); 40 41 x.setZero(); 42 // test the analyze/factorize API 43 solver.analyzePattern(A); 44 solver.factorize(A); 45 if (solver.info() != Success) 46 { 47 std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n"; 48 exit(0); 49 return; 50 } 51 x = solver.solve(b); 52 if (solver.info() != Success) 53 { 54 std::cerr << "sparse solver testing: solving failed\n"; 55 return; 56 } 57 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); 58 59 VERIFY(x.isApprox(refX,test_precision<Scalar>())); 60 61 // test Block as the result and rhs: 62 { 63 DenseRhs x(db.rows(), db.cols()); 64 DenseRhs b(db), oldb(db); 65 x.setZero(); 66 x.block(0,0,x.rows(),x.cols()) = solver.solve(b.block(0,0,b.rows(),b.cols())); 67 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); 68 VERIFY(x.isApprox(refX,test_precision<Scalar>())); 69 } 70} 71 72template<typename Solver, typename Rhs> 73void check_sparse_solving_real_cases(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const Rhs& refX) 74{ 75 typedef typename Solver::MatrixType Mat; 76 typedef typename Mat::Scalar Scalar; 77 typedef typename Mat::RealScalar RealScalar; 78 79 Rhs x(b.rows(), b.cols()); 80 81 solver.compute(A); 82 if (solver.info() != Success) 83 { 84 std::cerr << "sparse solver testing: factorization failed (check_sparse_solving_real_cases)\n"; 85 exit(0); 86 return; 87 } 88 x = solver.solve(b); 89 if (solver.info() != Success) 90 { 91 std::cerr << "sparse solver testing: solving failed\n"; 92 return; 93 } 94 95 RealScalar res_error; 96 // Compute the norm of the relative error 97 if(refX.size() != 0) 98 res_error = (refX - x).norm()/refX.norm(); 99 else 100 { 101 // Compute the relative residual norm 102 res_error = (b - A * x).norm()/b.norm(); 103 } 104 if (res_error > test_precision<Scalar>() ){ 105 std::cerr << "Test " << g_test_stack.back() << " failed in "EI_PP_MAKE_STRING(__FILE__) 106 << " (" << EI_PP_MAKE_STRING(__LINE__) << ")" << std::endl << std::endl; 107 abort(); 108 } 109 110} 111template<typename Solver, typename DenseMat> 112void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA) 113{ 114 typedef typename Solver::MatrixType Mat; 115 typedef typename Mat::Scalar Scalar; 116 typedef typename Mat::RealScalar RealScalar; 117 118 solver.compute(A); 119 if (solver.info() != Success) 120 { 121 std::cerr << "sparse solver testing: factorization failed (check_sparse_determinant)\n"; 122 return; 123 } 124 125 Scalar refDet = dA.determinant(); 126 VERIFY_IS_APPROX(refDet,solver.determinant()); 127} 128 129 130template<typename Solver, typename DenseMat> 131int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300) 132{ 133 typedef typename Solver::MatrixType Mat; 134 typedef typename Mat::Scalar Scalar; 135 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 136 137 int size = internal::random<int>(1,maxSize); 138 double density = (std::max)(8./(size*size), 0.01); 139 140 Mat M(size, size); 141 DenseMatrix dM(size, size); 142 143 initSparse<Scalar>(density, dM, M, ForceNonZeroDiag); 144 145 A = M * M.adjoint(); 146 dA = dM * dM.adjoint(); 147 148 halfA.resize(size,size); 149 halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M); 150 151 return size; 152} 153 154 155#ifdef TEST_REAL_CASES 156template<typename Scalar> 157inline std::string get_matrixfolder() 158{ 159 std::string mat_folder = TEST_REAL_CASES; 160 if( internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value ) 161 mat_folder = mat_folder + static_cast<string>("/complex/"); 162 else 163 mat_folder = mat_folder + static_cast<string>("/real/"); 164 return mat_folder; 165} 166#endif 167 168template<typename Solver> void check_sparse_spd_solving(Solver& solver) 169{ 170 typedef typename Solver::MatrixType Mat; 171 typedef typename Mat::Scalar Scalar; 172 typedef typename Mat::Index Index; 173 typedef SparseMatrix<Scalar,ColMajor> SpMat; 174 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 175 typedef Matrix<Scalar,Dynamic,1> DenseVector; 176 177 // generate the problem 178 Mat A, halfA; 179 DenseMatrix dA; 180 int size = generate_sparse_spd_problem(solver, A, halfA, dA); 181 182 // generate the right hand sides 183 int rhsCols = internal::random<int>(1,16); 184 double density = (std::max)(8./(size*rhsCols), 0.1); 185 SpMat B(size,rhsCols); 186 DenseVector b = DenseVector::Random(size); 187 DenseMatrix dB(size,rhsCols); 188 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); 189 190 for (int i = 0; i < g_repeat; i++) { 191 check_sparse_solving(solver, A, b, dA, b); 192 check_sparse_solving(solver, halfA, b, dA, b); 193 check_sparse_solving(solver, A, dB, dA, dB); 194 check_sparse_solving(solver, halfA, dB, dA, dB); 195 check_sparse_solving(solver, A, B, dA, dB); 196 check_sparse_solving(solver, halfA, B, dA, dB); 197 } 198 199 // First, get the folder 200#ifdef TEST_REAL_CASES 201 if (internal::is_same<Scalar, float>::value 202 || internal::is_same<Scalar, std::complex<float> >::value) 203 return ; 204 205 std::string mat_folder = get_matrixfolder<Scalar>(); 206 MatrixMarketIterator<Scalar> it(mat_folder); 207 for (; it; ++it) 208 { 209 if (it.sym() == SPD){ 210 Mat halfA; 211 PermutationMatrix<Dynamic, Dynamic, Index> pnull; 212 halfA.template selfadjointView<Solver::UpLo>() = it.matrix().template triangularView<Eigen::Lower>().twistedBy(pnull); 213 214 std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n"; 215 check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX()); 216 check_sparse_solving_real_cases(solver, halfA, it.rhs(), it.refX()); 217 } 218 } 219#endif 220} 221 222template<typename Solver> void check_sparse_spd_determinant(Solver& solver) 223{ 224 typedef typename Solver::MatrixType Mat; 225 typedef typename Mat::Scalar Scalar; 226 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 227 228 // generate the problem 229 Mat A, halfA; 230 DenseMatrix dA; 231 generate_sparse_spd_problem(solver, A, halfA, dA, 30); 232 233 for (int i = 0; i < g_repeat; i++) { 234 check_sparse_determinant(solver, A, dA); 235 check_sparse_determinant(solver, halfA, dA ); 236 } 237} 238 239template<typename Solver, typename DenseMat> 240int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300) 241{ 242 typedef typename Solver::MatrixType Mat; 243 typedef typename Mat::Scalar Scalar; 244 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 245 246 int size = internal::random<int>(1,maxSize); 247 double density = (std::max)(8./(size*size), 0.01); 248 249 A.resize(size,size); 250 dA.resize(size,size); 251 252 initSparse<Scalar>(density, dA, A, ForceNonZeroDiag); 253 254 return size; 255} 256 257template<typename Solver> void check_sparse_square_solving(Solver& solver) 258{ 259 typedef typename Solver::MatrixType Mat; 260 typedef typename Mat::Scalar Scalar; 261 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 262 typedef Matrix<Scalar,Dynamic,1> DenseVector; 263 264 int rhsCols = internal::random<int>(1,16); 265 266 Mat A; 267 DenseMatrix dA; 268 int size = generate_sparse_square_problem(solver, A, dA); 269 270 DenseVector b = DenseVector::Random(size); 271 DenseMatrix dB = DenseMatrix::Random(size,rhsCols); 272 A.makeCompressed(); 273 for (int i = 0; i < g_repeat; i++) { 274 check_sparse_solving(solver, A, b, dA, b); 275 check_sparse_solving(solver, A, dB, dA, dB); 276 } 277 278 // First, get the folder 279#ifdef TEST_REAL_CASES 280 if (internal::is_same<Scalar, float>::value 281 || internal::is_same<Scalar, std::complex<float> >::value) 282 return ; 283 284 std::string mat_folder = get_matrixfolder<Scalar>(); 285 MatrixMarketIterator<Scalar> it(mat_folder); 286 for (; it; ++it) 287 { 288 std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n"; 289 check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX()); 290 } 291#endif 292 293} 294 295template<typename Solver> void check_sparse_square_determinant(Solver& solver) 296{ 297 typedef typename Solver::MatrixType Mat; 298 typedef typename Mat::Scalar Scalar; 299 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 300 301 // generate the problem 302 Mat A; 303 DenseMatrix dA; 304 generate_sparse_square_problem(solver, A, dA, 30); 305 A.makeCompressed(); 306 for (int i = 0; i < g_repeat; i++) { 307 check_sparse_determinant(solver, A, dA); 308 } 309} 310