1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud (at) inria.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 12 template<typename Scalar> void 13 initSPD(double density, 14 Matrix<Scalar,Dynamic,Dynamic>& refMat, 15 SparseMatrix<Scalar>& sparseMat) 16 { 17 Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols()); 18 initSparse(density,refMat,sparseMat); 19 refMat = refMat * refMat.adjoint(); 20 for (int k=0; k<2; ++k) 21 { 22 initSparse(density,aux,sparseMat,ForceNonZeroDiag); 23 refMat += aux * aux.adjoint(); 24 } 25 sparseMat.setZero(); 26 for (int j=0 ; j<sparseMat.cols(); ++j) 27 for (int i=j ; i<sparseMat.rows(); ++i) 28 if (refMat(i,j)!=Scalar(0)) 29 sparseMat.insert(i,j) = refMat(i,j); 30 sparseMat.finalize(); 31 } 32 33 template<typename Scalar> void sparse_solvers(int rows, int cols) 34 { 35 double density = (std::max)(8./(rows*cols), 0.01); 36 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 37 typedef Matrix<Scalar,Dynamic,1> DenseVector; 38 // Scalar eps = 1e-6; 39 40 DenseVector vec1 = DenseVector::Random(rows); 41 42 std::vector<Vector2i> zeroCoords; 43 std::vector<Vector2i> nonzeroCoords; 44 45 // test triangular solver 46 { 47 DenseVector vec2 = vec1, vec3 = vec1; 48 SparseMatrix<Scalar> m2(rows, cols); 49 DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols); 50 51 // lower - dense 52 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); 53 VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2), 54 m2.template triangularView<Lower>().solve(vec3)); 55 56 // upper - dense 57 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); 58 VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2), 59 m2.template triangularView<Upper>().solve(vec3)); 60 VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2), 61 m2.conjugate().template triangularView<Upper>().solve(vec3)); 62 { 63 SparseMatrix<Scalar> cm2(m2); 64 //Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr 65 MappedSparseMatrix<Scalar> mm2(rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(), cm2.valuePtr()); 66 VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2), 67 mm2.conjugate().template triangularView<Upper>().solve(vec3)); 68 } 69 70 // lower - transpose 71 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); 72 VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2), 73 m2.transpose().template triangularView<Upper>().solve(vec3)); 74 75 // upper - transpose 76 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords); 77 VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2), 78 m2.transpose().template triangularView<Lower>().solve(vec3)); 79 80 SparseMatrix<Scalar> matB(rows, rows); 81 DenseMatrix refMatB = DenseMatrix::Zero(rows, rows); 82 83 // lower - sparse 84 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular); 85 initSparse<Scalar>(density, refMatB, matB); 86 refMat2.template triangularView<Lower>().solveInPlace(refMatB); 87 m2.template triangularView<Lower>().solveInPlace(matB); 88 VERIFY_IS_APPROX(matB.toDense(), refMatB); 89 90 // upper - sparse 91 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular); 92 initSparse<Scalar>(density, refMatB, matB); 93 refMat2.template triangularView<Upper>().solveInPlace(refMatB); 94 m2.template triangularView<Upper>().solveInPlace(matB); 95 VERIFY_IS_APPROX(matB, refMatB); 96 97 // test deprecated API 98 initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords); 99 VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2), 100 m2.template triangularView<Lower>().solve(vec3)); 101 } 102 } 103 104 void test_sparse_solvers() 105 { 106 for(int i = 0; i < g_repeat; i++) { 107 CALL_SUBTEST_1(sparse_solvers<double>(8, 8) ); 108 int s = internal::random<int>(1,300); 109 CALL_SUBTEST_2(sparse_solvers<std::complex<double> >(s,s) ); 110 CALL_SUBTEST_1(sparse_solvers<double>(s,s) ); 111 } 112 } 113