[petsc-users] Trying to use MatSeq*AIJSetColumnIndices

[Mikael Öhman]

In OOFEM, we do a "dry assembly" from which we compute the preallocation (
e.g. MatSeqSBAIJSetPreallocation ).

Since I already compute the column index in order to compute the "nnz"
input for preallocation, I figured I'd try to use the
"MatSeq*AIJSetColumnIndices" functions as well.

First, I'd like to make sure I'm using these in the correct way. I do the
following order of calls:

        MatCreate(PETSC_COMM_SELF, & mtrx);
        MatSetSizes(mtrx, leqs, leqs, geqs, geqs);
        MatSetType(mtrx, MATSEQSBAIJ);
        MatSeqSBAIJSetPreallocation( mtrx, blocksize, 0,
d_nnz_sym.givePointer() );
        // Trying to add this for faster assembly:
        MatSeqSBAIJSetColumnIndices( mtrx, indices.givePointer() );
        MatSetUp(mtrx);
        MatSetOption(mtrx, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);

I run a fairly large finite element analysis (600k DOFs) through VTUNE to
profile.



I have two issues:
1. There is no difference in the MatAssembly routine. I can't tell if I'm
using the routine wrong, or if preallocation is really all you need to
achieve maximum performance.
MatAssembly isn't terribly slow for me, but it does take 1/5th of the time
that KSPSolve takes (which seems a big high).

2. Setting the column indices for AIJ, BAIJ works without errors, but when
I try with SBAIJ with a blocksize > 1, then it segfaults in
"MatSeqSBAIJSetColumnIndices_SeqSBAIJ" at line 1067:

for (i=0; i<n; i++) baij->ilen[i] = baij->imax[i];

where the loop (n = mat->cmap->n) goes up to the number of columns (as
opposed to columns/blocksize). The corresponding BAIJ implementation uses
"baij->mbs" as the loop limit.
A bug perhaps?


Regards, Mikael
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150417/a11d7b68/attachment.html>