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