<div dir="ltr">On Fri, Mar 1, 2013 at 12:50 PM, Ramsey, James J CIV (US) <span dir="ltr"><<a href="mailto:james.j.ramsey14.civ@mail.mil" target="_blank">james.j.ramsey14.civ@mail.mil</a>></span> wrote:<br><div class="gmail_extra">
<div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">I'm trying to get a better preallocation scheme. I noticed an earlier e-mail <<a href="https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2010-September/006886.html" target="_blank">https://lists.mcs.anl.gov/mailman/htdig/petsc-users/2010-September/006886.html</a>> mentioning that "for finite elements with a normal cell division or finite differences with a normal vertex division, you need only count the nonzeros for rows you own since the partitions overlap on the boundary." Unfortunately, this doesn't quite work for me. Here's the latest attempt at a routine, with some details skipped:<br>
</blockquote><div><br></div><div style>I understand the code you have below, but we have chosen to do it a complementary way.</div><div style>The way I understand FEM allocation, schematically, you</div><div style><br></div>
<div style>  Loop over all global unknowns u_i (coefficient of phi_i)</div><div style>    Loop over global unknowns u_j (coefficient of phi_j) which have nonzero overlap with the function \phi_i</div><div style>      Increment dnnz or onnz</div>
<div style><br></div><div style>We have code to do this in PETSc right now. You need to specify your topology using DMPlex, and then just use DMCreateMatrix().</div><div style>There is no getting around communication or ghosting when doing this operation.</div>
<div style><br></div><div style>   Matt</div><div style> <br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
std::vector<PetscInt> d_nnz(nLocalRows), o_nnz(nLocalRows);<br>
std::vector<std::set<PetscInt> > d_nzInds(nLocalRows), o_nzInds(nLocalRows);<br>
<br>
for ( ... loop over elements ...) {<br>
<br>
   wrapNodeIds(*mesh, nIdVec);  // Accounts for periodic boundary conditions<br>
   mesh->nodes2globalDOF(nIdVec, globalDOFsPerElem);<br>
<br>
   std::vector<PetscInt>::iterator gdBegin = globalDOFsPerElem.begin();<br>
   std::vector<PetscInt>::iterator gdEnd = globalDOFsPerElem.end();<br>
<br>
   for (std::vector<PetscInt>::iterator itrI = gdBegin; itrI != gdEnd; ++itrI) {<br>
<br>
       PetscInt I = *itrI;<br>
<br>
       if ((I >= minDOFId) && (I <= maxDOFId)) {<br>
          PetscInt Ioffset = I - minDOFId;<br>
<br>
          for (std::vector<PetscInt>::iterator itrJ = gdBegin; itrJ != gdEnd; ++itrJ) {<br>
<br>
            PetscInt J = *itrJ;<br>
<br>
            if ((J >= minDOFId) && (J <= maxDOFId)) {<br>
              // J will not be inserted if it is already in d_nzInds[Ioffset]. This avoids duplicates.<br>
              d_nzInds[Ioffset].insert(J);<br>
            }<br>
            else {<br>
              o_nzInds[Ioffset].insert(J);<br>
            }<br>
<br>
          }<br>
<br>
     }<br>
}<br>
<br>
for (PetscInt i = 0; i < nLocalRows; ++i) {<br>
      d_nnz[i] = d_nzInds[i].size();<br>
      o_nnz[i] = o_nzInds[i].size();<br>
}<br>
<br>
For nodes in the interior, this seems to work well, but at the boundaries of the domain, the number of non-zeros is underestimated. It seems like parallel communication is unavoidable here, but I'm not sure how to do it, and I'm unsure if I really want to mess with relatively low-level MPI routines to get what I want.<br>

<br>
Any ideas?</blockquote></div><br><br clear="all"><div><br></div>-- <br>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>