<div dir="ltr"><div class="gmail_default" style="font-size:small">I worked on the assumptions in my previous email and I at least partially implemented the function to assign the couplings.  For row 0, which is the redundant field, I set dnz[0] to end-start, and onz[0] to the size of the matrix minus dnz[0].  For all other rows, I just increment the existing values of dnz[i] and onz[i], since the coupling to the redundant field adds one extra element beyond what's allocated for the DMDA stencil.</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">I see in the source that the FormCoupleLocations function is called once if the DM has PREALLOC_ONLY set to true, but twice otherwise.  I assume that the second call is for setting the nonzero structure.  Do I need to do this?</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">In any case, something is still not right.  Even with the extra elements preallocated, the first assembly of the matrix is very slow.  I ran a test problem on a single process with -info, and got this:</div><div class="gmail_default" style="font-size:small"><br></div><div class="gmail_default" style="font-size:small">







<p class="gmail-p1"><span class="gmail-s1">0] MatAssemblyEnd_SeqAIJ(): Matrix size: 1 X 1; storage space: 0 unneeded,1 used</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 1</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 1) < 0.6. Do not use CompressedRow routines.</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatSeqAIJCheckInode(): Found 1 nodes out of 1 rows. Not using Inode routines</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 9009 X 9009; storage space: 0 unneeded,629703 used</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 9009) < 0.6. Do not use CompressedRow routines.</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatSeqAIJCheckInode(): Found 3003 nodes of 9009. Limit used: 5. Using Inode routines</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 1 X 1; storage space: 0 unneeded,1 used</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 1</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 1) < 0.6. Do not use CompressedRow routines.</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatSeqAIJCheckInode(): Found 1 nodes out of 1 rows. Not using Inode routines</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 9009 X 9009; storage space: 0 unneeded,629703 used</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 9009) < 0.6. Do not use CompressedRow routines.</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatSeqAIJCheckInode(): Found 3003 nodes of 9009. Limit used: 5. Using Inode routines</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 9010 X 9010; storage space: 18018 unneeded,629704 used</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 81</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 9010) < 0.6. Do not use CompressedRow routines.</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatSeqAIJCheckInode(): Found 3004 nodes of 9010. Limit used: 5. Using Inode routines</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446774208 30818112</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] DMGetDMSNES(): Creating new DMSNES</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] DMGetDMKSP(): Creating new DMKSP</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1"><span class="gmail-Apple-converted-space">  </span>0 SNES Function norm 2.302381528359e+00<span class="gmail-Apple-converted-space"> </span></span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] PetscCommDuplicate(): Using internal PETSc communicator 139995446773056 30194064</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 9010 X 9010; storage space: 126132 unneeded,647722 used</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 9610</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 9010</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 9010) < 0.6. Do not use CompressedRow routines.</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatSeqAIJCheckInode(): Found 3004 nodes of 9010. Limit used: 5. Using Inode routines</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Matrix size: 9010 X 9010; storage space: 0 unneeded,647722 used</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 9010</span></p>
<p class="gmail-p1"><span class="gmail-s1">[0] MatCheckCompressedRow(): Found the ratio (num_zerorows 0)/(num_localrows 9010) < 0.6. Do not use CompressedRow routines.</span></p><p class="gmail-p1"><br></p><p class="gmail-p1">The 9610 mallocs during MatSetValues seem suspicious and are probably what's taking so long with larger problems.  601 of them are apparently in the call to set Jbh, and 9009 are in the call to set Jhb (b is the redundant field, h is the DMDA field).  If I run with more than one process, I get a segfault when a process which has rank greater than 0 sets dnz or onz in the FormCoupleLocations call.<br></p></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Oct 31, 2017 at 10:40 PM, zakaryah . <span dir="ltr"><<a href="mailto:zakaryah@gmail.com" target="_blank">zakaryah@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 class="gmail_default" style="font-size:small">Thanks Barry, that looks like exactly what I need.  I'm looking at pack.c and packm.c and I want to check my understanding of what my coupling function should do.  The relevant line in <strong style="color:rgb(0,0,0)"><font color="#4169E1">DMCreateMatrix_Composite_<wbr>AIJ </font></strong>seems to be:</div><div class="gmail_default" style="font-size:small"><span style="color:rgb(0,0,0)"><br></span></div><div class="gmail_default" style="font-size:small"><span style="color:rgb(0,0,0)">(*com->FormCoupleLocations)(<wbr>dm,NULL,dnz,onz,__rstart,__<wbr>nrows,__start,__end);</span><br></div><div class="gmail_default" style="font-size:small"><span style="color:rgb(0,0,0)"><br></span></div><div class="gmail_default" style="font-size:small"><span style="color:rgb(0,0,0)">and I infer that dnz and onz are the number of nonzero elements in the diagonal and off-diagonal submatrices, for each row of the DMComposite matrix.  I suppose I can just set each of these in a for loop, but can I use the arguments to FormCoupleLocations as the range for the loop?  Which ones - __rstart to __rstart+__nrows?  How can I determine the number of rows on each processor from within the function that I pass?  From the preallocation macros it looks like __start to __end describe the range of the columns of the diagonal submatrix - is that right?  It looks like the ranges will be specific to each processor.  Do I just set the values in dnz and onz, or do I need to reduce them?</span></div><div class="gmail_default" style="font-size:small"><span style="color:rgb(0,0,0)"><br></span></div><div class="gmail_default" style="font-size:small"><span style="color:rgb(0,0,0)">Thanks for all the help!  Maybe if I get things working I can carve out the core of the code to make an example program for DMRedundant/Composite.</span></div></div>
</blockquote></div><br></div>