Hi,

I wanted to create, fill and use symmetric sparse matrices.

Ideally I wanted to create them using

Code: Select all

class SymSparseMatrix : public Eigen::SparseMatrix<SelfAdjoint|LowerTriangular>
{
};


My fill commands ensure that only the lowertriangular part will be filled. usually I convert a selfadjoint lowertriangular DynamicSparseMatrix to a SparseMatrix.

So far, most things work fine: arithmetic operators with itself and other matrix types, sum(), norm(), etc.

Now I wanted to use it in conjunction with external bindings and I learned that my matrices work differently than expected.

When filling my selfadjoint matrix, what is actually happening???

I expected a lowertriangular memory layout. Is it possible that it creates a full matrix from two elements, if only one element is given in the lowertriangular part?

For example, when converting a general sparsematrix to my symmetric sparse matrix type, I expected that the number of nonzeros decreased - but it was the same.

Any help?

My idea is:
I did not want to define my symmetric matrix class as a general Eigen::SparseMatrix and overload all arithmetic operators by

Code: Select all

   SparseMatrix SymSparseMatrix::operator*(const Eigen::MatrixXd &x)const
   {
      return this->Base::marked< Eigen::ColMajor|Eigen::SelfAdjoint|Eigen::LowerTriangular >().operator*(x);



Isn't there a way to do this in a more elegant way?
(I am using Eigen 2.0.x)

Ok,

I defined my own sparse matrix class and redefined most operators that I need for now. At least this helps in the moment. It would be interesting, though, how Eigen 3.0 will be designed in this regard.

In Eigen 2.0.x it was possible to do:

(1) SparseMatrix<SelfAdjoint|LowerTriangular> matrix;

which I thought will give me the desired matrix (the naming scheme would have been intuitive - there is, however, no doc). But at least the matrix behaved correctly in the sense, that a full matrix was created from lower triangular data and all operators were applied as in the default sparse matrix case.

(2) y = matrix.marked<SelfAdjoint|LowerTrianular>().operator*(x);

which should interpret the contents of the matrix as being symmetric and given as lower triangular (in practice: It does not work correctly with 2.0.12)

From Eigen-devel doc I read about plans for a class type
(3) SelfAdjointView<SparseMatrix, LowerTriangular> matrix;

Is this the prototype of a full-featured sparse matrix class which can be used as any other sparse matrix type or is are any special handlings required?

Thanks in advance,
Sebastian

In the devel branch I'm using the same API than for dense objects, i.e., using a SelfAdjointView class and its respective .selfadjointView<...> function.

The SelfAdjointView class will only expose a limited but relevant API like products and solvers.