blob: f0e721fcea2cc6949ab8b4ac66a0f3434b6c7dc5 [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2012 Desire Nuentsa Wakam // Copyright (C) 2014 Gael Guennebaud // // This Source Code Form is subject to the terms of the Mozilla // Public License v. 2.0. If a copy of the MPL was not distributed #include "sparse.h" #include template int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 150) { eigen_assert(maxRows >= maxCols); typedef typename MatrixType::Scalar Scalar; int rows = internal::random(1,maxRows); int cols = internal::random(1,maxCols); double density = (std::max)(8./(rows*cols), 0.01); A.resize(rows,cols); dA.resize(rows,cols); initSparse(density, dA, A,ForceNonZeroDiag); A.makeCompressed(); int nop = internal::random(0, internal::random(0,1) > 0.5 ? cols/2 : 0); for(int k=0; k(0,cols-1); int j1 = internal::random(0,cols-1); Scalar s = internal::random(); A.col(j0) = s * A.col(j1); dA.col(j0) = s * dA.col(j1); } // if(rows void test_sparseqr_scalar() { typedef SparseMatrix MatrixType; typedef Matrix DenseMat; typedef Matrix DenseVector; MatrixType A; DenseMat dA; DenseVector refX,x,b; SparseQR > solver; generate_sparse_rectangular_problem(A,dA); b = dA * DenseVector::Random(A.cols()); solver.compute(A); // Q should be MxM VERIFY_IS_EQUAL(solver.matrixQ().rows(), A.rows()); VERIFY_IS_EQUAL(solver.matrixQ().cols(), A.rows()); // R should be MxN VERIFY_IS_EQUAL(solver.matrixR().rows(), A.rows()); VERIFY_IS_EQUAL(solver.matrixR().cols(), A.cols()); // Q and R can be multiplied DenseMat recoveredA = solver.matrixQ() * DenseMat(solver.matrixR().template triangularView()) * solver.colsPermutation().transpose(); VERIFY_IS_EQUAL(recoveredA.rows(), A.rows()); VERIFY_IS_EQUAL(recoveredA.cols(), A.cols()); // and in the full rank case the original matrix is recovered if (solver.rank() == A.cols()) { VERIFY_IS_APPROX(A, recoveredA); } if(internal::random(0,1)>0.5f) solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change. if (solver.info() != Success) { std::cerr << "sparse QR factorization failed\n"; exit(0); return; } x = solver.solve(b); if (solver.info() != Success) { std::cerr << "sparse QR factorization failed\n"; exit(0); return; } VERIFY_IS_APPROX(A * x, b); //Compare with a dense QR solver ColPivHouseholderQR dqr(dA); refX = dqr.solve(b); VERIFY_IS_EQUAL(dqr.rank(), solver.rank()); if(solver.rank()==A.cols()) // full rank VERIFY_IS_APPROX(x, refX); // else // VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() ); // Compute explicitly the matrix Q MatrixType Q, QtQ, idM; Q = solver.matrixQ(); //Check ||Q' * Q - I || QtQ = Q * Q.adjoint(); idM.resize(Q.rows(), Q.rows()); idM.setIdentity(); VERIFY(idM.isApprox(QtQ)); // Q to dense DenseMat dQ; dQ = solver.matrixQ(); VERIFY_IS_APPROX(Q, dQ); } void test_sparseqr() { for(int i=0; i()); CALL_SUBTEST_2(test_sparseqr_scalar >()); } }