[petsc-users] MatCreateMPIAIJWithSplitArrays for symmetric matrix?

[Daniel Langr]

Dne 8.11.2010 13:31, Jed Brown napsal(a):
> On Mon, Nov 8, 2010 at 13:07, Daniel Langr <daniel.langr at gmail.com
> <mailto:daniel.langr at gmail.com>> wrote:
>
>     Ok, that's a good idea and I will try it. I can estimate total
>     number of nonzeros (limited by the amount of memory), but I can
>     hardly estimate number of nonzeros per particular rows (d_nnz,
>     o_nnz). So I need to use d_nz and o_nz, right? And for symmetric
>     matrix, what do d_nz and o_nz mean? Is it number of nonzeros for
>     whole matrix row, or just upper triangular part? (I don't see this
>     in manual page.)
>
>
> It is only for the upper triangular part.  If you have no way to guess
> the number per row, then just set {d,o}_nz, but note that this is
> pessimistic because the last rows of the matrix necessarily have fewer
> nonzeros in the upper triangular part.
>
>     (That was my point, through MatMPISBAIJSetPreallocation I need to
>     estimate number of nonzeros per row. Not just number of nonzeros of
>     whole local part. If I construct CSR arrays directly, I don't need
>     this row-estimates.)
>
>
> Here's a little secret: if you assemble the matrix row-by-row, in-order,
> then you need only estimate the average number of nonzeros per row.  The
> implementation is not stupid, so early rows can blow over their initial
> estimates as long as later rows make up for it by being under the
> estimate.  For example, with a 5-point Laplacian, you can
> MatSeqSBAIJSetPreallocation(A,1,3,PETSC_NULL) and you will see something
> like
>
> Matrix Object:
>    type: seqsbaij
>    rows=400, cols=400
>    total: nonzeros=1160, allocated nonzeros=1200
>    total number of mallocs used during MatSetValues calls =0
>        block size is 1
>
> So no mallocs even though the early rows have more than 3 nonzeros in
> the upper triangular part.  If you assemble out-of-order, or in a
> scattered fashion (such as for a finite element method) then this will
> be costly.  Does this help?

Great, helps much, that is something I was hoping for :). As for now, we 
can generate matrix in-order. In the future it seems we will need to 
generate matrix in "transposed" order, e.g. row-by-row order for the 
lower triangular part. But then I suppose we will need to use 
MatMultTranspose() or MatCreateTranspose() if the performance won't 
suffer much (I have not tried it yet, so I have no idea about it for now.)

Thanks a lot
Daniel

> Jed