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) 14 { 15 typedef typename MatrixType::Scalar Scalar; 16 int rows = internal::random<int>(1,maxRows); 17 int cols = internal::random<int>(1,rows); 18 double density = (std::max)(8./(rows*cols), 0.01); 19 20 A.resize(rows,cols); 21 dA.resize(rows,cols); 22 initSparse<Scalar>(density, dA, A,ForceNonZeroDiag); 23 A.makeCompressed(); 24 int nop = internal::random<int>(0, internal::random<double>(0,1) > 0.5 ? cols/2 : 0); 25 for(int k=0; k<nop; ++k) 26 { 27 int j0 = internal::random<int>(0,cols-1); 28 int j1 = internal::random<int>(0,cols-1); 29 Scalar s = internal::random<Scalar>(); 30 A.col(j0) = s * A.col(j1); 31 dA.col(j0) = s * dA.col(j1); 32 } 33 34 // if(rows<cols) { 35 // A.conservativeResize(cols,cols); 36 // dA.conservativeResize(cols,cols); 37 // dA.bottomRows(cols-rows).setZero(); 38 // } 39 40 return rows; 41 } 42 43 template<typename Scalar> void test_sparseqr_scalar() 44 { 45 typedef SparseMatrix<Scalar,ColMajor> MatrixType; 46 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMat; 47 typedef Matrix<Scalar,Dynamic,1> DenseVector; 48 MatrixType A; 49 DenseMat dA; 50 DenseVector refX,x,b; 51 SparseQR<MatrixType, COLAMDOrdering<int> > solver; 52 generate_sparse_rectangular_problem(A,dA); 53 54 b = dA * DenseVector::Random(A.cols()); 55 solver.compute(A); 56 if(internal::random<float>(0,1)>0.5) 57 solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change. 58 if (solver.info() != Success) 59 { 60 std::cerr << "sparse QR factorization failed\n"; 61 exit(0); 62 return; 63 } 64 x = solver.solve(b); 65 if (solver.info() != Success) 66 { 67 std::cerr << "sparse QR factorization failed\n"; 68 exit(0); 69 return; 70 } 71 72 VERIFY_IS_APPROX(A * x, b); 73 74 //Compare with a dense QR solver 75 ColPivHouseholderQR<DenseMat> dqr(dA); 76 refX = dqr.solve(b); 77 78 VERIFY_IS_EQUAL(dqr.rank(), solver.rank()); 79 if(solver.rank()==A.cols()) // full rank 80 VERIFY_IS_APPROX(x, refX); 81 // else 82 // VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() ); 83 84 // Compute explicitly the matrix Q 85 MatrixType Q, QtQ, idM; 86 Q = solver.matrixQ(); 87 //Check ||Q' * Q - I || 88 QtQ = Q * Q.adjoint(); 89 idM.resize(Q.rows(), Q.rows()); idM.setIdentity(); 90 VERIFY(idM.isApprox(QtQ)); 91 } 92 void test_sparseqr() 93 { 94 for(int i=0; i<g_repeat; ++i) 95 { 96 CALL_SUBTEST_1(test_sparseqr_scalar<double>()); 97 CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >()); 98 } 99 } 100 101