ground/gcs/src/libs/eigen/test/spqr_support.cpp

   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@inria.fr>
   5 //
   6 // This Source Code Form is subject to the terms of the Mozilla
   7 // Public License v. 2.0. If a copy of the MPL was not distributed
   8 #include "sparse.h"
   9 #include <Eigen/SPQRSupport>
  10
  11
  12 template<typename MatrixType,typename DenseMat>
  13 int generate_sparse_rectangular_problem(MatrixType& A, DenseMat& dA, int maxRows = 300, int maxCols = 300)
  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,rows);
  19   double density = (std::max)(8./(rows*cols), 0.01);
  20
  21   A.resize(rows,rows);
  22   dA.resize(rows,rows);
  23   initSparse<Scalar>(density, dA, A,ForceNonZeroDiag);
  24   A.makeCompressed();
  25   return rows;
  26 }
  27
  28 template<typename Scalar> void test_spqr_scalar()
  29 {
  30   typedef SparseMatrix<Scalar,ColMajor> MatrixType;
  31   MatrixType A;
  32   Matrix<Scalar,Dynamic,Dynamic> dA;
  33   typedef Matrix<Scalar,Dynamic,1> DenseVector;
  34   DenseVector refX,x,b;
  35   SPQR<MatrixType> solver;
  36   generate_sparse_rectangular_problem(A,dA);
  37
  38   int m = A.rows();
  39   b = DenseVector::Random(m);
  40   solver.compute(A);
  41   if (solver.info() != Success)
  42   {
  43     std::cerr << "sparse QR factorization failed\n";
  44     exit(0);
  45     return;
  46   }
  47   x = solver.solve(b);
  48   if (solver.info() != Success)
  49   {
  50     std::cerr << "sparse QR factorization failed\n";
  51     exit(0);
  52     return;
  53   }
  54   //Compare with a dense solver
  55   refX = dA.colPivHouseholderQr().solve(b);
  56   VERIFY(x.isApprox(refX,test_precision<Scalar>()));
  57 }
  58 void test_spqr_support()
  59 {
  60   CALL_SUBTEST_1(test_spqr_scalar<double>());
  61   CALL_SUBTEST_2(test_spqr_scalar<std::complex<double> >());
  62 }