On Wed, Dec 23, 2015 at 12:34 AM, Timothée Nicolas < [email protected]> wrote:
> Ho, > > 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 *square *case > only (and by the way the whole allocation business seems to be faster in > this case than using DMCreateMatrix). > > 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 ! > > In the multiprocessor case, I get the "argument out of range" error like > this: > > 1]PETSC ERROR: --------------------- Error Message > -------------------------------------------------------------- > [1]PETSC ERROR: Argument out of range > [1]PETSC ERROR: New nonzero at (3833,342512) caused a malloc > Use MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE) to turn > off this check > [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html > for trouble shooting. > [1]PETSC ERROR: Petsc Release Version 3.6.0, Jun, 09, 2015 > [1]PETSC ERROR: ./miips on a arch-linux2-c-opt named helios91 by tnicolas > Wed Dec 23 15:15:52 2015 > [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" > [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 > [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 > [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 > [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 > > This suggests my allocation is incorrect, however using the -info option, > I see that > > [1] MatAssemblyEnd_SeqAIJ(): Matrix size: 209088 X 348480; storage space: > 20645544 unneeded,6170106 used > [1] MatAssemblyEnd_SeqAIJ(): Number of mallocs during MatSetValues() is 0 > [1] MatAssemblyEnd_SeqAIJ(): Maximum nonzeros in any row is 31 > > Which seems to be in contradiction with the above error. So I am thinking > it is a different problem. > 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. Matt > 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 : > http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatMPIBAIJSetPreallocation.html > : > > "*bs* - size of block, the blocks are ALWAYS square. One can use > MatSetBlockSizes > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetBlockSizes.html#MatSetBlockSizes>() > 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 MatCreateVecs > <http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateVecs.html#MatCreateVecs>() > " > > So I cannot call MatSetBlockSizes(A,3,5,ierr) > > 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). > > Best > > Timothee > > > > 2015-12-23 6:05 GMT+09:00 Barry Smith <[email protected]>: > >> >> 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. >> >> Barry >> >> > On Dec 22, 2015, at 12:02 AM, Timothée Nicolas < >> [email protected]> wrote: >> > >> > Hi, >> > >> > 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. >> > >> > 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 >> > >> > PetscErrorCode :: ierr >> > DM :: da >> > Mat :: A >> > PetscInt :: xm >> > ISLocalToGlobalMapping :: ltog >> > PetscInt :: dims(1), starts(1) >> > PetscInt :: dim >> > MatStencil :: row(4,1), col(4,1) >> > PetscScalar :: v >> > >> > call >> DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,iten,ione,ione,PETSC_NULL_INTEGER,da,ierr) >> > call >> DMDAGetCorners(da,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, >> & >> > & xm,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr) >> > >> > call MatCreate(PETSC_COMM_WORLD,A,ierr) >> > >> > call MatSetSizes(A,xm,xm,PETSC_DETERMINE,PETSC_DETERMINE,ierr) >> > >> > call DMGetLocalToGlobalMapping(da,ltog,ierr) >> > call MatSetLocalToGlobalMapping(A,ltog,ltog,ierr) >> > call >> DMDAGetGhostCorners(da,starts(1),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, >> & >> > & >> dims(1),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,ierr) >> > dim=ione >> > call MatSetStencil(A,dim,dims,starts,ione,ierr) >> > >> > row(MatStencil_i,1) = 0 >> > col(MatStencil_i,1) = 0 >> > >> > call MatSetValuesStencil(A,ione,row,ione,col,zero,INSERT_VALUES,ierr) >> > >> > But I keep getting a segmentation fault error. >> > >> > [0]PETSC ERROR: >> ------------------------------------------------------------------------ >> > [0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, >> probably memory access out of range >> > [0]PETSC ERROR: Try option -start_in_debugger or >> -on_error_attach_debugger >> > [0]PETSC ERROR: or see >> http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind >> > [0]PETSC ERROR: or try http://valgrind.org on GNU/linux and Apple Mac >> OS X to find memory corruption errors >> > [0]PETSC ERROR: configure using --with-debugging=yes, recompile, link, >> and run >> > [0]PETSC ERROR: to get more information on the crash. >> > [0]PETSC ERROR: --------------------- Error Message >> -------------------------------------------------------------- >> > [0]PETSC ERROR: Signal received >> > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >> for trouble shooting. >> > [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015 >> > [0]PETSC ERROR: ./miips on a arch-darwin-c-debug named >> iMac27Nicolas.nifs.ac.jp by timotheenicolas Tue Dec 22 14:57:46 2015 >> > [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >> --with-fc=gfortran --download-fblaslapack --download-mpich >> --with-debugging=no >> > [0]PETSC ERROR: #1 User provided function() line 0 in unknown file >> > application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0 >> > [cli_0]: aborting job: >> > application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0 >> > >> > >> =================================================================================== >> > = BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES >> > = PID 9812 RUNNING AT iMac27Nicolas.nifs.ac.jp >> > = EXIT CODE: 59 >> > = CLEANING UP REMAINING PROCESSES >> > = YOU CAN IGNORE THE BELOW CLEANUP MESSAGES >> > >> =================================================================================== >> > >> > If I replace everything with call DMCreateMatrix(da,A,ierr) it works >> fine. >> > >> > What information did I miss from the manual ? What am I doing wrong ? >> > >> > Best >> > >> > Timothee >> > >> > >> > >> > >> > >> > 2015-12-06 23:16 GMT+09:00 Timothée Nicolas <[email protected] >> >: >> > Wow, Great answer, thanks, I should probably be able to do it like >> this, I'll try my best ! >> > >> > 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. >> > >> > Thanks >> > >> > Timothée >> > >> > >> > >> > 2015-12-06 22:29 GMT+09:00 Dave May <[email protected]>: >> > >> > >> > On 6 December 2015 at 11:19, Timothée Nicolas < >> [email protected]> wrote: >> > Hi, >> > >> > 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. >> > >> > Okay - I misunderstood. >> > >> > >> > 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. >> > >> > 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 : >> > >> > >> > Do you need to explicitly assemble this rectangular operator? >> > Expressing this operator via a matrix-free representation would be >> relatively straight forward. >> > >> > 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) >> > >> > PetscInt m,n,M,N; >> > >> > >> DMDAGetInfo(dm3,NULL,&M,&N,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL); >> > nodes = M *N; >> > DMDAGetCorners(dm3,NULL,NULL,NULL,&m,&n,NULL); >> > lnodes = m * n; >> > MatCreate(PETSC_COMM_WORLD,&A); >> > MatSetSizes(A,lnodes*5,lnodes*3,nodes*5,nodes*3); >> > >> > The sparsity pattern is defined by the connectivity of the DMDA points >> AND the different block sizes. >> > 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). >> > >> > 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 >> > >> > PetscInt *nnz,*nnz_od; >> > >> > PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz); >> > PetscMalloc(sizeof(PetscInt)*m*n*5,&nnz_od); >> > >> > PetscMemzero(nnz,sizeof(PetscInt)*m*n*5); >> > PetscMemzero(nnz_od,sizeof(PetscInt)*m*n*5); >> > >> > for (j=0; j<n; j++) { >> > for (i=0; i<m; i++) { >> > PetscInt idx,d,cnti,cntg; >> > PetscBool is_ghost; >> > >> > idx = i + j * m; >> > >> > cnti = 0; >> > cntg = 0; >> > for (all_points_in_stencil) { /* User logic to define this loop */ >> > >> > is_ghost = PETSC_FALSE; >> > /*User logic goes here to test if stencil point is a ghost point */ >> > >> > /* accumulate counts for the stencil */ >> > if (is_ghost) { cntg++; } /* values living on a neighbour rank */ >> > else { cnti++; } /* values local to current rank */ >> > } >> > >> > /* assume all dofs in row and col space in block are coupled to each >> other */ >> > for (d=0; d<5; d++) { >> > nnz[5*idx+d] = 3 * cnti; >> > nnz_od[5*idx+d] = 3 * cntg; >> > } >> > >> > } >> > } >> > >> > Thanks, >> > Dave >> > >> > >> > >> > [0]PETSC ERROR: --------------------- Error Message >> -------------------------------------------------------------- >> > [0]PETSC ERROR: Nonconforming object sizes >> > [0]PETSC ERROR: Mat mat,Vec y: global dim 4900 2940 >> > [0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html >> for trouble shooting. >> > [0]PETSC ERROR: Petsc Release Version 3.6.1, Jul, 22, 2015 >> > [0]PETSC ERROR: ./miips on a arch-linux2-c-debug named Carl-9000 by >> timothee Sun Dec 6 19:17:20 2015 >> > [0]PETSC ERROR: Configure options --with-cc=gcc --with-cxx=g++ >> --with-fc=gfortran --download-fblaslapack --download-mpich >> > [0]PETSC ERROR: #1 MatMult() line 2215 in >> /home/timothee/Documents/petsc-3.6.1/src/mat/interface/matrix.c >> > >> > So I assume it means the matrix is assumed to be square. Is there a way >> to change this ? >> > >> > Best >> > >> > Timothée >> > >> > >> > >> > 2015-12-05 20:41 GMT+09:00 Dave May <[email protected]>: >> > >> > >> > On 5 December 2015 at 11:15, Timothée Nicolas < >> [email protected]> wrote: >> > Hi, >> > >> > 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 ? >> > >> > Do the DMDA's have exactly the same parallel layout and only differ by >> the number of DOFs? >> > >> > 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. >> > >> > Thanks, >> > Dave >> > >> > >> > 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. >> > >> > Thanks >> > >> > Best >> > >> > Timothée >> > >> > >> > >> > >> > >> >> > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener
