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
12 template<typename Scalar> void
initSPD(double density,Matrix<Scalar,Dynamic,Dynamic> & refMat,SparseMatrix<Scalar> & sparseMat)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
sparse_solvers(int rows,int cols)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
test_sparse_solvers()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