Home | History | Annotate | Download | only in test 
      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