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 halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M); 165 166 return size; 167 } 168 169 170 #ifdef TEST_REAL_CASES 171 template<typename Scalar> 172 inline std::string get_matrixfolder() 173 { 174 std::string mat_folder = TEST_REAL_CASES; 175 if( internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value ) 176 mat_folder = mat_folder + static_cast<std::string>("/complex/"); 177 else 178 mat_folder = mat_folder + static_cast<std::string>("/real/"); 179 return mat_folder; 180 } 181 #endif 182 183 template<typename Solver> void check_sparse_spd_solving(Solver& solver) 184 { 185 typedef typename Solver::MatrixType Mat; 186 typedef typename Mat::Scalar Scalar; 187 typedef SparseMatrix<Scalar,ColMajor> SpMat; 188 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 189 typedef Matrix<Scalar,Dynamic,1> DenseVector; 190 191 // generate the problem 192 Mat A, halfA; 193 DenseMatrix dA; 194 for (int i = 0; i < g_repeat; i++) { 195 int size = generate_sparse_spd_problem(solver, A, halfA, dA); 196 197 // generate the right hand sides 198 int rhsCols = internal::random<int>(1,16); 199 double density = (std::max)(8./(size*rhsCols), 0.1); 200 SpMat B(size,rhsCols); 201 DenseVector b = DenseVector::Random(size); 202 DenseMatrix dB(size,rhsCols); 203 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); 204 205 check_sparse_solving(solver, A, b, dA, b); 206 check_sparse_solving(solver, halfA, b, dA, b); 207 check_sparse_solving(solver, A, dB, dA, dB); 208 check_sparse_solving(solver, halfA, dB, dA, dB); 209 check_sparse_solving(solver, A, B, dA, dB); 210 check_sparse_solving(solver, halfA, B, dA, dB); 211 212 // check only once 213 if(i==0) 214 { 215 b = DenseVector::Zero(size); 216 check_sparse_solving(solver, A, b, dA, b); 217 } 218 } 219 220 // First, get the folder 221 #ifdef TEST_REAL_CASES 222 if (internal::is_same<Scalar, float>::value 223 || internal::is_same<Scalar, std::complex<float> >::value) 224 return ; 225 226 std::string mat_folder = get_matrixfolder<Scalar>(); 227 MatrixMarketIterator<Scalar> it(mat_folder); 228 for (; it; ++it) 229 { 230 if (it.sym() == SPD){ 231 Mat halfA; 232 PermutationMatrix<Dynamic, Dynamic, Index> pnull; 233 halfA.template selfadjointView<Solver::UpLo>() = it.matrix().template triangularView<Eigen::Lower>().twistedBy(pnull); 234 235 std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n"; 236 check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX()); 237 check_sparse_solving_real_cases(solver, halfA, it.rhs(), it.refX()); 238 } 239 } 240 #endif 241 } 242 243 template<typename Solver> void check_sparse_spd_determinant(Solver& solver) 244 { 245 typedef typename Solver::MatrixType Mat; 246 typedef typename Mat::Scalar Scalar; 247 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 248 249 // generate the problem 250 Mat A, halfA; 251 DenseMatrix dA; 252 generate_sparse_spd_problem(solver, A, halfA, dA, 30); 253 254 for (int i = 0; i < g_repeat; i++) { 255 check_sparse_determinant(solver, A, dA); 256 check_sparse_determinant(solver, halfA, dA ); 257 } 258 } 259 260 template<typename Solver, typename DenseMat> 261 int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300) 262 { 263 typedef typename Solver::MatrixType Mat; 264 typedef typename Mat::Scalar Scalar; 265 266 int size = internal::random<int>(1,maxSize); 267 double density = (std::max)(8./(size*size), 0.01); 268 269 A.resize(size,size); 270 dA.resize(size,size); 271 272 initSparse<Scalar>(density, dA, A, ForceNonZeroDiag); 273 274 return size; 275 } 276 277 template<typename Solver> void check_sparse_square_solving(Solver& solver) 278 { 279 typedef typename Solver::MatrixType Mat; 280 typedef typename Mat::Scalar Scalar; 281 typedef SparseMatrix<Scalar,ColMajor> SpMat; 282 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 283 typedef Matrix<Scalar,Dynamic,1> DenseVector; 284 285 int rhsCols = internal::random<int>(1,16); 286 287 Mat A; 288 DenseMatrix dA; 289 for (int i = 0; i < g_repeat; i++) { 290 int size = generate_sparse_square_problem(solver, A, dA); 291 292 A.makeCompressed(); 293 DenseVector b = DenseVector::Random(size); 294 DenseMatrix dB(size,rhsCols); 295 SpMat B(size,rhsCols); 296 double density = (std::max)(8./(size*rhsCols), 0.1); 297 initSparse<Scalar>(density, dB, B, ForceNonZeroDiag); 298 B.makeCompressed(); 299 check_sparse_solving(solver, A, b, dA, b); 300 check_sparse_solving(solver, A, dB, dA, dB); 301 check_sparse_solving(solver, A, B, dA, dB); 302 303 // check only once 304 if(i==0) 305 { 306 b = DenseVector::Zero(size); 307 check_sparse_solving(solver, A, b, dA, b); 308 } 309 } 310 311 // First, get the folder 312 #ifdef TEST_REAL_CASES 313 if (internal::is_same<Scalar, float>::value 314 || internal::is_same<Scalar, std::complex<float> >::value) 315 return ; 316 317 std::string mat_folder = get_matrixfolder<Scalar>(); 318 MatrixMarketIterator<Scalar> it(mat_folder); 319 for (; it; ++it) 320 { 321 std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n"; 322 check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX()); 323 } 324 #endif 325 326 } 327 328 template<typename Solver> void check_sparse_square_determinant(Solver& solver) 329 { 330 typedef typename Solver::MatrixType Mat; 331 typedef typename Mat::Scalar Scalar; 332 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 333 334 // generate the problem 335 Mat A; 336 DenseMatrix dA; 337 generate_sparse_square_problem(solver, A, dA, 30); 338 A.makeCompressed(); 339 for (int i = 0; i < g_repeat; i++) { 340 check_sparse_determinant(solver, A, dA); 341 } 342 } 343 344 template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver) 345 { 346 typedef typename Solver::MatrixType Mat; 347 typedef typename Mat::Scalar Scalar; 348 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 349 350 // generate the problem 351 Mat A; 352 DenseMatrix dA; 353 generate_sparse_square_problem(solver, A, dA, 30); 354 A.makeCompressed(); 355 for (int i = 0; i < g_repeat; i++) { 356 check_sparse_abs_determinant(solver, A, dA); 357 } 358 } 359 360
