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@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
12template<typename Scalar> void
13initSPD(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
33template<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
104void 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