eigen/test/sparse_solvers.cpp

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include "sparse.h"

template <typename Scalar>
void initSPD(double density, Matrix<Scalar, Dynamic, Dynamic>& refMat, SparseMatrix<Scalar>& sparseMat) {
  Matrix<Scalar, Dynamic, Dynamic> aux(refMat.rows(), refMat.cols());
  initSparse(density, refMat, sparseMat);
  refMat = refMat * refMat.adjoint();
  for (int k = 0; k < 2; ++k) {
    initSparse(density, aux, sparseMat, ForceNonZeroDiag);
    refMat += aux * aux.adjoint();
  }
  sparseMat.setZero();
  for (int j = 0; j < sparseMat.cols(); ++j)
    for (int i = j; i < sparseMat.rows(); ++i)
      if (refMat(i, j) != Scalar(0)) sparseMat.insert(i, j) = refMat(i, j);
  sparseMat.finalize();
}

template <typename Scalar>
void sparse_solvers(int rows, int cols) {
  double density = (std::max)(8. / (rows * cols), 0.01);
  typedef Matrix<Scalar, Dynamic, Dynamic> DenseMatrix;
  typedef Matrix<Scalar, Dynamic, 1> DenseVector;
  // Scalar eps = 1e-6;

  DenseVector vec1 = DenseVector::Random(rows);

  std::vector<Vector2i> zeroCoords;
  std::vector<Vector2i> nonzeroCoords;

  // test triangular solver
  {
    DenseVector vec2 = vec1, vec3 = vec1;
    SparseMatrix<Scalar> m2(rows, cols);
    DenseMatrix refMat2 = DenseMatrix::Zero(rows, cols);

    // lower - dense
    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
    VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
                     m2.template triangularView<Lower>().solve(vec3));

    // upper - dense
    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
    VERIFY_IS_APPROX(refMat2.template triangularView<Upper>().solve(vec2),
                     m2.template triangularView<Upper>().solve(vec3));
    VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
                     m2.conjugate().template triangularView<Upper>().solve(vec3));
    {
      SparseMatrix<Scalar> cm2(m2);
      // Index rows, Index cols, Index nnz, Index* outerIndexPtr, Index* innerIndexPtr, Scalar* valuePtr
      Map<SparseMatrix<Scalar> > mm2(rows, cols, cm2.nonZeros(), cm2.outerIndexPtr(), cm2.innerIndexPtr(),
                                     cm2.valuePtr());
      VERIFY_IS_APPROX(refMat2.conjugate().template triangularView<Upper>().solve(vec2),
                       mm2.conjugate().template triangularView<Upper>().solve(vec3));
    }

    // lower - transpose
    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
    VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Upper>().solve(vec2),
                     m2.transpose().template triangularView<Upper>().solve(vec3));

    // upper - transpose
    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular, &zeroCoords, &nonzeroCoords);
    VERIFY_IS_APPROX(refMat2.transpose().template triangularView<Lower>().solve(vec2),
                     m2.transpose().template triangularView<Lower>().solve(vec3));

    SparseMatrix<Scalar> matB(rows, rows);
    DenseMatrix refMatB = DenseMatrix::Zero(rows, rows);

    // lower - sparse
    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular);
    initSparse<Scalar>(density, refMatB, matB);
    refMat2.template triangularView<Lower>().solveInPlace(refMatB);
    m2.template triangularView<Lower>().solveInPlace(matB);
    VERIFY_IS_APPROX(matB.toDense(), refMatB);

    // upper - sparse
    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeUpperTriangular);
    initSparse<Scalar>(density, refMatB, matB);
    refMat2.template triangularView<Upper>().solveInPlace(refMatB);
    m2.template triangularView<Upper>().solveInPlace(matB);
    VERIFY_IS_APPROX(matB, refMatB);

    // test deprecated API
    initSparse<Scalar>(density, refMat2, m2, ForceNonZeroDiag | MakeLowerTriangular, &zeroCoords, &nonzeroCoords);
    VERIFY_IS_APPROX(refMat2.template triangularView<Lower>().solve(vec2),
                     m2.template triangularView<Lower>().solve(vec3));

    // test empty triangular matrix
    {
      m2.resize(0, 0);
      refMatB.resize(0, refMatB.cols());
      DenseMatrix res = m2.template triangularView<Lower>().solve(refMatB);
      VERIFY_IS_EQUAL(res.rows(), 0);
      VERIFY_IS_EQUAL(res.cols(), refMatB.cols());
      res = refMatB;
      m2.template triangularView<Lower>().solveInPlace(res);
      VERIFY_IS_EQUAL(res.rows(), 0);
      VERIFY_IS_EQUAL(res.cols(), refMatB.cols());
    }
  }
}

EIGEN_DECLARE_TEST(sparse_solvers) {
  for (int i = 0; i < g_repeat; i++) {
    CALL_SUBTEST_1(sparse_solvers<double>(8, 8));
    int s = internal::random<int>(1, 300);
    CALL_SUBTEST_2(sparse_solvers<std::complex<double> >(s, s));
    CALL_SUBTEST_1(sparse_solvers<double>(s, s));
  }
}