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// 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include "sparse.h" 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include <Eigen/SPQRSupport> 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType,typename DenseMat> 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezint generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 300) 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(maxRows >= maxCols); 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int rows = internal::random<int>(1,maxRows); 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int cols = internal::random<int>(1,rows); 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez double density = (std::max)(8./(rows*cols), 0.01); 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez A.resize(rows,rows); 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dA.resize(rows,rows); 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez initSparse<Scalar>(density, dA, A,ForceNonZeroDiag); 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez A.makeCompressed(); 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return rows; 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename Scalar> void test_spqr_scalar() 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef SparseMatrix<Scalar,ColMajor> MatrixType; 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType A; 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar,Dynamic,Dynamic> dA; 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar,Dynamic,1> DenseVector; 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez DenseVector refX,x,b; 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez SPQR<MatrixType> solver; 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez generate_sparse_rectangular_problem(A,dA); 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez int m = A.rows(); 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b = DenseVector::Random(m); 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez solver.compute(A); 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (solver.info() != Success) 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cerr << "sparse QR factorization failed\n"; 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez exit(0); 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x = solver.solve(b); 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (solver.info() != Success) 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::cerr << "sparse QR factorization failed\n"; 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez exit(0); 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez //Compare with a dense solver 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez refX = dA.colPivHouseholderQr().solve(b); 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VERIFY(x.isApprox(refX,test_precision<Scalar>())); 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezvoid test_spqr_support() 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST_1(test_spqr_scalar<double>()); 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez CALL_SUBTEST_2(test_spqr_scalar<std::complex<double> >()); 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 63