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 <[email protected]> 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 <[email protected]> wrote:
> Barry Smith <[email protected]> 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
