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