1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2011 Gael Guennebaud <g.gael@free.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 #include <Eigen/SparseCore>
12 
13 template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs>
check_sparse_solving(Solver & solver,const typename Solver::MatrixType & A,const Rhs & b,const DenseMat & dA,const DenseRhs & db)14 void check_sparse_solving(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const DenseMat& dA, const DenseRhs& db)
15 {
16   typedef typename Solver::MatrixType Mat;
17   typedef typename Mat::Scalar Scalar;
18 
19   DenseRhs refX = dA.lu().solve(db);
20   {
21     Rhs x(b.rows(), b.cols());
22     Rhs oldb = b;
23 
24     solver.compute(A);
25     if (solver.info() != Success)
26     {
27       std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n";
28       exit(0);
29       return;
30     }
31     x = solver.solve(b);
32     if (solver.info() != Success)
33     {
34       std::cerr << "sparse solver testing: solving failed\n";
35       return;
36     }
37     VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!");
38 
39     VERIFY(x.isApprox(refX,test_precision<Scalar>()));
40     x.setZero();
41     // test the analyze/factorize API
42     solver.analyzePattern(A);
43     solver.factorize(A);
44     if (solver.info() != Success)
45     {
46       std::cerr << "sparse solver testing: factorization failed (check_sparse_solving)\n";
47       exit(0);
48       return;
49     }
50     x = solver.solve(b);
51     if (solver.info() != Success)
52     {
53       std::cerr << "sparse solver testing: solving failed\n";
54       return;
55     }
56     VERIFY(oldb.isApprox(b) && "sparse solver testing: the rhs should not be modified!");
57 
58     VERIFY(x.isApprox(refX,test_precision<Scalar>()));
59   }
60 
61   // test dense Block as the result and rhs:
62   {
63     DenseRhs x(db.rows(), db.cols());
64     DenseRhs oldb(db);
65     x.setZero();
66     x.block(0,0,x.rows(),x.cols()) = solver.solve(db.block(0,0,db.rows(),db.cols()));
67     VERIFY(oldb.isApprox(db) && "sparse solver testing: the rhs should not be modified!");
68     VERIFY(x.isApprox(refX,test_precision<Scalar>()));
69   }
70 }
71 
72 template<typename Solver, typename Rhs>
check_sparse_solving_real_cases(Solver & solver,const typename Solver::MatrixType & A,const Rhs & b,const Rhs & refX)73 void check_sparse_solving_real_cases(Solver& solver, const typename Solver::MatrixType& A, const Rhs& b, const Rhs& refX)
74 {
75   typedef typename Solver::MatrixType Mat;
76   typedef typename Mat::Scalar Scalar;
77   typedef typename Mat::RealScalar RealScalar;
78 
79   Rhs x(b.rows(), b.cols());
80 
81   solver.compute(A);
82   if (solver.info() != Success)
83   {
84     std::cerr << "sparse solver testing: factorization failed (check_sparse_solving_real_cases)\n";
85     exit(0);
86     return;
87   }
88   x = solver.solve(b);
89   if (solver.info() != Success)
90   {
91     std::cerr << "sparse solver testing: solving failed\n";
92     return;
93   }
94 
95   RealScalar res_error;
96   // Compute the norm of the relative error
97   if(refX.size() != 0)
98     res_error = (refX - x).norm()/refX.norm();
99   else
100   {
101     // Compute the relative residual norm
102     res_error = (b - A * x).norm()/b.norm();
103   }
104   if (res_error > test_precision<Scalar>() ){
105     std::cerr << "Test " << g_test_stack.back() << " failed in "EI_PP_MAKE_STRING(__FILE__)
106     << " (" << EI_PP_MAKE_STRING(__LINE__) << ")" << std::endl << std::endl;
107     abort();
108   }
109 
110 }
111 template<typename Solver, typename DenseMat>
check_sparse_determinant(Solver & solver,const typename Solver::MatrixType & A,const DenseMat & dA)112 void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
113 {
114   typedef typename Solver::MatrixType Mat;
115   typedef typename Mat::Scalar Scalar;
116 
117   solver.compute(A);
118   if (solver.info() != Success)
119   {
120     std::cerr << "sparse solver testing: factorization failed (check_sparse_determinant)\n";
121     return;
122   }
123 
124   Scalar refDet = dA.determinant();
125   VERIFY_IS_APPROX(refDet,solver.determinant());
126 }
127 template<typename Solver, typename DenseMat>
check_sparse_abs_determinant(Solver & solver,const typename Solver::MatrixType & A,const DenseMat & dA)128 void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
129 {
130   using std::abs;
131   typedef typename Solver::MatrixType Mat;
132   typedef typename Mat::Scalar Scalar;
133 
134   solver.compute(A);
135   if (solver.info() != Success)
136   {
137     std::cerr << "sparse solver testing: factorization failed (check_sparse_abs_determinant)\n";
138     return;
139   }
140 
141   Scalar refDet = abs(dA.determinant());
142   VERIFY_IS_APPROX(refDet,solver.absDeterminant());
143 }
144 
145 template<typename Solver, typename DenseMat>
146 int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300)
147 {
148   typedef typename Solver::MatrixType Mat;
149   typedef typename Mat::Scalar Scalar;
150   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
151 
152   int size = internal::random<int>(1,maxSize);
153   double density = (std::max)(8./(size*size), 0.01);
154 
155   Mat M(size, size);
156   DenseMatrix dM(size, size);
157 
158   initSparse<Scalar>(density, dM, M, ForceNonZeroDiag);
159 
160   A = M * M.adjoint();
161   dA = dM * dM.adjoint();
162 
163   halfA.resize(size,size);
164   if(Solver::UpLo==(Lower|Upper))
165     halfA = A;
166   else
167     halfA.template selfadjointView<Solver::UpLo>().rankUpdate(M);
168 
169   return size;
170 }
171 
172 
173 #ifdef TEST_REAL_CASES
174 template<typename Scalar>
get_matrixfolder()175 inline std::string get_matrixfolder()
176 {
177   std::string mat_folder = TEST_REAL_CASES;
178   if( internal::is_same<Scalar, std::complex<float> >::value || internal::is_same<Scalar, std::complex<double> >::value )
179     mat_folder  = mat_folder + static_cast<std::string>("/complex/");
180   else
181     mat_folder = mat_folder + static_cast<std::string>("/real/");
182   return mat_folder;
183 }
184 #endif
185 
check_sparse_spd_solving(Solver & solver)186 template<typename Solver> void check_sparse_spd_solving(Solver& solver)
187 {
188   typedef typename Solver::MatrixType Mat;
189   typedef typename Mat::Scalar Scalar;
190   typedef SparseMatrix<Scalar,ColMajor> SpMat;
191   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
192   typedef Matrix<Scalar,Dynamic,1> DenseVector;
193 
194   // generate the problem
195   Mat A, halfA;
196   DenseMatrix dA;
197   for (int i = 0; i < g_repeat; i++) {
198     int size = generate_sparse_spd_problem(solver, A, halfA, dA);
199 
200     // generate the right hand sides
201     int rhsCols = internal::random<int>(1,16);
202     double density = (std::max)(8./(size*rhsCols), 0.1);
203     SpMat B(size,rhsCols);
204     DenseVector b = DenseVector::Random(size);
205     DenseMatrix dB(size,rhsCols);
206     initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
207 
208     check_sparse_solving(solver, A,     b,  dA, b);
209     check_sparse_solving(solver, halfA, b,  dA, b);
210     check_sparse_solving(solver, A,     dB, dA, dB);
211     check_sparse_solving(solver, halfA, dB, dA, dB);
212     check_sparse_solving(solver, A,     B,  dA, dB);
213     check_sparse_solving(solver, halfA, B,  dA, dB);
214 
215     // check only once
216     if(i==0)
217     {
218       b = DenseVector::Zero(size);
219       check_sparse_solving(solver, A, b, dA, b);
220     }
221   }
222 
223   // First, get the folder
224 #ifdef TEST_REAL_CASES
225   if (internal::is_same<Scalar, float>::value
226       || internal::is_same<Scalar, std::complex<float> >::value)
227     return ;
228 
229   std::string mat_folder = get_matrixfolder<Scalar>();
230   MatrixMarketIterator<Scalar> it(mat_folder);
231   for (; it; ++it)
232   {
233     if (it.sym() == SPD){
234       Mat halfA;
235       PermutationMatrix<Dynamic, Dynamic, Index> pnull;
236       halfA.template selfadjointView<Solver::UpLo>() = it.matrix().template triangularView<Eigen::Lower>().twistedBy(pnull);
237 
238       std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n";
239       check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX());
240       check_sparse_solving_real_cases(solver, halfA, it.rhs(), it.refX());
241     }
242   }
243 #endif
244 }
245 
check_sparse_spd_determinant(Solver & solver)246 template<typename Solver> void check_sparse_spd_determinant(Solver& solver)
247 {
248   typedef typename Solver::MatrixType Mat;
249   typedef typename Mat::Scalar Scalar;
250   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
251 
252   // generate the problem
253   Mat A, halfA;
254   DenseMatrix dA;
255   generate_sparse_spd_problem(solver, A, halfA, dA, 30);
256 
257   for (int i = 0; i < g_repeat; i++) {
258     check_sparse_determinant(solver, A,     dA);
259     check_sparse_determinant(solver, halfA, dA );
260   }
261 }
262 
263 template<typename Solver, typename DenseMat>
264 int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, DenseMat& dA, int maxSize = 300)
265 {
266   typedef typename Solver::MatrixType Mat;
267   typedef typename Mat::Scalar Scalar;
268 
269   int size = internal::random<int>(1,maxSize);
270   double density = (std::max)(8./(size*size), 0.01);
271 
272   A.resize(size,size);
273   dA.resize(size,size);
274 
275   initSparse<Scalar>(density, dA, A, ForceNonZeroDiag);
276 
277   return size;
278 }
279 
check_sparse_square_solving(Solver & solver)280 template<typename Solver> void check_sparse_square_solving(Solver& solver)
281 {
282   typedef typename Solver::MatrixType Mat;
283   typedef typename Mat::Scalar Scalar;
284   typedef SparseMatrix<Scalar,ColMajor> SpMat;
285   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
286   typedef Matrix<Scalar,Dynamic,1> DenseVector;
287 
288   int rhsCols = internal::random<int>(1,16);
289 
290   Mat A;
291   DenseMatrix dA;
292   for (int i = 0; i < g_repeat; i++) {
293     int size = generate_sparse_square_problem(solver, A, dA);
294 
295     A.makeCompressed();
296     DenseVector b = DenseVector::Random(size);
297     DenseMatrix dB(size,rhsCols);
298     SpMat B(size,rhsCols);
299     double density = (std::max)(8./(size*rhsCols), 0.1);
300     initSparse<Scalar>(density, dB, B, ForceNonZeroDiag);
301     B.makeCompressed();
302     check_sparse_solving(solver, A, b,  dA, b);
303     check_sparse_solving(solver, A, dB, dA, dB);
304     check_sparse_solving(solver, A, B,  dA, dB);
305 
306     // check only once
307     if(i==0)
308     {
309       b = DenseVector::Zero(size);
310       check_sparse_solving(solver, A, b, dA, b);
311     }
312   }
313 
314   // First, get the folder
315 #ifdef TEST_REAL_CASES
316   if (internal::is_same<Scalar, float>::value
317       || internal::is_same<Scalar, std::complex<float> >::value)
318     return ;
319 
320   std::string mat_folder = get_matrixfolder<Scalar>();
321   MatrixMarketIterator<Scalar> it(mat_folder);
322   for (; it; ++it)
323   {
324     std::cout<< " ==== SOLVING WITH MATRIX " << it.matname() << " ==== \n";
325     check_sparse_solving_real_cases(solver, it.matrix(), it.rhs(), it.refX());
326   }
327 #endif
328 
329 }
330 
check_sparse_square_determinant(Solver & solver)331 template<typename Solver> void check_sparse_square_determinant(Solver& solver)
332 {
333   typedef typename Solver::MatrixType Mat;
334   typedef typename Mat::Scalar Scalar;
335   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
336 
337   // generate the problem
338   Mat A;
339   DenseMatrix dA;
340   generate_sparse_square_problem(solver, A, dA, 30);
341   A.makeCompressed();
342   for (int i = 0; i < g_repeat; i++) {
343     check_sparse_determinant(solver, A, dA);
344   }
345 }
346 
check_sparse_square_abs_determinant(Solver & solver)347 template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver)
348 {
349   typedef typename Solver::MatrixType Mat;
350   typedef typename Mat::Scalar Scalar;
351   typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
352 
353   // generate the problem
354   Mat A;
355   DenseMatrix dA;
356   generate_sparse_square_problem(solver, A, dA, 30);
357   A.makeCompressed();
358   for (int i = 0; i < g_repeat; i++) {
359     check_sparse_abs_determinant(solver, A, dA);
360   }
361 }
362 
363