Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | |
| 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); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 20 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 21 | Rhs x(b.rows(), b.cols()); |
| 22 | Rhs oldb = b; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 23 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 60 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 61 | // test dense Block as the result and rhs: |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 62 | { |
| 63 | DenseRhs x(db.rows(), db.cols()); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 64 | DenseRhs oldb(db); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 65 | x.setZero(); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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!"); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 ) |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 160 | mat_folder = mat_folder + static_cast<std::string>("/complex/"); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 161 | else |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 162 | mat_folder = mat_folder + static_cast<std::string>("/real/"); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 178 | for (int i = 0; i < g_repeat; i++) { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 202 | } |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 203 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 265 | typedef SparseMatrix<Scalar,ColMajor> SpMat; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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; |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 273 | for (int i = 0; i < g_repeat; i++) { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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(); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 283 | check_sparse_solving(solver, A, b, dA, b); |
| 284 | check_sparse_solving(solver, A, dB, dA, dB); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 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 | } |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 293 | } |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 294 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | } |