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 (at) inria.fr> 5 // Copyright (C) 2014 Gael Guennebaud <gael.guennebaud (at) 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> 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 44 template<typename Scalar> void test_sparseqr_scalar() 45 { 46 typedef SparseMatrix<Scalar,ColMajor> MatrixType; 47 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMat; 48 typedef Matrix<Scalar,Dynamic,1> DenseVector; 49 MatrixType A; 50 DenseMat dA; 51 DenseVector refX,x,b; 52 SparseQR<MatrixType, COLAMDOrdering<int> > solver; 53 generate_sparse_rectangular_problem(A,dA); 54 55 b = dA * DenseVector::Random(A.cols()); 56 solver.compute(A); 57 if(internal::random<float>(0,1)>0.5f) 58 solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change. 59 if (solver.info() != Success) 60 { 61 std::cerr << "sparse QR factorization failed\n"; 62 exit(0); 63 return; 64 } 65 x = solver.solve(b); 66 if (solver.info() != Success) 67 { 68 std::cerr << "sparse QR factorization failed\n"; 69 exit(0); 70 return; 71 } 72 73 VERIFY_IS_APPROX(A * x, b); 74 75 //Compare with a dense QR solver 76 ColPivHouseholderQR<DenseMat> dqr(dA); 77 refX = dqr.solve(b); 78 79 VERIFY_IS_EQUAL(dqr.rank(), solver.rank()); 80 if(solver.rank()==A.cols()) // full rank 81 VERIFY_IS_APPROX(x, refX); 82 // else 83 // VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() ); 84 85 // Compute explicitly the matrix Q 86 MatrixType Q, QtQ, idM; 87 Q = solver.matrixQ(); 88 //Check ||Q' * Q - I || 89 QtQ = Q * Q.adjoint(); 90 idM.resize(Q.rows(), Q.rows()); idM.setIdentity(); 91 VERIFY(idM.isApprox(QtQ)); 92 93 // Q to dense 94 DenseMat dQ; 95 dQ = solver.matrixQ(); 96 VERIFY_IS_APPROX(Q, dQ); 97 } 98 void test_sparseqr() 99 { 100 for(int i=0; i<g_repeat; ++i) 101 { 102 CALL_SUBTEST_1(test_sparseqr_scalar<double>()); 103 CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >()); 104 } 105 } 106 107