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