17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Desire Nuentsa Wakam <desire.nuentsa_wakam@inria.fr> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2014 Gael Guennebaud <gael.guennebaud@inria.fr> 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include "sparse.h" 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include <Eigen/SparseQR> 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType,typename DenseMat> 132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangint generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 150) 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(maxRows >= maxCols); 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int rows = internal::random<int>(1,maxRows); 182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang int cols = internal::random<int>(1,maxCols); 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double density = (std::max)(8./(rows*cols), 0.01); 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez A.resize(rows,cols); 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dA.resize(rows,cols); 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez initSparse<Scalar>(density, dA, A,ForceNonZeroDiag); 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez A.makeCompressed(); 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int nop = internal::random<int>(0, internal::random<double>(0,1) > 0.5 ? cols/2 : 0); 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int k=0; k<nop; ++k) 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int j0 = internal::random<int>(0,cols-1); 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int j1 = internal::random<int>(0,cols-1); 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar s = internal::random<Scalar>(); 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez A.col(j0) = s * A.col(j1); 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dA.col(j0) = s * dA.col(j1); 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// if(rows<cols) { 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// A.conservativeResize(cols,cols); 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// dA.conservativeResize(cols,cols); 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// dA.bottomRows(cols-rows).setZero(); 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// } 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return rows; 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Scalar> void test_sparseqr_scalar() 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef SparseMatrix<Scalar,ColMajor> MatrixType; 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar,Dynamic,Dynamic> DenseMat; 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar,Dynamic,1> DenseVector; 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType A; 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez DenseMat dA; 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez DenseVector refX,x,b; 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SparseQR<MatrixType, COLAMDOrdering<int> > solver; 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez generate_sparse_rectangular_problem(A,dA); 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b = dA * DenseVector::Random(A.cols()); 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez solver.compute(A); 572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(internal::random<float>(0,1)>0.5f) 58615d816d068b4d0f5e8df601930b5f160bf7eda1Tim Murray solver.factorize(A); // this checks that calling analyzePattern is not needed if the pattern do not change. 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (solver.info() != Success) 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cerr << "sparse QR factorization failed\n"; 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez exit(0); 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x = solver.solve(b); 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (solver.info() != Success) 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cerr << "sparse QR factorization failed\n"; 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez exit(0); 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(A * x, b); 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Compare with a dense QR solver 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColPivHouseholderQR<DenseMat> dqr(dA); 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez refX = dqr.solve(b); 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_EQUAL(dqr.rank(), solver.rank()); 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(solver.rank()==A.cols()) // full rank 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY_IS_APPROX(x, refX); 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// else 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// VERIFY((dA * refX - b).norm() * 2 > (A * x - b).norm() ); 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Compute explicitly the matrix Q 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType Q, QtQ, idM; 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Q = solver.matrixQ(); 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Check ||Q' * Q - I || 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez QtQ = Q * Q.adjoint(); 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez idM.resize(Q.rows(), Q.rows()); idM.setIdentity(); 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(idM.isApprox(QtQ)); 922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Q to dense 942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang DenseMat dQ; 952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang dQ = solver.matrixQ(); 962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang VERIFY_IS_APPROX(Q, dQ); 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid test_sparseqr() 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(int i=0; i<g_repeat; ++i) 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST_1(test_sparseqr_scalar<double>()); 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST_2(test_sparseqr_scalar<std::complex<double> >()); 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 107