<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Wed, Dec 23, 2015 at 12:34 AM, Timothée Nicolas <span dir="ltr"><<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@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><div><div><div><div>Ho,<br></div><div><br>I found that MatSetUp or something like MatMPIAIJSetPreallocation have to be used beofre the MatSetValuesStencil. And if MatMPIAIJSetPreallocation then accordingly MatSetType must also be set to MPIAIJ. That solves the above problem. With this, and using my own routine to set the preallocation, as suggested by Dave, I can recover the same behavior for my matrices as when I set them with DMCreateMatrix, but in the <i><u>square</u> </i>case only (and by the way the whole allocation business seems to be faster in this case than using DMCreateMatrix).<br><br></div>However I still have problems in the rectangular case. The parallel case does not work and in one processor it works and gives the expected result only if I use ADD_VALUES in MatSetValuesStencil, even though the entries are supposed to be empty !<br><br></div>In the multiprocessor case, I get the "argument out of range" error like this:<br><br>1]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>[1]PETSC ERROR: Argument out of range<br>[1]PETSC ERROR: New nonzero at (3833,342512) caused a malloc<br>Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn off this check<br>[1]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>[1]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015 <br>[1]PETSC ERROR: ./miips on a arch-linux2-c-opt named helios91 by tnicolas Wed Dec 23 15:15:52 2015<br>[1]PETSC ERROR: Configure options --prefix=/csc/softs/anl/petsc-3.6.0/intel-15.0.0.090/bullxmpi-1.2.8.2/real --with-debugging=0 --with-x=0 --with-cc=mpicc --with-fc=mpif90 --with-cxx=mpicxx --with-fortran --known-mpi-shared-libraries=1 --with-scalar-type=real --with-precision=double --CFLAGS="-g -O3 -mavx -mkl" --CXXFLAGS="-g -O3 -mavx -mkl" --FFLAGS="-g -O3 -mavx -mkl"<br>[1]PETSC ERROR: #2265 MatSetValues_MPIAIJ() line 616 in /csc/releases/buildlog/anl/petsc-3.6.0/intel-15.0.0.090/bullxmpi-1.2.8.2/real/petsc-3.6.0/src/mat/impls/aij/mpi/mpiaij.c<br>[1]PETSC ERROR: #2266 MatSetValues() line 1175 in /csc/releases/buildlog/anl/petsc-3.6.0/intel-15.0.0.090/bullxmpi-1.2.8.2/real/petsc-3.6.0/src/mat/interface/matrix.c<br>[1]PETSC ERROR: #2267 MatSetValuesLocal() line 2023 in /csc/releases/buildlog/anl/petsc-3.6.0/intel-15.0.0.090/bullxmpi-1.2.8.2/real/petsc-3.6.0/src/mat/interface/matrix.c<br>[1]PETSC ERROR: #2268 MatSetValuesStencil() line 1423 in /csc/releases/buildlog/anl/petsc-3.6.0/intel-15.0.0.090/bullxmpi-1.2.8.2/real/petsc-3.6.0/src/mat/interface/matrix.c<br><br></div>This suggests my allocation is incorrect, however using the -info option, I see that<br><br>[1] MatAssemblyEnd_SeqAIJ(): Matrix size: 209088 X 348480; storage space: 20645544 unneeded,6170106 used<br>[1] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0<br>[1] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 31<br><br></div>Which seems to be in contradiction with the above error. So I am thinking it is a different problem.<br></div></div></div></div></div></div></blockquote><div><br></div><div>Make this problem simpler. That is the only way to tell what is going on. Do a 4x4 grid, so you have maybe 16 rows and can print everything out.</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><div><div><div><div></div>I am wondering if the problem does not come from the fact that ideally the matrix should be constituted from 3x5 blocks, but this does not seem to be allowed yet : <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIBAIJSetPreallocation.html" target="_blank">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIBAIJSetPreallocation.html</a> :<br><br>"<b>bs</b> - size of block, the blocks are ALWAYS square. One can use <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetBlockSizes.html#MatSetBlockSizes" target="_blank">MatSetBlockSizes</a>() to set a different row and column blocksize but the row
blocksize always defines the size of the blocks. The column blocksize sets the blocksize of the vectors obtained with <a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateVecs.html#MatCreateVecs" target="_blank">MatCreateVecs</a>()
"<br><br></div>So I cannot call MatSetBlockSizes(A,3,5,ierr)<br><br></div>and then MatMPIBAIJSetPreallocation because it would still assume square blocks (size 3) if I understand well. So I still don't see how to define a rectangular matrix which handily takes a vector from a DMDA and upon applying the operation, turns it to an other DMDA with same layout (except for number of DOF).<br><br></div>Best<br><br></div>Timothee<a name="151cd8c19f9ec943_MatMPIBAIJSetPreallocation"></a><div><div><div><div><div><div><div><br><br></div></div></div></div></div></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-12-23 6:05 GMT+09:00 Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
First run with valgrind and if that doesn't find anything then in the debugger and check the data structures where it crashes. There is nothing obviously wrong, some subtle mistake that can only be found with valgrind.<br>
<span><font color="#888888"><br>
Barry<br>
</font></span><div><div><br>
> On Dec 22, 2015, at 12:02 AM, Timothée Nicolas <<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>> wrote:<br>
><br>
> Hi,<br>
><br>
> I'd like to take advantage of the stencil structure of the DMDA to write into my matrix. I have already written all my routines so that I can use MatSetValuesStencil so it's most convenient for me, otherwise I will have to rethink all the structure.<br>
><br>
> Since I don't create my rectangular matrix with DMCreateMatrix but with MatCreate and MatSetSizes instead, the manual tells me I have to call first MatSetStencil and MatSetLocalToGlobalMapping before being allowed to call MatSetValuesStencil. So before going on the full example, I tried to do this on a simple square matrix but I can't get the thing to work. It must be an obvious mistake but can't find it. My code is<br>
><br>
> PetscErrorCode :: ierr<br>
> DM :: da<br>
> Mat :: A<br>
> PetscInt :: xm<br>
> ISLocalToGlobalMapping :: ltog<br>
> PetscInt :: dims(1), starts(1)<br>
> PetscInt :: dim<br>
> MatStencil :: row(4,1), col(4,1)<br>
> PetscScalar :: v<br>
><br>
> call DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,iten,ione,ione,PETSC_NULL_INTEGER,da,ierr)<br>
> call DMDAGetCorners(da,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &<br>
> & xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)<br>
><br>
> call MatCreate(PETSC_COMM_WORLD,A,ierr)<br>
><br>
> call MatSetSizes(A,xm,xm,PETSC_DETERMINE,PETSC_DETERMINE,ierr)<br>
><br>
> call DMGetLocalToGlobalMapping(da,ltog,ierr)<br>
> call MatSetLocalToGlobalMapping(A,ltog,ltog,ierr)<br>
> call DMDAGetGhostCorners(da,starts(1),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, &<br>
> & dims(1),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr)<br>
> dim=ione<br>
> call MatSetStencil(A,dim,dims,starts,ione,ierr)<br>
><br>
> row(MatStencil_i,1) = 0<br>
> col(MatStencil_i,1) = 0<br>
><br>
> call MatSetValuesStencil(A,ione,row,ione,col,zero,INSERT_VALUES,ierr)<br>
><br>
> But I keep getting a segmentation fault error.<br>
><br>
> [0]PETSC ERROR: ------------------------------------------------------------------------<br>
> [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range<br>
> [0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger<br>
> [0]PETSC ERROR: or see <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind</a><br>
> [0]PETSC ERROR: or try <a href="http://valgrind.org" rel="noreferrer" target="_blank">http://valgrind.org</a> on GNU/linux and Apple Mac OS X to find memory corruption errors<br>
> [0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, and run<br>
> [0]PETSC ERROR: to get more information on the crash.<br>
> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>
> [0]PETSC ERROR: Signal received<br>
> [0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>
> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015<br>
> [0]PETSC ERROR: ./miips on a arch-darwin-c-debug named <a href="http://iMac27Nicolas.nifs.ac.jp" rel="noreferrer" target="_blank">iMac27Nicolas.nifs.ac.jp</a> by timotheenicolas Tue Dec 22 14:57:46 2015<br>
> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich --with-debugging=no<br>
> [0]PETSC ERROR: #1 User provided function() line 0 in unknown file<br>
> application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0<br>
> [cli_0]: aborting job:<br>
> application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0<br>
><br>
> ===================================================================================<br>
> = BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES<br>
> = PID 9812 RUNNING AT <a href="http://iMac27Nicolas.nifs.ac.jp" rel="noreferrer" target="_blank">iMac27Nicolas.nifs.ac.jp</a><br>
> = EXIT CODE: 59<br>
> = CLEANING UP REMAINING PROCESSES<br>
> = YOU CAN IGNORE THE BELOW CLEANUP MESSAGES<br>
> ===================================================================================<br>
><br>
> If I replace everything with call DMCreateMatrix(da,A,ierr) it works fine.<br>
><br>
> What information did I miss from the manual ? What am I doing wrong ?<br>
><br>
> Best<br>
><br>
> Timothee<br>
><br>
><br>
><br>
><br>
><br>
> 2015-12-06 23:16 GMT+09:00 Timothée Nicolas <<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>>:<br>
> Wow, Great answer, thanks, I should probably be able to do it like this, I'll try my best !<br>
><br>
> I actually do want the assembled matrix now. I have coded the matrix-free version already and it works fine, but I want to use multigrid and I figured that in order to be be most efficient at the coarse level it is best to use -ksp_type preonly -pc_type lu, which requires a true matrix (the matrix I am talking about is the product between a matrix 3dof --> 5dof with a matrix 5dof --> 3dof, which ends up square). But I want to retain a memory light implementation so I use matrix-free formulation at all the levels but the last, which is also light by definition, so I should not suffer much in terms of memory.<br>
><br>
> Thanks<br>
><br>
> Timothée<br>
><br>
><br>
><br>
> 2015-12-06 22:29 GMT+09:00 Dave May <<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>>:<br>
><br>
><br>
> On 6 December 2015 at 11:19, Timothée Nicolas <<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>> wrote:<br>
> Hi,<br>
><br>
> The answer to the first question is yes. I was puzzled about the second part of the answer, but I realize I was not clear enough. I don't want to just transfer information between 2 DMDAs, I want to perform an operation on vectors with a rectangular matrix.<br>
><br>
> Okay - I misunderstood.<br>
><br>
><br>
> To be more explicit, the input vector of the operation is a vector with 3 components (hence dof=3), and the results of the operator is put in one vector equation + 2 scalar equation, which makes it dof = 5.<br>
><br>
> So, either I cannot fill my matrix (if the number of rows exceeds what is allowed by the DMDA), or I can fill my matrix but then the matrix and vectors are not compatible, as in this error :<br>
><br>
><br>
> Do you need to explicitly assemble this rectangular operator?<br>
> Expressing this operator via a matrix-free representation would be relatively straight forward.<br>
><br>
> However, since the DMDA's have the same parallel layout you are in good shape to define the rectangular matrix if you really want the assembled thing. Your row space would be defined by the DMDA with block size 5 and the column space would be defined by the DMDA with block size 3. For example the matrix would be created like in the following way (assuming a 2D DMDA)<br>
><br>
> PetscInt m,n,M,N;<br>
><br>
> DMDAGetInfo(dm3,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL);<br>
> nodes = M *N;<br>
> DMDAGetCorners(dm3,NULL,NULL,NULL,&m,&n,NULL);<br>
> lnodes = m * n;<br>
> MatCreate(PETSC_COMM_WORLD,&A);<br>
> MatSetSizes(A,lnodes*5,lnodes*3,nodes*5,nodes*3);<br>
><br>
> The sparsity pattern is defined by the connectivity of the DMDA points AND the different block sizes.<br>
> To define the non-zero structure, you need to traverse the part of the dmda defined by the natural indices given by si <= i < si+m, sj <= j < sj+n and check whether the the 5x3 block associated with each i,j is living within the diagonal or off-diagonal block of the matrix. For a given i,j the off-diagonal is defined by any i',j' associated your stencil which is not in the range si <= i' <= si + m, sj <= j' < sj + n. (e.g. it is a ghost node).<br>
><br>
> Your arrays counting the non-zero entries on the diagonal (nnz[]) and off diag (nnz_od[]) can be indexed using i + j * m. This corresponds to the indexing used by the DMDA. Your non-zero counting function would look something like the pseudo code below<br>
><br>
> PetscInt *nnz,*nnz_od;<br>
><br>
> PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz);<br>
> PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz_od);<br>
><br>
> PetscMemzero(nnz,sizeof(PetscInt)*m*n*5);<br>
> PetscMemzero(nnz_od,sizeof(PetscInt)*m*n*5);<br>
><br>
> for (j=0; j<n; j++) {<br>
> for (i=0; i<m; i++) {<br>
> PetscInt idx,d,cnti,cntg;<br>
> PetscBool is_ghost;<br>
><br>
> idx = i + j * m;<br>
><br>
> cnti = 0;<br>
> cntg = 0;<br>
> for (all_points_in_stencil) { /* User logic to define this loop */<br>
><br>
> is_ghost = PETSC_FALSE;<br>
> /*User logic goes here to test if stencil point is a ghost point */<br>
><br>
> /* accumulate counts for the stencil */<br>
> if (is_ghost) { cntg++; } /* values living on a neighbour rank */<br>
> else { cnti++; } /* values local to current rank */<br>
> }<br>
><br>
> /* assume all dofs in row and col space in block are coupled to each other */<br>
> for (d=0; d<5; d++) {<br>
> nnz[5*idx+d] = 3 * cnti;<br>
> nnz_od[5*idx+d] = 3 * cntg;<br>
> }<br>
><br>
> }<br>
> }<br>
><br>
> Thanks,<br>
> Dave<br>
><br>
><br>
><br>
> [0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>
> [0]PETSC ERROR: Nonconforming object sizes<br>
> [0]PETSC ERROR: Mat mat,Vec y: global dim 4900 2940<br>
> [0]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>
> [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015<br>
> [0]PETSC ERROR: ./miips on a arch-linux2-c-debug named Carl-9000 by timothee Sun Dec 6 19:17:20 2015<br>
> [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ --with-fc=gfortran --download-fblaslapack --download-mpich<br>
> [0]PETSC ERROR: #1 MatMult() line 2215 in /home/timothee/Documents/petsc-3.6.1/src/mat/interface/matrix.c<br>
><br>
> So I assume it means the matrix is assumed to be square. Is there a way to change this ?<br>
><br>
> Best<br>
><br>
> Timothée<br>
><br>
><br>
><br>
> 2015-12-05 20:41 GMT+09:00 Dave May <<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>>:<br>
><br>
><br>
> On 5 December 2015 at 11:15, Timothée Nicolas <<a href="mailto:timothee.nicolas@gmail.com" target="_blank">timothee.nicolas@gmail.com</a>> wrote:<br>
> Hi,<br>
><br>
> I need to create a rectangular matrix which takes vector structured by a DMDA and transforms them to vectors structured with an other DMDA. It does not look like there is a routine to do this since apparently DMCreateMatrix assumes square matrix. What is the clean way to do this ?<br>
><br>
> Do the DMDA's have exactly the same parallel layout and only differ by the number of DOFs?<br>
><br>
> If yes, and you only want to move entries between vectors associated with the two DMDA's, you can perform the transfer between vectors locally (rank-wise) in a very simple manner using the local number of points from DMDA and the known block sizes (DOFs) associated with your DMDA and the size of the subspace you ultimately want to construct.<br>
><br>
> Thanks,<br>
> Dave<br>
><br>
><br>
> The problem is I have a DMDA with 8 degrees of freedom (dof) and my matrix sends vectors living in a subspace with 3 dof to the other 5 dof subspace. I can define a square matrix living in the large, 8 dof subspace, but this is of course very inefficient in terms of memory and time.<br>
><br>
> Thanks<br>
><br>
> Best<br>
><br>
> Timothée<br>
><br>
><br>
><br>
><br>
><br>
<br>
</div></div></blockquote></div><br></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>