eigen2_sparse_solvers.cpp
Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra. Eigen itself is part of the KDE project.
00003 //
00004 // Copyright (C) 2008 Daniel Gomez Ferro <dgomezferro@gmail.com>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #include "sparse.h"
00026 
00027 template<typename Scalar> void
00028 initSPD(double density,
00029         Matrix<Scalar,Dynamic,Dynamic>& refMat,
00030         SparseMatrix<Scalar>& sparseMat)
00031 {
00032   Matrix<Scalar,Dynamic,Dynamic> aux(refMat.rows(),refMat.cols());
00033   initSparse(density,refMat,sparseMat);
00034   refMat = refMat * refMat.adjoint();
00035   for (int k=0; k<2; ++k)
00036   {
00037     initSparse(density,aux,sparseMat,ForceNonZeroDiag);
00038     refMat += aux * aux.adjoint();
00039   }
00040   sparseMat.startFill();
00041   for (int j=0 ; j<sparseMat.cols(); ++j)
00042     for (int i=j ; i<sparseMat.rows(); ++i)
00043       if (refMat(i,j)!=Scalar(0))
00044         sparseMat.fill(i,j) = refMat(i,j);
00045   sparseMat.endFill();
00046 }
00047 
00048 template<typename Scalar> void sparse_solvers(int rows, int cols)
00049 {
00050   double density = std::max(8./(rows*cols), 0.01);
00051   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
00052   typedef Matrix<Scalar,Dynamic,1> DenseVector;
00053   // Scalar eps = 1e-6;
00054 
00055   DenseVector vec1 = DenseVector::Random(rows);
00056 
00057   std::vector<Vector2i> zeroCoords;
00058   std::vector<Vector2i> nonzeroCoords;
00059 
00060   // test triangular solver
00061   {
00062     DenseVector vec2 = vec1, vec3 = vec1;
00063     SparseMatrix<Scalar> m2(rows, cols);
00064     DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);
00065 
00066     // lower
00067     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
00068     VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().solveTriangular(vec2),
00069                      m2.template marked<LowerTriangular>().solveTriangular(vec3));
00070 
00071     // lower - transpose
00072     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
00073     VERIFY_IS_APPROX(refMat2.template marked<LowerTriangular>().transpose().solveTriangular(vec2),
00074                      m2.template marked<LowerTriangular>().transpose().solveTriangular(vec3));
00075 
00076     // upper
00077     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
00078     VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().solveTriangular(vec2),
00079                      m2.template marked<UpperTriangular>().solveTriangular(vec3));
00080 
00081     // upper - transpose
00082     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
00083     VERIFY_IS_APPROX(refMat2.template marked<UpperTriangular>().transpose().solveTriangular(vec2),
00084                      m2.template marked<UpperTriangular>().transpose().solveTriangular(vec3));
00085   }
00086 
00087   // test LLT
00088   {
00089     // TODO fix the issue with complex (see SparseLLT::solveInPlace)
00090     SparseMatrix<Scalar> m2(rows, cols);
00091     DenseMatrix refMat2(rows, cols);
00092 
00093     DenseVector b = DenseVector::Random(cols);
00094     DenseVector refX(cols), x(cols);
00095 
00096     initSPD(density, refMat2, m2);
00097 
00098     refMat2.llt().solve(b, &refX);
00099     typedef SparseMatrix<Scalar,LowerTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
00100     if (!NumTraits<Scalar>::IsComplex)
00101     {
00102       x = b;
00103       SparseLLT<SparseSelfAdjointMatrix> (m2).solveInPlace(x);
00104       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: default");
00105     }
00106     #ifdef EIGEN_CHOLMOD_SUPPORT
00107     x = b;
00108     SparseLLT<SparseSelfAdjointMatrix,Cholmod>(m2).solveInPlace(x);
00109     VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: cholmod");
00110     #endif
00111     if (!NumTraits<Scalar>::IsComplex)
00112     {
00113       #ifdef EIGEN_TAUCS_SUPPORT
00114       x = b;
00115       SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,IncompleteFactorization).solveInPlace(x);
00116       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (IncompleteFactorization)");
00117       x = b;
00118       SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalMultifrontal).solveInPlace(x);
00119       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalMultifrontal)");
00120       x = b;
00121       SparseLLT<SparseSelfAdjointMatrix,Taucs>(m2,SupernodalLeftLooking).solveInPlace(x);
00122       VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LLT: taucs (SupernodalLeftLooking)");
00123       #endif
00124     }
00125   }
00126 
00127   // test LDLT
00128   if (!NumTraits<Scalar>::IsComplex)
00129   {
00130     // TODO fix the issue with complex (see SparseLDLT::solveInPlace)
00131     SparseMatrix<Scalar> m2(rows, cols);
00132     DenseMatrix refMat2(rows, cols);
00133 
00134     DenseVector b = DenseVector::Random(cols);
00135     DenseVector refX(cols), x(cols);
00136 
00137     //initSPD(density, refMat2, m2);
00138     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag|MakeUpperTriangular, 0, 0);
00139     refMat2 += refMat2.adjoint();
00140     refMat2.diagonal() *= 0.5;
00141 
00142     refMat2.ldlt().solve(b, &refX);
00143     typedef SparseMatrix<Scalar,UpperTriangular|SelfAdjoint> SparseSelfAdjointMatrix;
00144     x = b;
00145     SparseLDLT<SparseSelfAdjointMatrix> ldlt(m2);
00146     if (ldlt.succeeded())
00147       ldlt.solveInPlace(x);
00148     VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LDLT: default");
00149   }
00150 
00151   // test LU
00152   {
00153     static int count = 0;
00154     SparseMatrix<Scalar> m2(rows, cols);
00155     DenseMatrix refMat2(rows, cols);
00156 
00157     DenseVector b = DenseVector::Random(cols);
00158     DenseVector refX(cols), x(cols);
00159 
00160     initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag, &zeroCoords, &nonzeroCoords);
00161 
00162     LU<DenseMatrix> refLu(refMat2);
00163     refLu.solve(b, &refX);
00164     #if defined(EIGEN_SUPERLU_SUPPORT) || defined(EIGEN_UMFPACK_SUPPORT)
00165     Scalar refDet = refLu.determinant();
00166     #endif
00167     x.setZero();
00168     // // SparseLU<SparseMatrix<Scalar> > (m2).solve(b,&x);
00169     // // VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: default");
00170     #ifdef EIGEN_SUPERLU_SUPPORT
00171     {
00172       x.setZero();
00173       SparseLU<SparseMatrix<Scalar>,SuperLU> slu(m2);
00174       if (slu.succeeded())
00175       {
00176         if (slu.solve(b,&x)) {
00177           VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: SuperLU");
00178         }
00179         // std::cerr << refDet << " == " << slu.determinant() << "\n";
00180         if (count==0) {
00181           VERIFY_IS_APPROX(refDet,slu.determinant()); // FIXME det is not very stable for complex
00182         }
00183       }
00184     }
00185     #endif
00186     #ifdef EIGEN_UMFPACK_SUPPORT
00187     {
00188       // check solve
00189       x.setZero();
00190       SparseLU<SparseMatrix<Scalar>,UmfPack> slu(m2);
00191       if (slu.succeeded()) {
00192         if (slu.solve(b,&x)) {
00193           if (count==0) {
00194             VERIFY(refX.isApprox(x,test_precision<Scalar>()) && "LU: umfpack");  // FIXME solve is not very stable for complex
00195           }
00196         }
00197         VERIFY_IS_APPROX(refDet,slu.determinant());
00198         // TODO check the extracted data
00199         //std::cerr << slu.matrixL() << "\n";
00200       }
00201     }
00202     #endif
00203     count++;
00204   }
00205 
00206 }
00207 
00208 void test_eigen2_sparse_solvers()
00209 {
00210   for(int i = 0; i < g_repeat; i++) {
00211     CALL_SUBTEST_1( sparse_solvers<double>(8, 8) );
00212     CALL_SUBTEST_2( sparse_solvers<std::complex<double> >(16, 16) );
00213     CALL_SUBTEST_1( sparse_solvers<double>(101, 101) );
00214   }
00215 }


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:39