Re: [eigen] new API for Cholmod

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Thu, Nov 4, 2010 at 7:10 PM, Bill Greene <w.h.greene@xxxxxxxxx> wrote:
> Yes, I think including ordering routines is quite important.
> For any significant size problem, ordering is almost always
> required to avoid very long computational times.
>
> I've thought a bit about the ordering issue recently. Here are
> a couple of things you might want to consider.
>
> 1. The ordering algorithm in a particular sparse matrix package
>     is likely to be useful with other sparse matrix factoring routines (i.e.
>     the CSparse AMD would be useful with MUMPS or SPOOLES, etc.)
>     The Metis ordering package is often very effective and is distributed or
>     recommended in several sparse matrix packages-- probably all
>     sparse backends would find this useful.
>
> 2. Users who are particularly interested in performance, often like to
>     experiment with different ordering algorithms to find the optimum
>     one for their problem.
>
> So its probably better to allow the user to separate the ordering step
> from the analysis and factorization steps.

yes, this is definitely something I have in mind too. Need more time...

>>This still require some work to make it compatible with any Eigen's
>>sparse representation (SparseMatrix, DynamicSparseMatrix, row-column
>>major, with only the lower, or upper or both triangular part etc.).
>
> Yes, I can see how this would be quite a bit of work. Did you consider
> just stopping with an error message if the matrix is not a SparseMatrix
> in the correct form-- the upper triangle for CSparse? As a user I don't
> think that would be so objectionable. I could then choose to convert the
> matrix to the correct format, myself, possibly avoiding the memory cost
> of having two copies of the matrix.

sure, this what I should have already done...

> On Wed, Nov 3, 2010 at 10:16 AM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
> wrote:
>>
>> Hi,
>>
>> regarding CSparse, since it is LGPL, last week, I've started a port of
>> its AMD ordering routine. This is the most interesting routine to port
>> since writing the decomposition on top of it is relatively easy.
>>
>> So basically, I finally have a usable LDLt decomposition which
>> automatically reduce the fill-in and does not require any external
>> library.
>>
>> This still require some work to make it compatible with any Eigen's
>> sparse representation (SparseMatrix, DynamicSparseMatrix, row-column
>> major, with only the lower, or upper or both triangular part etc.).
>>
>> gael
>>
>> On Sun, Oct 31, 2010 at 5:46 PM, Bill Greene <w.h.greene@xxxxxxxxx> wrote:
>> >>I'd be interested in CSparse for instance :
>> >>The good news is that it is CCS, but the Eigen doc[0] warns about a
>> >> difference:
>> >>"unlike several other sparse libraries (e.g., CSparse), the coefficients
>> >> of
>> >> our CCS matrices are always sorted. " Am I right to believe that this
>> >> won't
>> >> be a stumbling block ?
>> >
>> > There is no problem extracting the internal arrays of an Eigen CCS
>> > SparseMatrix and passing them to the CSparse routines. CSparse does not
>> > require the columns
>> > to be sorted by row-index but is quite happy to accept sparse matrices
>> > stored this way.
>> >
>> > CSparse (like CHOLMOD) *does* require the sparse matrices to be stored
>> > as
>> > CCS (rather than CRS). I downloaded the 10/29 zip file for Eigen3 and
>> > did
>> > not see in the
>> > code a test to make sure the matrix is CCS-- but I could easily have
>> > missed
>> > that.
>> >
>> > CSparse (unlike CHOLMOD) does not accept a sparse matrix with only the
>> > coefficients in the lower triangle stored.
>> > So something like
>> > Eigen::CSparseDecomposition<Eigen::SparseMatrix<double>,Eigen::Lower>
>> > chol;
>> > would need to be rejected.
>> >
>> > Which brings up a design question. Should the SparseMatrix class itself
>> > allow the user to set and check the characteristics of the matrix-- i.e.
>> > symmetric-lower,
>> > symmetric-upper, unsymmetric? I know the dense matrix classes treat
>> > these
>> > characteristics as a "view" of the matrix. I'm sure this issue has been
>> > discussed
>> > previously but I haven't seen this discussion, myself.
>> >
>> > Here is one other small observation about Eigen::CholmodDecomposition().
>> > The
>> > design lets you extract the success or failure of the last method via in
>> > the
>> > info()
>> > method. Internally, the class maintains m_analysisIsOk and
>> > m_factorizationIsOk values. Shouldn't the user be able to retrieve these
>> > anytime also?
>> >
>> > Bill Greene
>> >
>> >
>> > On Fri, Oct 29, 2010 at 12:23 PM,  <bernard.hugueney@xxxxxxxxxx> wrote:
>> >>
>> >> Hi,
>> >>>
>> >>> i,
>> >>>
>> >>> I've started to rework our sparse solvers/backends with Cholmod. The
>> >>> old API is still available.
>> >>> The first difference is that there is no SparssLLT or SparseLDLT or so
>> >>> base classes which was really confusing to use.
>> >>> Instead, you directly use the relevant wrapper class,
>> >>> CholmodDecomposition<MatrixType,UpLo>:
>> >>> [snip]
>> >>> You can also configure the underlying decomposition with:
>> >>> chol.setMode(mode);
>> >>> [snip]
>> >>> It is not that much changes but I think this approach is much cleaner
>> >>> and easier to use.
>> >>> I plane to write all solvers/wrappers using the same model, so it is
>> >>> perfect time to give your opinion on what would be the best API for
>> >>> them. The constraint is that we must be as close as possible to the
>> >>> dense API, but it is perfectly reasonable to add convenient methods
>> >>> with more runtime facilities....
>> >>
>> >> Cool, this is great news to me !
>> >> As previously discussed, I believe that the runtime of most (all ?)
>> >>  member function calls on sparse matrices and the lack of compile-time
>> >>  unrolling through fixed size makes static polymorphism pointless for
>> >>  such use cases, so I'll be glad to get the ease of use of runtime
>> >>  polymorphism instead.
>> >>
>> >> Wrt API, I think it should be possible to closely follow the dense API
>> >> with the only addition of your setMode(XXX) with XXX being the run-time
>> >> equivalent of the compile-time template parameters.
>> >>
>> >>> Then we could imagine to have a meta solver built on top of them with
>> >>> runtime registration of backends, priority rules, etc. but that's
>> >>> another story.
>> >>
>> >> The main difference is that the availability of backends could be
>> >> checked
>> >> at runtime instead of compile-time. That would be an asset when
>> >> shipping binaries using system libraries (e.g. libsparsesuite,
>> >> libumfpack on Debian). Do you have plans/ideas to support this ?
>> >>
>> >> Speaking about backends, I am wondering how hard is it to write a new
>> >> wrapper.
>> >> I'd be interested in CSparse for instance :
>> >> The good news is that it is CCS, but the Eigen doc[0] warns about a
>> >> difference:
>> >> "unlike several other sparse libraries (e.g., CSparse), the
>> >> coefficients
>> >> of
>> >> our CCS matrices are always sorted. " Am I right to believe that this
>> >> won't
>> >> be a stumbling block ?
>> >>
>> >> One of my reasons to be interested by CSparse is that the LU
>> >> decomposition
>> >> is rank revealing : is there another way to compute the rank with
>> >> currently
>> >> supported Eigen backends ?
>> >>
>> >> (I would also be insterested in SPRSBLKLLT)
>> >>
>> >> On the subject of backend choice, maybe it would be useful to add a few
>> >> references (I not qualified to evaluated the, but I found [1] anf [2]
>> >> to
>> >> be
>> >> interesting.), in the doc (maybe @[3]).
>> >>
>> >> Kind regards,
>> >>
>> >> Bernard
>> >> [0]
>> >> http://eigen.tuxfamily.org/index.php?title=SparseMatrix#Current_state
>> >> [1] http://www.cise.ufl.edu/research/sparse/codes/codes.pdf
>> >> [2] ftp://ftp.numerical.rl.ac.uk/pub/reports/ghsRAL200505.pdf
>> >> [3]
>> >>
>> >>
>> >> http://eigen.tuxfamily.org/dox-devel/TutorialSparse.html#TutorialSparseDirectSolvers
>> >>
>> >>
>> >>
>> >
>> >
>>
>>
>
>