<div dir="ltr"><div><div><div><div><div>In OOFEM, we do a "dry assembly" from which we compute the preallocation ( e.g. MatSeqSBAIJSetPreallocation ). <br><br></div>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.<br></div></div><br></div>First, I'd like to make sure I'm using these in the correct way. I do the following order of calls:<br><br>        MatCreate(PETSC_COMM_SELF, & mtrx);<br>        MatSetSizes(mtrx, leqs, leqs, geqs, geqs);<br>        MatSetType(mtrx, MATSEQSBAIJ);<br>        MatSeqSBAIJSetPreallocation( mtrx, blocksize, 0, d_nnz_sym.givePointer() );<br></div><div>        // Trying to add this for faster assembly:<br></div><div>        MatSeqSBAIJSetColumnIndices( mtrx, indices.givePointer() );<br>        MatSetUp(mtrx);<br>        MatSetOption(mtrx, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);<br><br></div><div>I run a fairly large finite element analysis (600k DOFs) through VTUNE to profile.<br><br><br><br>I have two issues:<br>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.<br></div><div>MatAssembly isn't terribly slow for me, but it does take 1/5th of the time that KSPSolve takes (which seems a big high).<br><br></div><div>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<span class=""> "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" at line 1067:<br></span><pre><span class="">for</span> <span class="">(</span><span class="">i</span><span class="">=</span><span class="">0</span><span class="">;</span> <span class="">i</span><span class=""><</span><span class="">n</span><span class="">;</span> <span class="">i</span><span class="">++</span><span class="">)</span> <span class="">baij</span><span class="">-></span><span class="">ilen</span><span class="">[</span><span class="">i</span><span class="">]</span> <span class="">=</span> <span class="">baij</span><span class="">-></span><span class="">imax</span><span class="">[</span><span class="">i</span><span class="">];</span></pre>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.<br></div><div>A bug perhaps?<br></div><div><br><br></div><div>Regards, Mikael<br></div></div>