Hi Gael, Hi Benoît,

Thank you for Eigen, it really looks very promising.

I'd find very useful to have std::complex support for Eigen::Cholmod backend. Do you plan to add it soon?

I'm new to c++ templating, so I'll try to catch up with this and Eigen and provide help later if possible.

Cheers
Romain

I'd find very useful to have std::complex support for Eigen::Cholmod backend. Do you plan to add it soon?

That's a question for Gael....

I'm new to c++ templating, so I'll try to catch up with this and Eigen and provide help later if possible.

Sounds great!

Hi,

actually the support is already there, the problem is in the sparse triangular solver with an adjoint expression (see SparseLLT::solve()). I don't know if that will be easy to solve elegantly, but in the mean time I can try to evaluate the adjoint expression. I'll update this thread if I have any news.

ggael wrote:Hi,

actually the support is already there, the problem is in the sparse triangular solver with an adjoint expression (see SparseLLT::solve()). I don't know if that will be easy to solve elegantly, but in the mean time I can try to evaluate the adjoint expression. I'll update this thread if I have any news.

Hi,

Here is an update on using eigen/cholmod on complex sparse matrices. Here are the steps I followed, written in a newbie style...

0. My problem was to solve Ax=B, where A is a given sparse complex SPD matrix and B is a given dense complex vector. The idea was to use Eigen and a Cholesky decomposition plus direct solver. Typical size of the problem is above 1e6 (signal processing application)

1. I started by deriving sparse_cholesky.cpp (from bench directory), in order to understand a bit how to handle matrices, vectors, sparse objects, how to fill them, etc. Thanks to Gael/Benoit, native cholmod calls are also included in sparse_cholesky.cpp, so I could find my way to achieve what was needed using the native cholmod_* calls.

2. I've added MatrixBase::asCholmodDense() method in MatrixBaseAddons.h: it converts the rhs vector B (Matrix,...> ) into a cholmod_dense structure so I can pass it to the solver.

As explained in the doc, this line in your code is to tell MatrixBase to include some extensions

Code: Select all

#define EIGEN_MATRIXBASE_PLUGIN "MatrixBaseAddons.h"

MatrixBaseAddons.h: included by MatrixBase

Code: Select all

/**
 * asCholmodDense()
 *
 */
cholmod_dense asCholmodDense()
{
  cholmod_dense res;
  res.nrow   = rows();
  res.ncol   = cols();
  res.nzmax  = res.nrow * res.ncol;
  res.d      = res.nrow;
  res.x      = derived().data();
  // res.z unused

  if (ei_is_same_type::ret)
  {
    res.xtype = CHOLMOD_REAL;
    res.dtype = 1;
  }
  else if (ei_is_same_type::ret)
  {
    res.xtype = CHOLMOD_REAL;
    res.dtype = 0;
  }
  else if (ei_is_same_type >::ret)
  {
    res.xtype = CHOLMOD_COMPLEX;
    res.dtype = 1;
  }
  else if (ei_is_same_type >::ret)
  {
    res.xtype = CHOLMOD_COMPLEX;
    res.dtype = 0;
  }
  else
  {
    ei_assert(false && "Scalar type not supported by CHOLMOD");
  }

  return res;
}


3. I've added a SparseLLT::cholmodSolve(...B, ...&x) method in CholmodSupport.h. It simply uses the fast direct solver included in the cholmod library: cholmod_solve(...) rather than the triangular solver that is currently included in Eigen. I chose the cholmod solver for having complex support quickly (I know Gael is working on improving the triangular solver, but I had not much time, and it's interesting anyway).

In CholmodSupport.h:

Code: Select all

// in SparseLLT declaration
template
void cholmodSolve(MatrixBase& b, MatrixBase* res) const;

[.../...]

/**
 * cholmodSolve()
 *
 */
template
template
void SparseLLT::cholmodSolve(MatrixBase &b, MatrixBase* res) const
{
  /* FIXME: I don't know why, but converting cholmod_factor into a cholmod_sparse
     prevents cholmod_solve from working ("invalid xtype" error).
     So matrixL() is cannot be used prior to the cholmod solver.
   */
  //matrixL();
  cholmod_dense cdb = b.asCholmodDense();
  cholmod_dense* x;
  x = cholmod_solve(CHOLMOD_LDLt, m_cholmodFactor, &cdb, &m_cholmod);
  //cholmod_write_dense(stdout, x, 0, &m_cholmod);

  double* d = reinterpret_cast(x->x);
  for (int j=0; jderived().data()[j] = std::complex(d[2*j], d[2*j+1]);
  }

  cholmod_free_dense(&x, &m_cholmod);
}


Three remarks
a. I have made an ugly conversion from cholmod_dense to Matrix, ...> that needs to be templated. I'll be glad for any help on this!

b. Calling cholmod_factor_to_dense(...) prior to cholmod_solve prevents the latter to work properly ("invalid xtype" error): it looks like the m_cholmodFactor is subtly modified during cholmod_factor_to_dense(). Thus I cannot call matrixL() to return convert the m_choldmodFactor into a sparse matrix before any call to cholmodSolve... Nevermind, m_cholmodFactor is computed correctly during the cholesky factorization process, so I don't need to actually get it as a SparseMatrix.

c. To my understanding, my code actually does a LDLt decomposition (which is a non sense in the SparseLLT class!), maybe the CHOLMOD_LDLt may be changed to something else (CHOLMOD_A, etc) with the same results. I should test this soon.

I can now decompose and solve my complex sparse problems quickly, and still use Eigen's nice-nice API. My current tests show correct results and I'm happy with the speed (close to the native calls). I hope this can be of any help.

Cheers and many thanks to Eigen team!
Romain

hi, thanks for the explanations.

in short, from revision 902226:

* enable complex support for the CHOLMOD LLT backend using CHOLMOD's triangular solver
* quick fix for complex support in SparseLLT::solve

To answer your concerns:

a) I simply use a global function to internally convert a dense matrix to a cholmod matrix. I'll probably change the others for global functions too.

b) indeed, cholmod_solve does the convertion to a sparse matrix itself, and unfortunately, cholmod_factor_to_sparse modify the factor. that's really unfortunate, I'll handle that later

c) don't worry, there is no CHOLMOD_LLt because CHOLMOD knows if we performed a LLT or LDLT factorization, therefore using CHOLMOD_LDLt is fine. In SparseLLT::compute() cholmod is set to return a LLT factorization.