[136] | 1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
| 2 | // for linear algebra. Eigen itself is part of the KDE project.
|
---|
| 3 | //
|
---|
| 4 | // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
|
---|
| 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.startFill();
|
---|
| 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.fill(i,j) = refMat(i,j);
|
---|
| 30 | sparseMat.endFill();
|
---|
| 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
|
---|
| 52 | initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
|
---|
| 53 | VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2),
|
---|
| 54 | m2.template marked<LowerTriangular>().solveTriangular(vec3));
|
---|
| 55 |
|
---|
| 56 | // lower - transpose
|
---|
| 57 | initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
|
---|
| 58 | VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().transpose().solveTriangular(vec2),
|
---|
| 59 | m2.template marked<LowerTriangular>().transpose().solveTriangular(vec3));
|
---|
| 60 |
|
---|
| 61 | // upper
|
---|
| 62 | initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
|
---|
| 63 | VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().solveTriangular(vec2),
|
---|
| 64 | m2.template marked<UpperTriangular>().solveTriangular(vec3));
|
---|
| 65 |
|
---|
| 66 | // upper - transpose
|
---|
| 67 | initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
|
---|
| 68 | VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().transpose().solveTriangular(vec2),
|
---|
| 69 | m2.template marked<UpperTriangular>().transpose().solveTriangular(vec3));
|
---|
| 70 | }
|
---|
| 71 |
|
---|
| 72 | // test LLT
|
---|
| 73 | {
|
---|
| 74 | // TODO fix the issue with complex (see SparseLLT::solveInPlace)
|
---|
| 75 | SparseMatrix<Scalar> m2(rows, cols);
|
---|
| 76 | DenseMatrix refMat2(rows, cols);
|
---|
| 77 |
|
---|
| 78 | DenseVector b = DenseVector::Random(cols);
|
---|
| 79 | DenseVector refX(cols), x(cols);
|
---|
| 80 |
|
---|
| 81 | initSPD(density, refMat2, m2);
|
---|
| 82 |
|
---|
| 83 | refMat2.llt().solve(b, &refX);
|
---|
| 84 | typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
|
---|
| 85 | if (!NumTraits<Scalar>::IsComplex)
|
---|
| 86 | {
|
---|
| 87 | x = b;
|
---|
| 88 | SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
|
---|
| 89 | VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
|
---|
| 90 | }
|
---|
| 91 | #ifdef EIGEN_CHOLMOD_SUPPORT
|
---|
| 92 | x = b;
|
---|
| 93 | SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x);
|
---|
| 94 | VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
|
---|
| 95 | #endif
|
---|
| 96 | if (!NumTraits<Scalar>::IsComplex)
|
---|
| 97 | {
|
---|
| 98 | #ifdef EIGEN_TAUCS_SUPPORT
|
---|
| 99 | x = b;
|
---|
| 100 | SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
|
---|
| 101 | VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
|
---|
| 102 | x = b;
|
---|
| 103 | SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
|
---|
| 104 | VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
|
---|
| 105 | x = b;
|
---|
| 106 | SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
|
---|
| 107 | VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
|
---|
| 108 | #endif
|
---|
| 109 | }
|
---|
| 110 | }
|
---|
| 111 |
|
---|
| 112 | // test LDLT
|
---|
| 113 | if (!NumTraits<Scalar>::IsComplex)
|
---|
| 114 | {
|
---|
| 115 | // TODO fix the issue with complex (see SparseLDLT::solveInPlace)
|
---|
| 116 | SparseMatrix<Scalar> m2(rows, cols);
|
---|
| 117 | DenseMatrix refMat2(rows, cols);
|
---|
| 118 |
|
---|
| 119 | DenseVector b = DenseVector::Random(cols);
|
---|
| 120 | DenseVector refX(cols), x(cols);
|
---|
| 121 |
|
---|
| 122 | //initSPD(density, refMat2, m2);
|
---|
| 123 | initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
|
---|
| 124 | refMat2 += refMat2.adjoint();
|
---|
| 125 | refMat2.diagonal() *= 0.5;
|
---|
| 126 |
|
---|
| 127 | refMat2.ldlt().solve(b, &refX);
|
---|
| 128 | typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
|
---|
| 129 | x = b;
|
---|
| 130 | SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
|
---|
| 131 | if (ldlt.succeeded())
|
---|
| 132 | ldlt.solveInPlace(x);
|
---|
| 133 | VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
|
---|
| 134 | }
|
---|
| 135 |
|
---|
| 136 | // test LU
|
---|
| 137 | {
|
---|
| 138 | static int count = 0;
|
---|
| 139 | SparseMatrix<Scalar> m2(rows, cols);
|
---|
| 140 | DenseMatrix refMat2(rows, cols);
|
---|
| 141 |
|
---|
| 142 | DenseVector b = DenseVector::Random(cols);
|
---|
| 143 | DenseVector refX(cols), x(cols);
|
---|
| 144 |
|
---|
| 145 | initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
|
---|
| 146 |
|
---|
| 147 | LU<DenseMatrix> refLu(refMat2);
|
---|
| 148 | refLu.solve(b, &refX);
|
---|
| 149 | #if defined(EIGEN_SUPERLU_SUPPORT) || defined(EIGEN_UMFPACK_SUPPORT)
|
---|
| 150 | Scalar refDet = refLu.determinant();
|
---|
| 151 | #endif
|
---|
| 152 | x.setZero();
|
---|
| 153 | // // SparseLU<SparseMatrix<Scalar> > (m2).solve(b,&x);
|
---|
| 154 | // // VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: default");
|
---|
| 155 | #ifdef EIGEN_SUPERLU_SUPPORT
|
---|
| 156 | {
|
---|
| 157 | x.setZero();
|
---|
| 158 | SparseLU<SparseMatrix<Scalar>,SuperLU> slu(m2);
|
---|
| 159 | if (slu.succeeded())
|
---|
| 160 | {
|
---|
| 161 | if (slu.solve(b,&x)) {
|
---|
| 162 | VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU");
|
---|
| 163 | }
|
---|
| 164 | // std::cerr << refDet << " == " << slu.determinant() << "\n";
|
---|
| 165 | if (count==0) {
|
---|
| 166 | VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
|
---|
| 167 | }
|
---|
| 168 | }
|
---|
| 169 | }
|
---|
| 170 | #endif
|
---|
| 171 | #ifdef EIGEN_UMFPACK_SUPPORT
|
---|
| 172 | {
|
---|
| 173 | // check solve
|
---|
| 174 | x.setZero();
|
---|
| 175 | SparseLU<SparseMatrix<Scalar>,UmfPack> slu(m2);
|
---|
| 176 | if (slu.succeeded()) {
|
---|
| 177 | if (slu.solve(b,&x)) {
|
---|
| 178 | if (count==0) {
|
---|
| 179 | VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: umfpack"); // FIXME solve is not very stable for complex
|
---|
| 180 | }
|
---|
| 181 | }
|
---|
| 182 | VERIFY_IS_APPROX(refDet,slu.determinant());
|
---|
| 183 | // TODO check the extracted data
|
---|
| 184 | //std::cerr << slu.matrixL() << "\n";
|
---|
| 185 | }
|
---|
| 186 | }
|
---|
| 187 | #endif
|
---|
| 188 | count++;
|
---|
| 189 | }
|
---|
| 190 |
|
---|
| 191 | }
|
---|
| 192 |
|
---|
| 193 | void test_eigen2_sparse_solvers()
|
---|
| 194 | {
|
---|
| 195 | for(int i = 0; i < g_repeat; i++) {
|
---|
| 196 | CALL_SUBTEST_1( sparse_solvers<double>(8, 8) );
|
---|
| 197 | CALL_SUBTEST_2( sparse_solvers<std::complex<double> >(16, 16) );
|
---|
| 198 | CALL_SUBTEST_1( sparse_solvers<double>(101, 101) );
|
---|
| 199 | }
|
---|
| 200 | }
|
---|