Factorization of ChunkSparseMatrix<double>

[Corbin Foucart]

Hello everyone,

I'm looking to LU factorize a ChunkSparseMatrix<double> instance (the same type as is used in step-51). However, it seems that ChunkSparseMatrix does not inherit from SparseMatrix, so I receive the predictable:

undefined reference to 'void dealii::SparseDirectUMFPACK::factorize<dealii::ChunkSparseMatrix<double> >(dealii::ChunkSparseMatrix<double> const&)'

I have a few ideas around this, but I'm not sure which is simplest to do (or whether I'm missing an easy templating/casting trick that would avoid the problem):

1. Implement a way to copy the ChunkSparseMatrix to a SparseMatrix as an intermediate before factorization
2. Implement SparseDirectUMFPACK::factorize(const Matrix &matrix) directly (this is perhaps tougher, since the function is templated by Matrix, and I can't seem to find the explicit template instantiations in the source code)
3. Do away with use of ChunkSparseMatrix entirely, and simply assemble into a SparseMatrix throughout (since I'm factorizing it anyway--perhaps this is the truest to the intentions of the library)

Any guidance would be appreciated!

[Wolfgang Bangerth]

Corbin,

> I'm looking to LU factorize a ChunkSparseMatrix<double> instance (the same
> type as is used in step-51). However, it seems that ChunkSparseMatrix does not
> inherit from SparseMatrix, so I receive the predictable:
>
> undefined reference to 'void
> dealii::SparseDirectUMFPACK::factorize<dealii::ChunkSparseMatrix<double>
> >(dealii::ChunkSparseMatrix<double> const&)'

Yes: This isn't implemented, because nobody so far saw a need to do so. The
ChunkSparseMatrix was intended to be used to make matrix-vector products
faster by using whole chunks that can be multiplied using BLAS operations.
Nobody thought of writing an interface to SparseDirectUMFPACK because if you
use the latter class, you are building a direct decomposition of the matrix
and the speed matrix-vector products are almost certainly no longer your
concern (because you're not using an iterative solver any more).

That doesn't mean that writing such an interface isn't useful: It's just that
nobody has seen a need for it in the past.

> I have a few ideas around this, but I'm not sure which is simplest to do (or
> whether I'm missing an easy templating/casting trick that would avoid the
> problem):
>

> 1. Implement a way to copy the ChunkSparseMatrix to a SparseMatrix as an
> intermediate before factorization

Yes, that would work.


> 2. Implement SparseDirectUMFPACK::factorize(const Matrix &matrix) directly

> (this is perhaps tougher, since the function is templated by Matrix, and I
> can't seem to find the explicit template instantiations in the source code)

If you look for it, there are functions that copy SparseMatrix and
BlockSparseMatrix objects to the necessary UMFPACK data structures. It
wouldn't be very difficult to add similar functions for ChunkSparseMatrix as
well, and we'd gladly take a patch of course!


> 3. Do away with use of ChunkSparseMatrix entirely, and simply assemble into a

> SparseMatrix throughout (since I'm factorizing it anyway--perhaps this is
> the truest to the intentions of the library)

That too would work. I think the question, if you intend to go with this
option, is why you chose ChunkSparseMatrix to begin with, and what effect it
would have to not use it any more.

Best
W.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/