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
13template<typename Solver, typename Rhs, typename DenseMat, typename DenseRhs>
14void 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
72template<typename Solver, typename Rhs>
73void 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}
111template<typename Solver, typename DenseMat>
112void 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}
127template<typename Solver, typename DenseMat>
128void 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
145template<typename Solver, typename DenseMat>
146int 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
174template<typename Scalar>
175inline 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
186template<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
246template<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
263template<typename Solver, typename DenseMat>
264int 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
280template<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
331template<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
347template<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