February 2012 - petsc-users - mailman3.cels.anl.gov

ParMETIS_V3_PartGeomKway
by Mohammad Mirzadeh 16 Feb '12

16 Feb '12

Hi guys, I'm wondering if there is any implementation for ParMETIS_V3_PartGeomKway()? All I can find is the implementation for ParMETIS_V3_PartKway and I'm wondering if including the vertex positions could help me get a better partitioning? Also I have a general question. Is minimizing number of edge cuts essentially the same as minimizing communication? I understand that they are related, but communication is proportional to the number of ghost points which is not exactly equal (its actually less than) to the number of edge cuts. So then is it possible that this could actually result in larger number of ghost points, at least for some of processors? Thanks, Mohammad

4 9

Multigrid as a preconditioner
by coco＠dmi.unict.it 16 Feb '12

16 Feb '12

Dear list, I would like to parallelize a multigrid code by Petsc. I do not want to use the DMMG infrastructure, since it will be replaced in the next PETSc release. Therefore I preferred to use the multigrid as a preconditioner. In practice, I use the Richardson iteration, choosing the same matrix of the linear system as a preconditioner, so that I think the Richardson iteration should converge in only one iteration, and effectively it is like solving the whole linear system by the multigrid. As a first test, I tried to use a one-grid multigrid (then not a truly multigrid). I just set the coarse solver (by a standard KSPSolve), and it should be enough, because the multigrid starts already from the coarsest grid and then it does not need the smoother and the transfer operators. Unfortunately, the iteration scheme (which should be the Richardson scheme) converges with a reason=4 (KSP_CONVERGED_ITS) to a wrong solution. On the other hand, if I solve the whole problem with the standard KSPSolve (then withouth setting the multigrid as a preconditioner ...), it converges to the right solution with a reason=2. I thought that the two methods should be the exactly the same method, and I do not understand why they provide different convergence results. Here is the relevant code: // Set the matrix of the linear system Mat Mcc; ierr=MatCreate(PETSC_COMM_WORLD,&Mcc); CHKERRQ(ierr); ierr=MatSetType(Mcc, MATMPIAIJ); CHKERRQ(ierr); ierr=MatSetSizes(Mcc,PETSC_DECIDE,PETSC_DECIDE,1000,1000); CHKERRQ(ierr); ierr=setMatrix(Mcc); //It is a routine that set the values of the matrix Mcc // Set the ksp solver with the multigrid as a preconditioner KSP ksp, KspSolver; ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr); ierr = KSPSetType(ksp,KSPRICHARDSON); ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr); ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); ierr = PCMGSetLevels(pc,1,&PETSC_COMM_WORLD);CHKERRQ(ierr); ierr = PCMGSetType(pc,PC_MG_MULTIPLICATIVE);CHKERRQ(ierr); ierr = PCMGGetCoarseSolve(pc,&kspCoarseSolve);CHKERRQ(ierr); ierr = KSPSetOperators(kspCoarseSolve,Mcc,Mcc,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); ierr = KSPSetTolerances(kspCoarseSolve,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); ierr = KSPSetOperators(ksp,Mcc,Mcc,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); ierr = KSPSetTolerances(ksp,1.e-12,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); ierr = KSPSolve(ksp,RHS,U);CHKERRQ(ierr); // Solve with the standard KSPSolve KSP ksp1; ierr = KSPCreate(PETSC_COMM_WORLD,&ksp1);CHKERRQ(ierr); ierr = KSPSetOperators(ksp1,Mcc,Mcc,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); ierr = KSPSetTolerances(ksp1,1.e-12/(2*nn123),PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);CHKERRQ(ierr); ierr = KSPSolve(ksp1,RHS,U1);CHKERRQ(ierr); At the end, the Vector U and U1 are different. Thank you. Best regards, Armando

2 1

preconditioned eigensolver in slepc
by Verena Kuhlemann 16 Feb '12

16 Feb '12

Hi, I am trying to use LOBPCG via Slepc and I am a little confused on how to set the preconditioner and operator matrices. I have a shell matrix A and a matrix P that should be used for the preconditioning. Here is the relevant part of my code: EPSCreate(PETSC_COMM_WORLD,&eps); EPSSetOperators(eps, A, PETSC_NULL); EPSSetProblemType(eps, EPS_HEP); EPSSetType(eps, EPSBLOPEX); EPSGetST(eps, &st); STPrecondSetMatForPC(st, P); STGetKSP(st,&ksp); KSPGetPC(ksp, &pc); KSPSetOperators(ksp,A, P, DIFFERENT_NONZERO_PATTERN); //do I have to set this?? If I don't the program complains that a matrix is not set PCSetType(pc, PCJACOBI); EPSSetUp(eps); Is that correct or do I need to do something else? Thanks, Verena

2 3

Ghost point communication for Red-Black Gauss-Siedel
by Sylvain Barbot 15 Feb '12

15 Feb '12

Hi All, I'm implementing an elasticity multi-grid solver with Petsc with matrix shells. I am using the Red-Black Gauss-Siedel smoother. In the middle of the Red-Black passes, I need to update the ghost points between processors. To do that effectively, I'd like to update only the ghost points between the different processors, instead of the whole array. The goal is to do the most possible operations in place in local array lx, and having to update global vector x only once. The working matrix multiply smoothing operation is shown below. I'm using Petsc 3.1. I read about http://www.mcs.anl.gov/petsc/petsc-3.1/docs/manualpages/Vec/VecGhostUpdateB…, but I'm not entirely clear whether this does what I want or not. SUBROUTINE matrixsmooth(A,b,omega,flag,shift,its,level,x,ierr) USE types USE heap USE petscsys USE petscda USE petscis USE petscvec USE petscmat USE petscpc USE petscksp IMPLICIT NONE Mat, INTENT(IN) :: A Vec, INTENT(IN) :: b Vec, INTENT(INOUT) :: x PetscReal, INTENT(IN) :: omega,shift MatSORType, INTENT(IN) :: flag PetscInt, INTENT(IN) :: its, level ! iterations, multi-grid level PetscErrorCode, INTENT(OUT) :: ierr Vec :: lx,lb PetscInt :: istart,iend PetscScalar, POINTER :: px(:) ! pointer to solution array PetscScalar, POINTER :: pb(:) ! pointer to body-force array INTEGER :: rank,isize,i,k INTEGER :: sw,off,p INTEGER :: i2,i3,i2i,i3i, & i000,i0p0,i00p, & i0m0,i00m TYPE(DALOCALINFOF90) :: info CALL MPI_COMM_RANK(PETSC_COMM_WORLD,rank,ierr) CALL MPI_COMM_SIZE(PETSC_COMM_WORLD,isize,ierr) CALL VecGetOwnershipRange(x,istart,iend,ierr) ! allocate memory for local vector with ghost points CALL DACreateLocalVector(c%daul(1+level),lx,ierr) CALL VecDuplicate(lx,lb,ierr) ! retrieve forcing term b with ghost points CALL DAGlobalToLocalBegin(c%daul(1+level),b,INSERT_VALUES,lb,ierr) CALL DAGlobalToLocalEnd(c%daul(1+level),b,INSERT_VALUES,lb,ierr) ! obtain a pointer to local data CALL VecGetArrayF90(lx,px,ierr); CALL VecGetArrayF90(lb,pb,ierr); ! geometry info about local vector and padding CALL DAGetLocalInfoF90(c%daul(1+level),info,ierr); ! retrieve stencil width (ghost-point padding) CALL DAGetInfo(c%daul(1+level),PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER, & sw,PETSC_NULL_INTEGER, & PETSC_NULL_INTEGER,ierr); ! offset due to node topology off=MOD(info%xm+info%ym,2) ! fast smoothing (relaxation) its=number of iteration DO k=1,its ! Red-Black Gauss-Siedel scheme DO p=0,1 ! retrieve initial guess x with ghost points CALL DAGlobalToLocalBegin(c%daul(1+level),x,INSERT_VALUES,lx,ierr) CALL DAGlobalToLocalEnd(c%daul(1+level),x,INSERT_VALUES,lx,ierr) ! smoothing (relaxation) DO i=istart+p+off,iend-1,info%dof*2 i3=(i-istart)/(info%xm*info%dof) i2=(i-istart-i3*info%xm*info%dof)/info%dof i3=i3+i3i ! i3 in ( 0 .. sx3-1 ) i2=i2+i2i ! i2 in ( 0 .. sx2-1 ) i000=1+((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+0)*info%dof i0p0=1+((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i+1)*info%dof i00p=1+((sw+i3-i3i+1)*(info%xm+2*sw)+sw+i2-i2i+0)*info%dof i0m0=1+((sw+i3-i3i+0)*(info%xm+2*sw)+sw+i2-i2i-1)*info%dof i00m=1+((sw+i3-i3i-1)*(info%xm+2*sw)+sw+i2-i2i+0)*info%dof px(i000+IU1)=px(i000+IU1)+ & (pb(i000+IU1)-((+(px(i0p0+IU1)-px(i000+IU1) )-(px(i000+IU1)-px(i0m0+IU1) )) & +(+(px(i00p+IU1)-px(i000+IU1) )-(px(i000+IU1)-px(i00m+IU1) )))) / (-4._8) END DO ! publish new values of x for its global vector CALL DALocalToGlobal(c%daul(1+level),lx,INSERT_VALUES,x,ierr) END DO END DO ! dispose of the local vectors CALL VecRestoreArrayF90(lx,px,ierr) CALL VecRestoreArrayF90(lb,pb,ierr) CALL DARestoreLocalVector(c%daul(1+level),lx,ierr) CALL DARestoreLocalVector(c%daul(1+level),lb,ierr) END SUBROUTINE matrixsmooth Could this be changed to something like: ... INITIALIZE... CALL DAGlobalToLocalBegin(c%daul(1+c%level),x,INSERT_VALUES,lx,ierr) CALL DAGlobalToLocalEnd(c%daul(1+c%level),x,INSERT_VALUES,lx,ierr) ! fast smoothing (relaxation) its=number of iteration DO k=1,its ! Red-Black Gauss-Siedel scheme DO p=0,1 ! smoothing (relaxation) DO i=istart+p+off,iend-1,info%dof*2 ...STENCIL OPERATION... END DO ...UPDATE GHOST POINTS... END DO END DO ! publish new values of x for its global vector CALL DALocalToGlobal(c%daul(1+level),lx,INSERT_VALUES,x,ierr) ...CLEAN MEMORY... Any recommendation? Best wishes, Sylvain Barbot

3 6

PETSc should not be used to attempt to provide a “parallel linear solver” ??
by Bibrak Qamar 15 Feb '12

15 Feb '12

Hello all, I am new to PETSc so here comes a newbie question. On the first page of the manual it says that "• PETSc should not be used to attempt to provide a “parallel linear solver” in an otherwise sequential code. Certainly all parts of a previously sequential code need not be parallelized but the matrix generation portion must be parallelized to expect any kind of reasonable performance. Do not expect to generate your matrix sequentially and then “use PETSc” to solve the linear system in parallel. " Please help me understand what is meant by it. What if I read the matrix from a file (sequential) and distribute it to other processes and then call any linear solver, how is it going to affect the performance of the code? because reading or generating the matrix is a one time cost for one timestep, or may be that CAUTION is for multiple timesteps? Thnaks Bibrak

2 1

MatMult doesn't give the correct result
by Longmin RAN 15 Feb '12

15 Feb '12

Dear all, I tried to use MatMult on a 2*2 matrix but could not get the right answer. With MatView I saw the Matrix was correctly defined as : row 0: (0, 0.707107) row 1: (1, 0.707107) The vector to be multiplied is 0 -0.0141421 But MatMult gave the result vector as 0 -0.01 Has anyone a clue to the problem? Thank you for you help in advance. Cheers, Longmin

1 1

valgrind
by Mohammad Mirzadeh 15 Feb '12

15 Feb '12

Hi, This might just not be important but I'm still curious about it. When I run valgrind on this: int main (int argc, char **argv){ PetscInitialize(&argc, &argv,(char *)0,help); PetscFinalize(); return 0; } I get, ==7262== HEAP SUMMARY: ==7262== in use at exit: 300 bytes in 11 blocks ==7262== total heap usage: 203 allocs, 192 frees, 643,759 bytes allocated ==7262== ==7262== LEAK SUMMARY: ==7262== definitely lost: 60 bytes in 1 blocks ==7262== indirectly lost: 240 bytes in 10 blocks ==7262== possibly lost: 0 bytes in 0 blocks ==7262== still reachable: 0 bytes in 0 blocks ==7262== suppressed: 0 bytes in 0 blocks ==7262== Rerun with --leak-check=full to see details of leaked memory Is this normal? Thanks, Mohammad

2 1

Null space of discrete Laplace with periodic boundary conditions
by Thomas Witkowski 14 Feb '12

14 Feb '12

I discretize the Laplace operator (using finite element) on the unit square equipped with periodic boundary conditions on all four edges. Is it correct that the null space is still constant? I wounder, because when I run the same code on a sphere (so a 2D surface embedded in 3D), the resulting matrix is non-singular. I thought, that both cases should be somehow equal with respect to the null space? Thomas

3 4

-mat_partitioning_view
by Mohammad Mirzadeh 13 Feb '12

13 Feb '12

Hi, When using -mat_partitioning_view flag, PETSc reports number of edge cuts. Is this the number of edge cuts in the graph for the new partitioning (i.e. using new global numbering)? Also, does the partitioning favors minimizing the edge cuts over having equal number of points per processor? If so, how can this be altered? Thanks, Mohammad

2 2

MatSetValues() for MATMPICUSP
by recrusader 11 Feb '12

11 Feb '12

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

3 17