sparse_cholesky.cpp File Reference

#include <iostream>
#include <Eigen/Sparse>
#include "BenchSparseUtil.h"
#include <Eigen/Cholesky>

Macros
#define	NOGMM
#define	NOMTL
#define	SIZE   10
#define	DENSITY   0.01
#define	REPEAT   1
#define	MINDENSITY   0.0004
#define	NBTRIES   10
#define	BENCH(X)

Typedefs
typedef SparseMatrix< Scalar, SelfAdjoint|LowerTriangular >	EigenSparseSelfAdjointMatrix

Functions
void	fillSpdMatrix (float density, int rows, int cols, EigenSparseSelfAdjointMatrix &dst)
template<int Backend>
void	doEigen (const char *name, const EigenSparseSelfAdjointMatrix &sm1, int flags=0)
int	main (int argc, char *argv[])

Macro Definition Documentation

BENCH

#define BENCH

(

X

)

Value:

  timer.reset();                          \
  for (int _j = 0; _j < NBTRIES; ++_j) {  \
    timer.start();                        \
    for (int _k = 0; _k < REPEAT; ++_k) { \
      X                                   \
    }                                     \
    timer.stop();                         \
  }

DENSITY

#define DENSITY   0.01

MINDENSITY

#define MINDENSITY   0.0004

NBTRIES

#define NBTRIES   10

NOGMM

#define NOGMM

NOMTL

#define NOMTL

REPEAT

#define REPEAT   1

SIZE

#define SIZE   10

Typedef Documentation

EigenSparseSelfAdjointMatrix

typedef SparseMatrix<Scalar, SelfAdjoint | LowerTriangular> EigenSparseSelfAdjointMatrix

Function Documentation

doEigen()

template<int Backend>

void doEigen	(	const char *	name,
		const EigenSparseSelfAdjointMatrix &	sm1,
		int	flags = 0
	)

                                                                                       {
  std::cout << name << "..." << std::flush;
  BenchTimer timer;
  timer.start();
  SparseLLT<EigenSparseSelfAdjointMatrix, Backend> chol(sm1, flags);
  timer.stop();
  std::cout << ":\t" << timer.value() << endl;
 
  std::cout << "  nnz: " << sm1.nonZeros() << " => " << chol.matrixL().nonZeros() << "\n";
  //   std::cout << "sparse\n" << chol.matrixL() << "%\n";
}

References plotDoE::name, and Eigen::SparseMatrix< Scalar_, Options_, StorageIndex_ >::nonZeros().

fillSpdMatrix()

void fillSpdMatrix	(	float	density,
		int	rows,
		int	cols,
		EigenSparseSelfAdjointMatrix &	dst
	)

                                                                                         {
  dst.startFill(rows * cols * density);
  for (int j = 0; j < cols; j++) {
    dst.fill(j, j) = internal::random<Scalar>(10, 20);
    for (int i = j + 1; i < rows; i++) {
      Scalar v = (internal::random<float>(0, 1) < density) ? internal::random<Scalar>() : 0;
      if (v != 0) dst.fill(i, j) = v;
    }
  }
  dst.endFill();
}

References cols, UniformPSDSelfTest::density, i, j, rows, and v.

Referenced by main().

main()

int main	(	int argc	,
		char *	argv[]
	)

                                 {
  int rows = SIZE;
  int cols = SIZE;
  float density = DENSITY;
  BenchTimer timer;
 
  VectorXf b = VectorXf::Random(cols);
  VectorXf x = VectorXf::Random(cols);
 
  bool densedone = false;
 
  // for (float density = DENSITY; density>=MINDENSITY; density*=0.5)
  //   float density = 0.5;
  {
    EigenSparseSelfAdjointMatrix sm1(rows, cols);
    std::cout << "Generate sparse matrix (might take a while)...\n";
    fillSpdMatrix(density, rows, cols, sm1);
    std::cout << "DONE\n\n";
 
// dense matrices
#ifdef DENSEMATRIX
    if (!densedone) {
      densedone = true;
      std::cout << "Eigen Dense\t" << density * 100 << "%\n";
      DenseMatrix m1(rows, cols);
      eiToDense(sm1, m1);
      m1 = (m1 + m1.transpose()).eval();
      m1.diagonal() *= 0.5;
 
      //       BENCH(LLT<DenseMatrix> chol(m1);)
      //       std::cout << "dense:\t" << timer.value() << endl;
 
      BenchTimer timer;
      timer.start();
      LLT<DenseMatrix> chol(m1);
      timer.stop();
      std::cout << "dense:\t" << timer.value() << endl;
      int count = 0;
      for (int j = 0; j < cols; ++j)
        for (int i = j; i < rows; ++i)
          if (!internal::isMuchSmallerThan(internal::abs(chol.matrixL()(i, j)), 0.1)) count++;
      std::cout << "dense: "
                << "nnz = " << count << "\n";
      //       std::cout << "dense:\n" << m1 << "\n\n" << chol.matrixL() << endl;
    }
#endif
 
    // eigen sparse matrices
    doEigen<Eigen::DefaultBackend>("Eigen/Sparse", sm1, Eigen::IncompleteFactorization);
 
#ifdef EIGEN_CHOLMOD_SUPPORT
    doEigen<Eigen::Cholmod>("Eigen/Cholmod", sm1, Eigen::IncompleteFactorization);
#endif
 
#ifdef EIGEN_TAUCS_SUPPORT
    doEigen<Eigen::Taucs>("Eigen/Taucs", sm1, Eigen::IncompleteFactorization);
#endif
 
#if 0
    // TAUCS
    {
      taucs_ccs_matrix A = sm1.asTaucsMatrix();
 
      //BENCH(taucs_ccs_matrix* chol = taucs_ccs_factor_llt(&A, 0, 0);)
//       BENCH(taucs_supernodal_factor_to_ccs(taucs_ccs_factor_llt_ll(&A));)
//       std::cout << "taucs:\t" << timer.value() << endl;
 
      taucs_ccs_matrix* chol = taucs_ccs_factor_llt(&A, 0, 0);
 
      for (int j=0; j<cols; ++j)
      {
        for (int i=chol->colptr[j]; i<chol->colptr[j+1]; ++i)
          std::cout << chol->values.d[i] << " ";
      }
    }
 
    // CHOLMOD
#ifdef EIGEN_CHOLMOD_SUPPORT
    {
      cholmod_common c;
      cholmod_start (&c);
      cholmod_sparse A;
      cholmod_factor *L;
 
      A = sm1.asCholmodMatrix();
      BenchTimer timer;
//       timer.reset();
      timer.start();
      std::vector<int> perm(cols);
//       std::vector<int> set(ncols);
      for (int i=0; i<cols; ++i)
        perm[i] = i;
//       c.nmethods = 1;
//       c.method[0] = 1;
 
      c.nmethods = 1;
      c.method [0].ordering = CHOLMOD_NATURAL;
      c.postorder = 0;
      c.final_ll = 1;
 
      L = cholmod_analyze_p(&A, &perm[0], &perm[0], cols, &c);
      timer.stop();
      std::cout << "cholmod/analyze:\t" << timer.value() << endl;
      timer.reset();
      timer.start();
      cholmod_factorize(&A, L, &c);
      timer.stop();
      std::cout << "cholmod/factorize:\t" << timer.value() << endl;
 
      cholmod_sparse* cholmat = cholmod_factor_to_sparse(L, &c);
 
      cholmod_print_factor(L, "Factors", &c);
 
      cholmod_print_sparse(cholmat, "Chol", &c);
      cholmod_write_sparse(stdout, cholmat, 0, 0, &c);
//
//       cholmod_print_sparse(&A, "A", &c);
//       cholmod_write_sparse(stdout, &A, 0, 0, &c);
 
 
//       for (int j=0; j<cols; ++j)
//       {
//           for (int i=chol->colptr[j]; i<chol->colptr[j+1]; ++i)
//             std::cout << chol->values.s[i] << " ";
//       }
    }
#endif
 
#endif
  }
 
  return 0;
}

References abs(), b, calibrate::c, cols, UniformPSDSelfTest::density, DENSITY, eiToDense(), eval(), fillSpdMatrix(), i, Eigen::internal::isMuchSmallerThan(), j, L, m1, Eigen::LLT< MatrixType_, UpLo_ >::matrixL(), rows, SIZE, fix_broken_doxygen_formulae::stdout, and plotDoE::x.