Dear PETSc developers, I use Libmesh with PETSc for my FEM simulation. The program works when using MATSEQCUSP. However, When I test MATMPICUSP with 2 GPU cards, I met the following errors: " [1]PETSC ERROR: --------------------- Error Message ------------------------------------ [1]PETSC ERROR: Argument out of range! [1]PETSC ERROR: New nonzero at (1,5) caused a malloc! [1]PETSC ERROR: ------------------------------------------------------------------------ [1]PETSC ERROR: Petsc Development HG revision: b4086d236fb35071ea565635e24ea99d70deeaac HG Date: Sat Jan 21 21:50:39 2012 -0600 [1]PETSC ERROR: See docs/changes/index.html for recent updates. [1]PETSC ERROR: See docs/faq.html for hints about trouble shooting. [1]PETSC ERROR: See docs/index.html for manual pages. [1]PETSC ERROR: ------------------------------------------------------------------------ [1]PETSC ERROR: /work/01820/ylu/libmesh_svn01232012_petscme/examples/myproj_sp1/myproj-opt on a linux named c300-202.ls4.tacc.utexas.edu by ylu Fri Feb 10 13:52:39 2012 [1]PETSC ERROR: Libraries linked from /work/01820/ylu/petsc-dev01222012_real/linux/lib [1]PETSC ERROR: Configure run at Thu Feb 9 10:42:41 2012 [1]PETSC ERROR: Configure options --with-clanguage=C++ --with-debugging=1 --with-shared-libraries=1 --with-mpi-compilers=1 --with-mpi-dir=/opt/apps/intel11_1/mvapich2/1.6 --with-blas-lapack-dir=/opt/apps/intel/11.1/mkl/lib/em64t/ --with-cuda=1 --with-cusp=1 --with-thrust=1 --with-cuda-dir=/opt/apps/cuda/4.0/cuda/ --with-cusp-dir=/work/01820/ylu --with-valgrind-dir=/opt/apps/valgrind/3.6.0/ [1]PETSC ERROR: ------------------------------------------------------------------------ [1]PETSC ERROR: MatSetValues_MPIAIJ() line 505 in src/mat/impls/aij/mpi/mpiaij.c [1]PETSC ERROR: MatSetValues() line 1119 in src/mat/interface/matrix.c [1]PETSC ERROR: add_matrix() line 468 in "unknowndirectory/"src/numerics/petsc_matrix.C application called MPI_Abort(comm=0x84000004, 63) - process 1 [0]PETSC ERROR: --------------------- Error Message ------------------------------------ [0]PETSC ERROR: Argument out of range! [0]PETSC ERROR: New nonzero at (2,54) caused a malloc! [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: Petsc Development HG revision: b4086d236fb35071ea565635e24ea99d70deeaac HG Date: Sat Jan 21 21:50:39 2012 -0600 [0]PETSC ERROR: See docs/changes/index.html for recent updates. [0]PETSC ERROR: See docs/faq.html for hints about trouble shooting. [0]PETSC ERROR: See docs/index.html for manual pages. [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: /work/01820/ylu/libmesh_svn01232012_petscme/examples/myproj_sp1/myproj-opt on a linux named c300-202.ls4.tacc.utexas.edu by ylu Fri Feb 10 13:52:39 2012 [0]PETSC ERROR: Libraries linked from /work/01820/ylu/petsc-dev01222012_real/linux/lib [0]PETSC ERROR: Configure run at Thu Feb 9 10:42:41 2012 [0]PETSC ERROR: Configure options --with-clanguage=C++ --with-debugging=1 --with-shared-libraries=1 --with-mpi-compilers=1 --with-mpi-dir=/opt/apps/intel11_1/mvapich2/1.6 --with-blas-lapack-dir=/opt/apps/intel/11.1/mkl/lib/em64t/ --with-cuda=1 --with-cusp=1 --with-thrust=1 --with-cuda-dir=/opt/apps/cuda/4.0/cuda/ --with-cusp-dir=/work/01820/ylu --with-valgrind-dir=/opt/apps/valgrind/3.6.0/ [0]PETSC ERROR: ------------------------------------------------------------------------ [0]PETSC ERROR: MatSetValues_MPIAIJ() line 505 in src/mat/impls/aij/mpi/mpiaij.c [0]PETSC ERROR: MatSetValues() line 1119 in src/mat/interface/matrix.c [0]PETSC ERROR: add_matrix() line 468 in "unknowndirectory/"src/numerics/petsc_matrix.C application called MPI_Abort(comm=0xC4000000, 63) - process 0 " It seem MatSetValues sets some values out of range. In parallel FEM simulation, I partition the mesh into 2 submeshes. Test is performed with 2GPU cards/2 CPU cores. Each CPU is responsible for one submesh for matrix assembly (I think MatSetValues should set the values within the submesh, otherwise it should generate errors like the above). However, I don't understand why there are errors. Is anything wrong? Thank you very much, Yujie
On Fri, Feb 10, 2012 at 14:24, recrusader <[email protected]> wrote:
I use Libmesh with PETSc for my FEM simulation. The program works when using MATSEQCUSP. However, When I test MATMPICUSP with 2 GPU cards, I met the following errors: " [1]PETSC ERROR: --------------------- Error Message ------------------------------------ [1]PETSC ERROR: Argument out of range! [1]PETSC ERROR: New nonzero at (1,5) caused a malloc!
Can you try with plain MATMPIAIJ? I suspect you are not preallocating correctly. *Preallocation routines now automatically set MAT_NEW_NONZERO_ALLOCATION_ERR, if you intentionally preallocate less than necessary then use MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) to disable the error generation.* http://www.mcs.anl.gov/petsc/documentation/changes/dev.html
Dear Jed, MATMPIAIJ works ;(. Best, Yujie On Fri, Feb 10, 2012 at 2:26 PM, Jed Brown <[email protected]> wrote:
On Fri, Feb 10, 2012 at 14:24, recrusader <[email protected]> wrote:
I use Libmesh with PETSc for my FEM simulation. The program works when using MATSEQCUSP. However, When I test MATMPICUSP with 2 GPU cards, I met the following errors: " [1]PETSC ERROR: --------------------- Error Message ------------------------------------ [1]PETSC ERROR: Argument out of range! [1]PETSC ERROR: New nonzero at (1,5) caused a malloc!
Can you try with plain MATMPIAIJ? I suspect you are not preallocating correctly.
*Preallocation routines now automatically set MAT_NEW_NONZERO_ALLOCATION_ERR, if you intentionally preallocate less than necessary then use MatSetOption(mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) to disable the error generation.*
On Fri, Feb 10, 2012 at 14:35, recrusader <[email protected]> wrote:
MATMPIAIJ works ;(.
Can you send a reproducible test case? (Instructions for how to run a Libmesh example that would cause this problem is okay, producing it with a PETSc example would be better.)
Since MATMPIAIJ works (I didn't change anything. Just set the vec and mat types to mpicusp and mpiaijcusp for GPU), I think the problem is likely from MatSetValues_MPIAIJ(). Which PETSc examples can test this function? I can try to use it to test. Thanks a lot, Jed. Best, Yujie On Fri, Feb 10, 2012 at 2:44 PM, Jed Brown <[email protected]> wrote:
On Fri, Feb 10, 2012 at 14:35, recrusader <[email protected]> wrote:
MATMPIAIJ works ;(.
Can you send a reproducible test case? (Instructions for how to run a Libmesh example that would cause this problem is okay, producing it with a PETSc example would be better.)
On Fri, Feb 10, 2012 at 14:49, recrusader <[email protected]> wrote:
Since MATMPIAIJ works (I didn't change anything. Just set the vec and mat types to mpicusp and mpiaijcusp for GPU), I think the problem is likely from MatSetValues_MPIAIJ().
Which PETSc examples can test this function?
Try src/ksp/ksp/examples/tutorials/ex43.c and src/snes/examples/tutorials/ex48.c
Dear Jed, The first example works. However, the example uses MatSetValuesStencil() not MatSetValues(). Are they same? Thanks a lot, Yujie On Fri, Feb 10, 2012 at 2:51 PM, Jed Brown <[email protected]> wrote:
On Fri, Feb 10, 2012 at 14:49, recrusader <[email protected]> wrote:
Since MATMPIAIJ works (I didn't change anything. Just set the vec and mat types to mpicusp and mpiaijcusp for GPU), I think the problem is likely from MatSetValues_MPIAIJ().
Which PETSc examples can test this function?
Try src/ksp/ksp/examples/tutorials/ex43.c and src/snes/examples/tutorials/ex48.c
On Fri, Feb 10, 2012 at 3:12 PM, recrusader <[email protected]> wrote:
Dear Jed,
The first example works. However, the example uses MatSetValuesStencil() not MatSetValues(). Are they same?
MatSetValuesStencil() calls MatSetValues(). I suspect that your MPIAIJ matrix does not have the option set to throw an error when inserting a new nonzero, and your MPICUSP matrix does. Matt
Thanks a lot, Yujie
On Fri, Feb 10, 2012 at 2:51 PM, Jed Brown <[email protected]> wrote:
On Fri, Feb 10, 2012 at 14:49, recrusader <[email protected]> wrote:
Since MATMPIAIJ works (I didn't change anything. Just set the vec and mat types to mpicusp and mpiaijcusp for GPU), I think the problem is likely from MatSetValues_MPIAIJ().
Which PETSc examples can test this function?
Try src/ksp/ksp/examples/tutorials/ex43.c and src/snes/examples/tutorials/ex48.c
-- 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
Dear Matt, I added the print codes in libmesh after creating the matrix as follows: " ierr = MatCreateMPIAIJ (libMesh::COMM_WORLD, m_local, n_local, m_global, n_global, PETSC_NULL, (int*) &n_nz[0], PETSC_NULL, (int*) &n_oz[0], &_mat); CHKERRABORT(libMesh::COMM_WORLD,ierr); MatSetOption(_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); //by Yujie std::cout<<"MatSetOption"<<std::endl;" I run the same codes in CPU and GPU modes (the same parameters except that GPU uses '-vec_type mpicusp -mat_type mpiaijcusp'). I can find "MatSetOption" output from both the modes. Does that mean that the codes set the options for both the modes? Thank you very much. Best, Yujie On Fri, Feb 10, 2012 at 3:14 PM, Matthew Knepley <[email protected]> wrote:
On Fri, Feb 10, 2012 at 3:12 PM, recrusader <[email protected]> wrote:
Dear Jed,
The first example works. However, the example uses MatSetValuesStencil() not MatSetValues(). Are they same?
MatSetValuesStencil() calls MatSetValues(). I suspect that your MPIAIJ matrix does not have the option set to throw an error when inserting a new nonzero, and your MPICUSP matrix does.
Matt
Thanks a lot, Yujie
On Fri, Feb 10, 2012 at 2:51 PM, Jed Brown <[email protected]> wrote:
On Fri, Feb 10, 2012 at 14:49, recrusader <[email protected]> wrote:
Since MATMPIAIJ works (I didn't change anything. Just set the vec and mat types to mpicusp and mpiaijcusp for GPU), I think the problem is likely from MatSetValues_MPIAIJ().
Which PETSc examples can test this function?
Try src/ksp/ksp/examples/tutorials/ex43.c and src/snes/examples/tutorials/ex48.c
-- 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
On Fri, Feb 10, 2012 at 11:10 PM, recrusader <[email protected]> wrote:
Dear Matt,
I added the print codes in libmesh after creating the matrix as follows: " ierr = MatCreateMPIAIJ (libMesh::COMM_WORLD, m_local, n_local, m_global, n_global, PETSC_NULL, (int*) &n_nz[0], PETSC_NULL, (int*) &n_oz[0], &_mat); CHKERRABORT(libMesh::COMM_WORLD,ierr);
MatSetOption(_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); //by Yujie std::cout<<"MatSetOption"<<std::endl;"
I run the same codes in CPU and GPU modes (the same parameters except that GPU uses '-vec_type mpicusp -mat_type mpiaijcusp'). I can find "MatSetOption" output from both the modes. Does that mean that the codes set the options for both the modes? Thank you very much.
Yes, so you should have no problem with allocation errors. Partial reports like this help no one. It would be somewhat helpful to include a stack trace, to verify that after this change you see an error. If might actually enable us to find your error if you sent a small test code which failed. Thanks, Matt
Best, Yujie
On Fri, Feb 10, 2012 at 3:14 PM, Matthew Knepley <[email protected]>wrote:
On Fri, Feb 10, 2012 at 3:12 PM, recrusader <[email protected]> wrote:
Dear Jed,
The first example works. However, the example uses MatSetValuesStencil() not MatSetValues(). Are they same?
MatSetValuesStencil() calls MatSetValues(). I suspect that your MPIAIJ matrix does not have the option set to throw an error when inserting a new nonzero, and your MPICUSP matrix does.
Matt
Thanks a lot, Yujie
On Fri, Feb 10, 2012 at 2:51 PM, Jed Brown <[email protected]> wrote:
On Fri, Feb 10, 2012 at 14:49, recrusader <[email protected]> wrote:
Since MATMPIAIJ works (I didn't change anything. Just set the vec and mat types to mpicusp and mpiaijcusp for GPU), I think the problem is likely from MatSetValues_MPIAIJ().
Which PETSc examples can test this function?
Try src/ksp/ksp/examples/tutorials/ex43.c and src/snes/examples/tutorials/ex48.c
-- 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
-- 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
Dear Matt, in http://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/ksp/examples/tutorials/ex43.c... Why are two matrices set to MATAIJ? if I set MATAIJCUSP, is them changed? 1468: DMCreateMatrix<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix> (da_Stokes,MATAIJ<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ> ,&A); 1469: DMCreateMatrix <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix>(da_Stokes,MATAIJ <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ>,&B); Thanks a lot. Best, Yujie On Fri, Feb 10, 2012 at 11:30 PM, Matthew Knepley <[email protected]> wrote:
On Fri, Feb 10, 2012 at 11:10 PM, recrusader <[email protected]> wrote:
Dear Matt,
I added the print codes in libmesh after creating the matrix as follows: " ierr = MatCreateMPIAIJ (libMesh::COMM_WORLD, m_local, n_local, m_global, n_global, PETSC_NULL, (int*) &n_nz[0], PETSC_NULL, (int*) &n_oz[0], &_mat); CHKERRABORT(libMesh::COMM_WORLD,ierr);
MatSetOption(_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); //by Yujie std::cout<<"MatSetOption"<<std::endl;"
I run the same codes in CPU and GPU modes (the same parameters except that GPU uses '-vec_type mpicusp -mat_type mpiaijcusp'). I can find "MatSetOption" output from both the modes. Does that mean that the codes set the options for both the modes? Thank you very much.
Yes, so you should have no problem with allocation errors. Partial reports like this help no one. It would be somewhat helpful to include a stack trace, to verify that after this change you see an error. If might actually enable us to find your error if you sent a small test code which failed.
Thanks,
Matt
Best, Yujie
On Fri, Feb 10, 2012 at 3:14 PM, Matthew Knepley <[email protected]>wrote:
On Fri, Feb 10, 2012 at 3:12 PM, recrusader <[email protected]>wrote:
Dear Jed,
The first example works. However, the example uses MatSetValuesStencil() not MatSetValues(). Are they same?
MatSetValuesStencil() calls MatSetValues(). I suspect that your MPIAIJ matrix does not have the option set to throw an error when inserting a new nonzero, and your MPICUSP matrix does.
Matt
Thanks a lot, Yujie
On Fri, Feb 10, 2012 at 2:51 PM, Jed Brown <[email protected]>wrote:
On Fri, Feb 10, 2012 at 14:49, recrusader <[email protected]>wrote:
Since MATMPIAIJ works (I didn't change anything. Just set the vec and mat types to mpicusp and mpiaijcusp for GPU), I think the problem is likely from MatSetValues_MPIAIJ().
Which PETSc examples can test this function?
Try src/ksp/ksp/examples/tutorials/ex43.c and src/snes/examples/tutorials/ex48.c
-- 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
-- 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
On Sat, Feb 11, 2012 at 00:00, recrusader <[email protected]> wrote:
in http://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/ksp/examples/tutorials/ex43.c...
Why are two matrices set to MATAIJ? if I set MATAIJCUSP, is them changed?
1468: DMCreateMatrix<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix> (da_Stokes,MATAIJ<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ> ,&A);
1469: DMCreateMatrix <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix>(da_Stokes,MATAIJ <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ>,&B);
You have to follow the code down to where those matrices are used. 1478: AssembleA_Stokes(A,da_Stokes,da_prop,properties); 1479: AssembleA_PCStokes(B,da_Stokes,da_prop,properties); This uses a "block preconditioner", the B matrix is intentionally different from the A matrix in order to handle the saddle point.
Dear Jed, I am sorry that I cannot generate the errors using PETSc by itself. However, it easy to generate the errors with libmesh. There is an example from libmesh, that is http://libmesh.sourceforge.net/ex4.php I can run it using the command "ibrun -n 2 -o 0 ./introduction_ex4-dbg -d 3 -pc_type none -ksp_type gmres -vec_type mpicusp -mat_type mpiaijcusp -ksp_view -ksp_monitor -log_summary -malloc_debug -cuda_show_devices" (You should replace ibrun -n 2 -o 0 using mpiexec -n 2) I get the same errors. However, when using "ibrun -n 2 -o 0 ./introduction_ex4-dbg -d 3 -pc_type none -ksp_type gmres -ksp_view -ksp_monitor -log_summary". I can run it successfully. Thank you very much. Best, Yujie On Sat, Feb 11, 2012 at 4:43 AM, Jed Brown <[email protected]> wrote:
On Sat, Feb 11, 2012 at 00:00, recrusader <[email protected]> wrote:
in http://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/ksp/examples/tutorials/ex43.c...
Why are two matrices set to MATAIJ? if I set MATAIJCUSP, is them changed?
1468: DMCreateMatrix<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix> (da_Stokes,MATAIJ<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ> ,&A);
1469: DMCreateMatrix <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix>(da_Stokes,MATAIJ <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ>,&B);
You have to follow the code down to where those matrices are used.
1478: AssembleA_Stokes(A,da_Stokes,da_prop,properties); 1479: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
This uses a "block preconditioner", the B matrix is intentionally different from the A matrix in order to handle the saddle point.
Dear Jed, When I removed 'if (NONEW == -2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at (%D,%D) caused a malloc",ROW,COL); \' from the following function in src/mat/impls/aij/seq/aij.h. It works in GPU mode. Do you have any comments? Thanks a lot. #define MatSeqXAIJReallocateAIJ(Amat,AM,BS2,NROW,ROW,COL,RMAX,AA,AI,AJ,RP,AP,AIMAX,NONEW,datatype) \ if (NROW >= RMAX) {\ Mat_SeqAIJ *Ain = (Mat_SeqAIJ*)Amat->data;\ /* there is no extra room in row, therefore enlarge */ \ PetscInt CHUNKSIZE = 15,new_nz = AI[AM] + CHUNKSIZE,len,*new_i=0,*new_j=0; \ datatype *new_a; \ \ /*if (NONEW == -2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at (%D,%D) caused a malloc",ROW,COL); \ */ \ /* malloc new storage space */ \ ierr = PetscMalloc3(BS2*new_nz,datatype,&new_a,new_nz,PetscInt,&new_j,AM+1,PetscInt,&new_i);CHKERRQ(ierr);\ \ /* copy over old data into new slots */ \ for (ii=0; ii<ROW+1; ii++) {new_i[ii] = AI[ii];} \ for (ii=ROW+1; ii<AM+1; ii++) {new_i[ii] = AI[ii]+CHUNKSIZE;} \ ierr = PetscMemcpy(new_j,AJ,(AI[ROW]+NROW)*sizeof(PetscInt));CHKERRQ(ierr); \ len = (new_nz - CHUNKSIZE - AI[ROW] - NROW); \ ierr = PetscMemcpy(new_j+AI[ROW]+NROW+CHUNKSIZE,AJ+AI[ROW]+NROW,len*sizeof(PetscInt));CHKERRQ(ierr); \ ierr = PetscMemcpy(new_a,AA,BS2*(AI[ROW]+NROW)*sizeof(datatype));CHKERRQ(ierr); \ ierr = PetscMemzero(new_a+BS2*(AI[ROW]+NROW),BS2*CHUNKSIZE*sizeof(datatype));CHKERRQ(ierr);\ ierr = PetscMemcpy(new_a+BS2*(AI[ROW]+NROW+CHUNKSIZE),AA+BS2*(AI[ROW]+NROW),BS2*len*sizeof(datatype));CHKERRQ(ierr); \ /* free up old matrix storage */ \ ierr = MatSeqXAIJFreeAIJ(A,&Ain->a,&Ain->j,&Ain->i);CHKERRQ(ierr);\ AA = new_a; \ Ain->a = (MatScalar*) new_a; \ AI = Ain->i = new_i; AJ = Ain->j = new_j; \ Ain->singlemalloc = PETSC_TRUE; \ \ RP = AJ + AI[ROW]; AP = AA + BS2*AI[ROW]; \ RMAX = AIMAX[ROW] = AIMAX[ROW] + CHUNKSIZE; \ Ain->maxnz += BS2*CHUNKSIZE; \ Ain->reallocs++; \ } \ Best, Yujie On 2/11/12, Jed Brown <[email protected]> wrote:
On Sat, Feb 11, 2012 at 00:00, recrusader <[email protected]> wrote:
in http://www.mcs.anl.gov/petsc/petsc-dev/src/ksp/ksp/examples/tutorials/ex43.c...
Why are two matrices set to MATAIJ? if I set MATAIJCUSP, is them changed?
1468: DMCreateMatrix<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix> (da_Stokes,MATAIJ<http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ> ,&A);
1469: DMCreateMatrix <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/DM/DMCreateMatrix.html#DMCreateMatrix>(da_Stokes,MATAIJ <http://www.mcs.anl.gov/petsc/petsc-dev/docs/manualpages/Mat/MATAIJ.html#MATAIJ>,&B);
You have to follow the code down to where those matrices are used.
1478: AssembleA_Stokes(A,da_Stokes,da_prop,properties); 1479: AssembleA_PCStokes(B,da_Stokes,da_prop,properties);
This uses a "block preconditioner", the B matrix is intentionally different from the A matrix in order to handle the saddle point.
On Sat, Feb 11, 2012 at 10:52, recrusader <[email protected]> wrote:
When I removed 'if (NONEW == -2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at (%D,%D) caused a malloc",ROW,COL); \' from the following function in src/mat/impls/aij/seq/aij.h. It works in GPU mode. Do you have any comments? Thanks a lot.
If you want that effect, you can MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); or -mat_new_nonzero_allocation_err 0. The more serious problem is that preallocation information seems to be getting lost. If you don't fix that, assembly will be horrendously slow. Are you sure the Mat type is being set *before* the call to MatMPIAIJSetPreallocation()?
As you suggest, I added the print codes in libmesh after creating the matrix and do MatSetOption as follows: " ierr = MatCreateMPIAIJ (libMesh::COMM_WORLD, m_local, n_local, m_global, n_global, PETSC_NULL, (int*) &n_nz[0], PETSC_NULL, (int*) &n_oz[0], &_mat); CHKERRABORT(libMesh::COMM_WORLD,ierr); MatSetOption(_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); //by Yujie std::cout<<"MatSetOption"<<std::endl;" I run the same codes in CPU and GPU modes (the same parameters except that GPU uses '-vec_type mpicusp -mat_type mpiaijcusp'). I can find "MatSetOption" output from both the modes. Does that mean that the codes set the options for both the modes? Thank you very much. Best, Yujie On Sat, Feb 11, 2012 at 10:56 AM, Jed Brown <[email protected]> wrote:
On Sat, Feb 11, 2012 at 10:52, recrusader <[email protected]> wrote:
When I removed 'if (NONEW == -2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"New nonzero at (%D,%D) caused a malloc",ROW,COL); \' from the following function in src/mat/impls/aij/seq/aij.h. It works in GPU mode. Do you have any comments? Thanks a lot.
If you want that effect, you can
MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
or -mat_new_nonzero_allocation_err 0.
The more serious problem is that preallocation information seems to be getting lost. If you don't fix that, assembly will be horrendously slow. Are you sure the Mat type is being set *before* the call to MatMPIAIJSetPreallocation()?
On Sat, Feb 11, 2012 at 10:58, recrusader <[email protected]> wrote:
" ierr = MatCreateMPIAIJ (libMesh::COMM_WORLD,
m_local, n_local, m_global, n_global, PETSC_NULL, (int*) &n_nz[0], PETSC_NULL, (int*) &n_oz[0], &_mat); CHKERRABORT(libMesh::COMM_WORLD,ierr);
MatSetOption(_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); //by Yujie std::cout<<"MatSetOption"<<std::endl;"
I run the same codes in CPU and GPU modes (the same parameters except that GPU uses '-vec_type mpicusp -mat_type mpiaijcusp'). I can find "MatSetOption" output from both the modes. Does that mean that the codes set the options for both the modes?
Libmesh is calling MatSetFromOptions() after MatCreateMPIAIJ() which means the preallocation information will be lost if the type is changed. The code should be written differently ierr = MatCreate(comm,A);CHKERRQ(ierr); ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); ierr = MatSetType(*A,MATAIJ);CHKERRQ(ierr); ierr = MatSetOptionsPrefix(*A,optional_prefix);CHKERRQ(ierr); ierr = MatSetFromOptions(*A);CHKERRQ(ierr); ierr = MatMPIAIJSetPreallocation(*A,d_nz,d_nnz,o_nz,o_nnz);CHKERRQ(ierr); ierr = MatSeqAIJSetPreallocation(*A,d_nz,d_nnz);CHKERRQ(ierr); I can talk to the libmesh developers about making this change.
That's great. Thank you very much, Jed :). Best, Yujie On Sat, Feb 11, 2012 at 11:06 AM, Jed Brown <[email protected]> wrote:
On Sat, Feb 11, 2012 at 10:58, recrusader <[email protected]> wrote:
" ierr = MatCreateMPIAIJ (libMesh::COMM_WORLD,
m_local, n_local, m_global, n_global, PETSC_NULL, (int*) &n_nz[0], PETSC_NULL, (int*) &n_oz[0], &_mat); CHKERRABORT(libMesh::COMM_WORLD,ierr);
MatSetOption(_mat,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE); //by Yujie std::cout<<"MatSetOption"<<std::endl;"
I run the same codes in CPU and GPU modes (the same parameters except that GPU uses '-vec_type mpicusp -mat_type mpiaijcusp'). I can find "MatSetOption" output from both the modes. Does that mean that the codes set the options for both the modes?
Libmesh is calling MatSetFromOptions() after MatCreateMPIAIJ() which means the preallocation information will be lost if the type is changed. The code should be written differently
ierr = MatCreate(comm,A);CHKERRQ(ierr); ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr); ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); ierr = MatSetType(*A,MATAIJ);CHKERRQ(ierr); ierr = MatSetOptionsPrefix(*A,optional_prefix);CHKERRQ(ierr); ierr = MatSetFromOptions(*A);CHKERRQ(ierr); ierr = MatMPIAIJSetPreallocation(*A,d_nz,d_nnz,o_nz,o_nnz);CHKERRQ(ierr); ierr = MatSeqAIJSetPreallocation(*A,d_nz,d_nnz);CHKERRQ(ierr);
I can talk to the libmesh developers about making this change.
participants (3)
-
Jed Brown -
Matthew Knepley -
recrusader