[petsc-users] Trying to use MatSeq*AIJSetColumnIndices

[Matthew Knepley]

On Fri, Apr 17, 2015 at 2:11 PM, Mikael Öhman <micketeer at gmail.com> wrote:

> 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).
>

Preallocation is all you need. If you want to assure that the pattern is
stable under assembly, the best thing to do is insert zeros. That
is what we do by default for the DM matrices.


> 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?
>

Yes, that looks like a bug

   Matt


>
> Regards, Mikael
>



-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20150417/d5c4d4f4/attachment.html>