[petsc-users] MatMPIBAIJSetPreallocationCSR i and j indices

[Barry Smith]

   Matteo,

      This is 20+ years of PETSc history, it centers around the imperfect*, yet extremely important, principle of "data encapsulation". In our case, to first order, this means that scientific simulation codes that need solvers should NOT have to know about the specific storage formats formats of the sparse matrices and should not be written around any particular storage format. The compressed sparse row (CSR) and block CSR formats are specific storage formats, useful for some operations (like matrix vector product) but cumbersome for other operations (like inserting a value into the matrix). If you are calling MatCreateMPIBAIJWithArrays() or MatMPIBAIJSetPreallocationCSR() this means that YOUR application code is directly storing the sparse matrix in a particular data format; hence no data encapsulation. 

    The only reason the routines MatMPIBAIJSetPreallocationCSR() and friends and relations exist in PETSc is because some potential users begged us to add them to PETSc so they could take an old code that already built the matrix in CSR format and use the PETSc solvers without having to rewrite a bunch of their code. Being nice guys (except to each other) and desperate for users :-) we provided these convenience routines.

   We don't want people who are writing new codes, or seriously modifying old codes, to use these routines because we don't think that users should be building matrices directly in  (block) CSR, or any other specific, format. Rather we think they should use one of the many MatSetValues friends and relations to pass the matrix values to PETSc, or doing matrix-free etc. Partly this ideological, but it is also practical, I claim building parallel matrices from, say finite element methods, with CSR directly is painful and not needed since PETSc can handle the "merging" of the element stiffness matrices into the "global matrix" for you, we already wrote the code, why should each person write the code themselves? In addition, if another format is best for implementation (say in multicore or on GPUs) by abstracting away the specific storage from the user code, they can switch to the new format without changing their code. If you are going to write this CSR code, why not write your own GMRES, your own Newton's method, your own ODE integrator? We believe in numerical libraries to gather up all the code that yes, in theory, each person can write themselves, but write it only once in the best way we know how and continue to improve it based on input from all new users. This allows people to make far more scientific simulation progress instead of spending a lot of time reinventing the wheel. Now since PETSc is a library, you can choose to use whatever part you want and to write the other parts yourself, so you can do your own assembly to block CSR yourself if you want, it is just something we don't think would be the most productive use of your time.

  Barry

* No computer science abstraction/idea is perfect, but some are useful.

On Sep 29, 2013, at 5:22 PM, Matteo Parsani <parsani.matteo at gmail.com> wrote:

> Hi Barry and Jed,
> just a curiosity and maybe a good point to learn something:
> 
> 1- why are those inconvenient "convenience routines" ?
> 
> 2 - why those routines should not be used when writing a fresh code or a large re-factoring?
> 
> 
> Thanks in advance.
> 
> 
> 
> On Sun, Sep 29, 2013 at 3:52 PM, Jed Brown <jedbrown at mcs.anl.gov> wrote:
> Barry Smith <bsmith at mcs.anl.gov> writes:
> 
> >    We can do all this, but how much effort do we really want to make
> >    for supporting these inconvenient "convenience routines"?
> 
> Not much, but this is a tiny amount of code.
> 
> >    BTW: We should really label these function with
> >
> >    Level: Recommended only for legacy applications being quickly ported to PETSc
> >
> >     I really only want people who already have a working non-PETSc
> >     code that uses CSR directly already to use these routines so they
> >     can quickly try out PETSc, I never want someone writing fresh code
> >     or doing a large refactorization to EVER use these routines.
> 
> Agreed.
> 
> >> On second thought, why don't we just document that
> >> MatMPIBAIJSetPreallocationCSR respects MAT_ROW_ORIENTED within each
> >> block and let the user choose which they want to use?
> >
> >     Also there is no way to get this flag into MatCreateMPIBAIJWithArrays()
> 
> The user calls MatSetOption before making that call and changes it
> afterward if they want.
> 
> 
> 
> -- 
> Matteo