<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Fri, Apr 17, 2015 at 2:11 PM, Mikael Öhman <span dir="ltr"><<a href="mailto:micketeer@gmail.com" target="_blank">micketeer@gmail.com</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><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></div></div></blockquote><div><br></div><div>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</div><div>is what we do by default for the DM matrices.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div></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> "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" at line 1067:<br></span><pre><span>for</span> <span>(</span><span>i</span><span>=</span><span>0</span><span>;</span> <span>i</span><span><</span><span>n</span><span>;</span> <span>i</span><span>++</span><span>)</span> <span>baij</span><span>-></span><span>ilen</span><span>[</span><span>i</span><span>]</span> <span>=</span> <span>baij</span><span>-></span><span>imax</span><span>[</span><span>i</span><span>];</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?</div></div></blockquote><div><br></div><div>Yes, that looks like a bug</div><div><br></div><div>   Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div><br></div><div>Regards, Mikael<br></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>