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 | }
|
---|