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 (at) 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 13 template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs> 14 void 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 x.setZero(); 41 // test the analyze/factorize API 42 solver.analyzePattern(A); 43 solver.factorize(A); 44 if (solver.info() != Success) 45 { 46 std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n"; 47 exit(0); 48 return; 49 } 50 x = solver.solve(b); 51 if (solver.info() != Success) 52 { 53 std::cerr << "sparse solver testing: solving failed\n"; 54 return; 55 } 56 VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!"); 57 58 VERIFY(x.isApprox(refX,test_precision<Scalar>())); 59 } 60 61 // test dense Block as the result and rhs: 62 { 63 DenseRhs x(db.rows(), db.cols()); 64 DenseRhs oldb(db); 65 x.setZero(); 66 x.block(0,0,x.rows(),x.cols()) = solver.solve(db.block(0,0,db.rows(),db.cols())); 67 VERIFY(oldb.isApprox(db) && "sparse solver testing: the rhs should not be modified!"); 68 VERIFY(x.isApprox(refX,test_precision<Scalar>())); 69 } 70 } 71 72 template<typename Solver, typename Rhs> 73 void 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 } 111 template<typename Solver, typename DenseMat> 112 void 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 117 solver.compute(A); 118 if (solver.info() != Success) 119 { 120 std::cerr << "sparse solver testing: factorization failed (check_sparse_determinant)\n"; 121 return; 122 } 123 124 Scalar refDet = dA.determinant(); 125 VERIFY_IS_APPROX(refDet,solver.determinant()); 126 } 127 template<typename Solver, typename DenseMat> 128 void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA) 129 { 130 using std::abs; 131 typedef typename Solver::MatrixType Mat; 132 typedef typename Mat::Scalar Scalar; 133 134 solver.compute(A); 135 if (solver.info() != Success) 136 { 137 std::cerr << "sparse solver testing: factorization failed (check_sparse_abs_determinant)\n"; 138 return; 139 } 140 141 Scalar refDet = abs(dA.determinant()); 142 VERIFY_IS_APPROX(refDet,solver.absDeterminant()); 143 } 144 145 template<typename Solver, typename DenseMat> 146 int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300) 147 { 148 typedef typename Solver::MatrixType Mat; 149 typedef typename Mat::Scalar Scalar; 150 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 151 152 int size = internal::random<int>(1,maxSize); 153 double density = (std::max)(8./(size*size), 0.01); 154 155 Mat M(size, size); 156 DenseMatrix dM(size, size); 157 158 initSparse<Scalar>(density, dM, M, ForceNonZeroDiag); 159 160 A = M * M.adjoint(); 161 dA = dM * dM.adjoint(); 162 163 halfA.resize(size,size); 164 if(Solver::UpLo==(Lower|Upper)) 165 halfA = A; 166 else 167 halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M); 168 169 return size; 170 } 171 172 173 #ifdef TEST_REAL_CASES 174 template<typename Scalar> 175 inline std::string get_matrixfolder() 176 { 177 std::string mat_folder = TEST_REAL_CASES; 178 if( internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value ) 179 mat_folder = mat_folder + static_cast<std::string>("/complex/"); 180 else 181 mat_folder = mat_folder + static_cast<std::string>("/real/"); 182 return mat_folder; 183 } 184 #endif 185 186 template<typename Solver> void check_sparse_spd_solving(Solver& solver) 187 { 188 typedef typename Solver::MatrixType Mat; 189 typedef typename Mat::Scalar Scalar; 190 typedef SparseMatrix<Scalar,ColMajor> SpMat; 191 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 192 typedef Matrix<Scalar,Dynamic,1> DenseVector; 193 194 // generate the problem 195 Mat A, halfA; 196 DenseMatrix dA; 197 for (int i = 0; i < g_repeat; i++) { 198 int size = generate_sparse_spd_problem(solver, A, halfA, dA); 199 200 // generate the right hand sides 201 int rhsCols = internal::random<int>(1,16); 202 double density = (std::max)(8./(size*rhsCols), 0.1); 203 SpMat B(size,rhsCols); 204 DenseVector b = DenseVector::Random(size); 205 DenseMatrix dB(size,rhsCols); 206 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); 207 208 check_sparse_solving(solver, A, b, dA, b); 209 check_sparse_solving(solver, halfA, b, dA, b); 210 check_sparse_solving(solver, A, dB, dA, dB); 211 check_sparse_solving(solver, halfA, dB, dA, dB); 212 check_sparse_solving(solver, A, B, dA, dB); 213 check_sparse_solving(solver, halfA, B, dA, dB); 214 215 // check only once 216 if(i==0) 217 { 218 b = DenseVector::Zero(size); 219 check_sparse_solving(solver, A, b, dA, b); 220 } 221 } 222 223 // First, get the folder 224 #ifdef TEST_REAL_CASES 225 if (internal::is_same<Scalar, float>::value 226 || internal::is_same<Scalar, std::complex<float> >::value) 227 return ; 228 229 std::string mat_folder = get_matrixfolder<Scalar>(); 230 MatrixMarketIterator<Scalar> it(mat_folder); 231 for (; it; ++it) 232 { 233 if (it.sym() == SPD){ 234 Mat halfA; 235 PermutationMatrix<Dynamic, Dynamic, Index> pnull; 236 halfA.template selfadjointView<Solver::UpLo>() = it.matrix().template triangularView<Eigen::Lower>().twistedBy(pnull); 237 238 std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n"; 239 check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX()); 240 check_sparse_solving_real_cases(solver, halfA, it.rhs(), it.refX()); 241 } 242 } 243 #endif 244 } 245 246 template<typename Solver> void check_sparse_spd_determinant(Solver& solver) 247 { 248 typedef typename Solver::MatrixType Mat; 249 typedef typename Mat::Scalar Scalar; 250 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 251 252 // generate the problem 253 Mat A, halfA; 254 DenseMatrix dA; 255 generate_sparse_spd_problem(solver, A, halfA, dA, 30); 256 257 for (int i = 0; i < g_repeat; i++) { 258 check_sparse_determinant(solver, A, dA); 259 check_sparse_determinant(solver, halfA, dA ); 260 } 261 } 262 263 template<typename Solver, typename DenseMat> 264 int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300) 265 { 266 typedef typename Solver::MatrixType Mat; 267 typedef typename Mat::Scalar Scalar; 268 269 int size = internal::random<int>(1,maxSize); 270 double density = (std::max)(8./(size*size), 0.01); 271 272 A.resize(size,size); 273 dA.resize(size,size); 274 275 initSparse<Scalar>(density, dA, A, ForceNonZeroDiag); 276 277 return size; 278 } 279 280 template<typename Solver> void check_sparse_square_solving(Solver& solver) 281 { 282 typedef typename Solver::MatrixType Mat; 283 typedef typename Mat::Scalar Scalar; 284 typedef SparseMatrix<Scalar,ColMajor> SpMat; 285 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 286 typedef Matrix<Scalar,Dynamic,1> DenseVector; 287 288 int rhsCols = internal::random<int>(1,16); 289 290 Mat A; 291 DenseMatrix dA; 292 for (int i = 0; i < g_repeat; i++) { 293 int size = generate_sparse_square_problem(solver, A, dA); 294 295 A.makeCompressed(); 296 DenseVector b = DenseVector::Random(size); 297 DenseMatrix dB(size,rhsCols); 298 SpMat B(size,rhsCols); 299 double density = (std::max)(8./(size*rhsCols), 0.1); 300 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); 301 B.makeCompressed(); 302 check_sparse_solving(solver, A, b, dA, b); 303 check_sparse_solving(solver, A, dB, dA, dB); 304 check_sparse_solving(solver, A, B, dA, dB); 305 306 // check only once 307 if(i==0) 308 { 309 b = DenseVector::Zero(size); 310 check_sparse_solving(solver, A, b, dA, b); 311 } 312 } 313 314 // First, get the folder 315 #ifdef TEST_REAL_CASES 316 if (internal::is_same<Scalar, float>::value 317 || internal::is_same<Scalar, std::complex<float> >::value) 318 return ; 319 320 std::string mat_folder = get_matrixfolder<Scalar>(); 321 MatrixMarketIterator<Scalar> it(mat_folder); 322 for (; it; ++it) 323 { 324 std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n"; 325 check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX()); 326 } 327 #endif 328 329 } 330 331 template<typename Solver> void check_sparse_square_determinant(Solver& solver) 332 { 333 typedef typename Solver::MatrixType Mat; 334 typedef typename Mat::Scalar Scalar; 335 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 336 337 // generate the problem 338 Mat A; 339 DenseMatrix dA; 340 generate_sparse_square_problem(solver, A, dA, 30); 341 A.makeCompressed(); 342 for (int i = 0; i < g_repeat; i++) { 343 check_sparse_determinant(solver, A, dA); 344 } 345 } 346 347 template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver) 348 { 349 typedef typename Solver::MatrixType Mat; 350 typedef typename Mat::Scalar Scalar; 351 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 352 353 // generate the problem 354 Mat A; 355 DenseMatrix dA; 356 generate_sparse_square_problem(solver, A, dA, 30); 357 A.makeCompressed(); 358 for (int i = 0; i < g_repeat; i++) { 359 check_sparse_abs_determinant(solver, A, dA); 360 } 361 } 362 363
