1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr>
5 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 #include "sparse.h"
10 #include <Eigen/SparseQR>
11
12 template<typename MatrixType,typename DenseMat>
generate_sparse_rectangular_problem(MatrixType & A,DenseMat & dA,int maxRows=300,int maxCols=150)13 int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 150)
14 {
15 eigen_assert(maxRows >= maxCols);
16 typedef typename MatrixType::Scalar Scalar;
17 int rows = internal::random<int>(1,maxRows);
18 int cols = internal::random<int>(1,maxCols);
19 double density = (std::max)(8./(rows*cols), 0.01);
20
21 A.resize(rows,cols);
22 dA.resize(rows,cols);
23 initSparse<Scalar>(density, dA, A,ForceNonZeroDiag);
24 A.makeCompressed();
25 int nop = internal::random<int>(0, internal::random<double>(0,1) > 0.5 ? cols/2 : 0);
26 for(int k=0; k<nop; ++k)
27 {
28 int j0 = internal::random<int>(0,cols-1);
29 int j1 = internal::random<int>(0,cols-1);
30 Scalar s = internal::random<Scalar>();
31 A.col(j0) = s * A.col(j1);
32 dA.col(j0) = s * dA.col(j1);
33 }
34
35 // if(rows<cols) {
36 // A.conservativeResize(cols,cols);
37 // dA.conservativeResize(cols,cols);
38 // dA.bottomRows(cols-rows).setZero();
39 // }
40
41 return rows;
42 }
43
test_sparseqr_scalar()44 template<typename Scalar> void test_sparseqr_scalar()
45 {
46 typedef typename NumTraits<Scalar>::Real RealScalar;
47 typedef SparseMatrix<Scalar,ColMajor> MatrixType;
48 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMat;
49 typedef Matrix<Scalar,Dynamic,1> DenseVector;
50 MatrixType A;
51 DenseMat dA;
52 DenseVector refX,x,b;
53 SparseQR<MatrixType, COLAMDOrdering<int> > solver;
54 generate_sparse_rectangular_problem(A,dA);
55
56 b = dA * DenseVector::Random(A.cols());
57 solver.compute(A);
58
59 // Q should be MxM
60 VERIFY_IS_EQUAL(solver.matrixQ().rows(), A.rows());
61 VERIFY_IS_EQUAL(solver.matrixQ().cols(), A.rows());
62
63 // R should be MxN
64 VERIFY_IS_EQUAL(solver.matrixR().rows(), A.rows());
65 VERIFY_IS_EQUAL(solver.matrixR().cols(), A.cols());
66
67 // Q and R can be multiplied
68 DenseMat recoveredA = solver.matrixQ()
69 * DenseMat(solver.matrixR().template triangularView<Upper>())
70 * solver.colsPermutation().transpose();
71 VERIFY_IS_EQUAL(recoveredA.rows(), A.rows());
72 VERIFY_IS_EQUAL(recoveredA.cols(), A.cols());
73
74 // and in the full rank case the original matrix is recovered
75 if (solver.rank() == A.cols())
76 {
77 VERIFY_IS_APPROX(A, recoveredA);
78 }
79
80 if(internal::random<float>(0,1)>0.5f)
81 solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change.
82 if (solver.info() != Success)
83 {
84 std::cerr << "sparse QR factorization failed\n";
85 exit(0);
86 return;
87 }
88 x = solver.solve(b);
89 if (solver.info() != Success)
90 {
91 std::cerr << "sparse QR factorization failed\n";
92 exit(0);
93 return;
94 }
95
96 // Compare with a dense QR solver
97 ColPivHouseholderQR<DenseMat> dqr(dA);
98 refX = dqr.solve(b);
99
100 bool rank_deficient = A.cols()>A.rows() || dqr.rank()<A.cols();
101 if(rank_deficient)
102 {
103 // rank deficient problem -> we might have to increase the threshold
104 // to get a correct solution.
105 RealScalar th = RealScalar(20)*dA.colwise().norm().maxCoeff()*(A.rows()+A.cols()) * NumTraits<RealScalar>::epsilon();
106 for(Index k=0; (k<16) && !test_isApprox(A*x,b); ++k)
107 {
108 th *= RealScalar(10);
109 solver.setPivotThreshold(th);
110 solver.compute(A);
111 x = solver.solve(b);
112 }
113 }
114
115 VERIFY_IS_APPROX(A * x, b);
116
117 // For rank deficient problem, the estimated rank might
118 // be slightly off, so let's only raise a warning in such cases.
119 if(rank_deficient) ++g_test_level;
120 VERIFY_IS_EQUAL(solver.rank(), dqr.rank());
121 if(rank_deficient) --g_test_level;
122
123 if(solver.rank()==A.cols()) // full rank
124 VERIFY_IS_APPROX(x, refX);
125 // else
126 // VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() );
127
128 // Compute explicitly the matrix Q
129 MatrixType Q, QtQ, idM;
130 Q = solver.matrixQ();
131 //Check ||Q' * Q - I ||
132 QtQ = Q * Q.adjoint();
133 idM.resize(Q.rows(), Q.rows()); idM.setIdentity();
134 VERIFY(idM.isApprox(QtQ));
135
136 // Q to dense
137 DenseMat dQ;
138 dQ = solver.matrixQ();
139 VERIFY_IS_APPROX(Q, dQ);
140 }
EIGEN_DECLARE_TEST(sparseqr)141 EIGEN_DECLARE_TEST(sparseqr)
142 {
143 for(int i=0; i<g_repeat; ++i)
144 {
145 CALL_SUBTEST_1(test_sparseqr_scalar<double>());
146 CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >());
147 }
148 }
149
150