<div dir="auto">Look at the function that begins at line 1050 as the link should redirect to</div><div class="gmail_extra"><br><div class="gmail_quote">Il 10 Lug 2017 11:18 AM, "Florian Lindner" <<a href="mailto:mailinglists@xgm.de">mailinglists@xgm.de</a>> ha scritto:<br type="attribution"><blockquote class="quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hey Stefano,<br>
<div class="quoted-text"><br>
Am 10.07.2017 um 16:36 schrieb Stefano Zampini:<br>
> Florian,<br>
><br>
> Perhaps you might want to take a look at how this is done for MatIS<br>
><br>
> <a href="https://bitbucket.org/petsc/petsc/src/f9d5775f43f69cbce5a7014a6ce3b24cc0e1214a/src/mat/impls/is/matis.c?at=master&fileviewer=file-view-default#matis.c-1050" rel="noreferrer" target="_blank">https://bitbucket.org/petsc/<wbr>petsc/src/<wbr>f9d5775f43f69cbce5a7014a6ce3b2<wbr>4cc0e1214a/src/mat/impls/is/<wbr>matis.c?at=master&fileviewer=<wbr>file-view-default#matis.c-1050</a><br>
<br>
</div>to what piece of code your are specifically refering to?<br>
<br>
Line 655, MatSetValuesLocal_SubMat_IS?<br>
<br>
and I should use ISLocalToGlobalMappingApply?<br>
<br>
but this, if I understand correctly, would involve such a call for every index.<br>
<br>
Right now, my preallocation loop looks like:<br>
<br>
<br>
    PetscInt d_nnz[local_rows], o_nnz[local_rows];<br>
<div class="quoted-text"><br>
    PetscInt mapSize;<br>
    ISLocalToGlobalMappingGetSize(<wbr>ISmapping, &mapSize);<br>
    const PetscInt *mapIndizes;<br>
</div>    ISLocalToGlobalMappingGetIndic<wbr>es(ISmapping, &mapIndizes); // is called just once<br>
<br>
    for (int row = row_range_start; row < row_range_end; row++) {<br>
      PetscInt relative_row = row - row_range_start;<br>
      d_nnz[relative_row] = 0;<br>
      o_nnz[relative_row] = 0;<br>
      int colStart = row - nz_ratio/2 * size;<br>
      int colEnd   = row + nz_ratio/2 * size;<br>
      colStart = (colStart < 0)    ? 0    : colStart;<br>
      colEnd   = (colEnd   > size) ? size : colEnd;<br>
      for (int col = 0; col < global_cols; col++) {<br>
        if (col >= colStart and col < colEnd) {<br>
<br>
          int petsc_col = mapIndizes[col]; // should be cheap<br>
          if (petsc_col >= col_range_start and petsc_col < col_range_end) {<br>
<div class="quoted-text">            d_nnz[relative_row]++;<br>
          } else {<br>
            o_nnz[relative_row]++;<br>
          }<br>
<br>
        }<br>
      }<br>
    }<br>
</div>    ierr = MatMPIAIJSetPreallocation(A, 0, d_nnz, 0, o_nnz); CHKERRQ(ierr);<br>
<br>
Best,<br>
Florian<br>
<br>
><br>
> Stefano<br>
<div class="elided-text">><br>
> Il 10 Lug 2017 10:23 AM, "Florian Lindner" <<a href="mailto:mailinglists@xgm.de">mailinglists@xgm.de</a> <mailto:<a href="mailto:mailinglists@xgm.de">mailinglists@xgm.de</a>>> ha scritto:<br>
><br>
>     Hey,<br>
><br>
>     one more question about preallocation:<br>
><br>
>     I can determine if a column index is diagonal or off-diagonal using that code<br>
><br>
>     if (col >= col_range_start and col < col_range_end)<br>
>         d_nnz[relative_row]++;<br>
>     else<br>
>         o_nnz[relative_row]++;<br>
><br>
><br>
>     My code, however uses index sets from which a ISlocalToGlobalMapping created:<br>
><br>
>       // Creates a mapping from permutation and use that for the cols<br>
>       ierr = ISCreateGeneral(PETSC_COMM_<wbr>WORLD, local_cols, permutation.data(), PETSC_COPY_VALUES, &ISlocal); CHKERRV(ierr);<br>
>       ierr = ISSetPermutation(ISlocal); CHKERRV(ierr);<br>
>       ierr = ISAllGather(ISlocal, &ISglobal); CHKERRV(ierr); // Gather the IS from all processors<br>
>       ierr = ISLocalToGlobalMappingCreateIS<wbr>(ISglobal, &ISmapping); CHKERRV(ierr); // Make it a mapping<br>
><br>
>       // Create an identity mapping and use that for the rows of A.<br>
>       ierr = ISCreateStride(PETSC_COMM_<wbr>WORLD, local_rows, row_range_start, 1, &ISidentity); CHKERRV(ierr);<br>
>       ierr = ISSetIdentity(ISidentity); CHKERRV(ierr);<br>
>       ierr = ISAllGather(ISidentity, &ISidentityGlobal); CHKERRV(ierr);<br>
>       ierr = ISLocalToGlobalMappingCreateIS<wbr>(ISidentityGlobal, &ISidentityMapping); CHKERRV(ierr);<br>
><br>
>       ierr = MatSetLocalToGlobalMapping(A, ISidentityMapping, ISmapping); CHKERRV(ierr);<br>
><br>
>     since SetPreallocation routines define the diagonal / off-diagonal blocks from the petsc ordering, I have to translate<br>
>     the col to a petsc_col.<br>
><br>
>     What is the best / fastest way to do that?<br>
><br>
>     Is that the way to go?<br>
><br>
>       PetscInt mapSize;<br>
>       ISLocalToGlobalMappingGetSize(<wbr>ISmapping, &mapSize);<br>
>       const PetscInt *mapIndizes;<br>
>       ISLocalToGlobalMappingGetIndic<wbr>es(ISmapping, &mapIndizes);<br>
><br>
>     Thanks,<br>
>     Florian<br>
><br>
><br>
><br>
>     Am 07.07.2017 um 17:31 schrieb Florian Lindner:<br>
>     > Hello,<br>
>     ><br>
>     > I'm having some struggle understanding the preallocation for MPIAIJ matrices, especially when a value is in<br>
>     off-diagonal<br>
>     > vs. diagonal block.<br>
>     ><br>
>     > The small example program is at <a href="https://pastebin.com/67dXnGm3" rel="noreferrer" target="_blank">https://pastebin.com/67dXnGm3</a><br>
>     ><br>
>     > In general it should be parallel, but right now I just run it in serial.<br>
>     ><br>
>     > According to my understanding of<br>
>     ><br>
>     > <a href="http://www.mcs.anl.gov/petsc/petsc-3.7/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-3.7/docs/manualpages/<wbr>Mat/MatMPIAIJSetPreallocation.<wbr>html</a><br>
>     <<a href="http://www.mcs.anl.gov/petsc/petsc-3.7/docs/manualpages/Mat/MatMPIAIJSetPreallocation.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/<wbr>petsc-3.7/docs/manualpages/<wbr>Mat/MatMPIAIJSetPreallocation.<wbr>html</a>><br>
>     ><br>
>     > a entry is in the diagonal submatrix, if its row is in the OwnershipRange and its column is in OwnershipRangeColumn.<br>
>     > That also means that in a serial run, there is only a diagonal submatrix.<br>
>     ><br>
>     > However, having MAT_NEW_NONZERO_ALLOCATION_ERR set, I get an error when<br>
>     ><br>
>     > Inserting 6 elements in row 2, though I have exactly<br>
>     ><br>
>     > 2 o_nnz = 0, d_nnz = 6 (means 6 elements allocated in the diagonal submatrix of row 2)<br>
>     ><br>
>     > Error is:<br>
>     ><br>
>     > [0]PETSC ERROR: Argument out of range<br>
>     > [0]PETSC ERROR: New nonzero at (2,5) caused a malloc<br>
>     > Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_<wbr>ERR, PETSC_FALSE) to turn off this check<br>
>     ><br>
>     ><br>
>     > What is wrong with my understanding?<br>
>     ><br>
>     > Thanks,<br>
>     > Florian<br>
>     ><br>
><br>
</div></blockquote></div><br></div>