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

Multiple Sparse Matrix Vector products
by Nun ion 19 Feb '12

19 Feb '12

Hello i have a conceptual idea of a sparse matvec implementation where i have multiple matrices, how would i go about implementing something such as for i = ... for k = ... w_{ik} = K_i * u_k end end Where each of the K_i are sparse matrices... the K_i are various stiffness matrices whose size can range (although they are all the same size). The u_k are reused Thanks! Mark

2 1

Specifying different ksp_types for multiple linear systems
by Max Rudolph 19 Feb '12

19 Feb '12

In my thermomechanical convection code, I set up and solve two linear systems, the first for the Stokes system and the second for an energy equation. Currently, these are separate matrices and during each timestep I first create a KSP object for each system, then solve, then destroy the KSP contexts. I would like to try out the -pc_type fieldsplit for the Stokes system, but is it possible to use -pc_type fieldsplit for the Stokes system and different pc_type and ksp_type for the energy equation? Is there perhaps a way to name the KSP associated with the Stokes system and then refer to this label from the command line (e.g. -pc_type_0 fieldsplit -pc_type_1 lu) ? Thanks very much for your help. Max Rudolph

2 1

KSP_PCApply Documentation?
by Bibrak Qamar 18 Feb '12

18 Feb '12

Hello, KSP_PCApply I found the the above function used in KSPSolve_CG but couldn't find its documentation. I will appreciate any help in this direction. thanks Bibrak

2 2

Any updates in /src/mat/examples/tests/ex81.c
by Xuefei Yuan 18 Feb '12

18 Feb '12

Hello, all, This routine is to convert the matrix to the HB format, but it seems this code is not right. Any update on it? Thanks, Rebecca

1 0

Re: [petsc-users] petsc-users Digest, Vol 38, Issue 41
by coco＠dmi.unict.it 18 Feb '12

18 Feb '12

Thank you for the answer. > Date: Fri, 17 Feb 2012 10:11:51 -0600 > From: Matthew Knepley <knepley(a)gmail.com> > Subject: Re: [petsc-users] Multigrid as a preconditioner > To: PETSc users list <petsc-users(a)mcs.anl.gov> > Message-ID: > <CAMYG4GkKE6doSQ1FgHCSr1deZRYS0kpvLQF8daqv0au==tOgyw(a)mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > On Fri, Feb 17, 2012 at 6:38 AM, <coco(a)dmi.unict.it> wrote: > >> Thank you very much for the answer, but some other doubts remain. >> >> Date: Thu, 16 Feb 2012 13:39:11 -0600 >>> From: Matthew Knepley <knepley(a)gmail.com> >>> Subject: Re: [petsc-users] Multigrid as a preconditioner >>> To: PETSc users list <petsc-users(a)mcs.anl.gov> >>> Message-ID: >>> <CAMYG4Gmt_mn_**hPVKOrbC452geco3Zz26SF05pGy6t7** >>> RPAq+0Zg(a)mail.gmail.com<CAMYG4Gmt_mn_hPVKOrbC452geco3Zz26SF05pGy6t7RPAq%2B0Zg(a)mail.gmail.com> >>>> >>> Content-Type: text/plain; charset="iso-8859-1" >>> >>> On Thu, Feb 16, 2012 at 12:54 PM, <coco(a)dmi.unict.it> wrote: >>> >>> 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. >>>> >>>> >>> Your understanding of the Richardson iteration is flawed. You can consult >>> Yousef Saad's book for the standard definition and anaysis. >>> >>> >> I think my explanation was not so clear. What I would like to do is to use >> a preconditioned Richardson iteration: >> >> x^(n+1) = x^n + P^(-1) (f-A x^n) >> >> Choosing P=A, I should expect to obtain the exact solution at the first >> iteration. Then, the whole linear system is solved by the preconditioner >> method that I chose. Is it what Petsc would do? >> > > I am not sure what you mean by "Is it what Petsc would do?". PETSc does > what you tell it to do. If you want it > to solve in one iteration, tell it to use LU, -ksp_type richardson -pc_type > lu. > Indeed I would like to solve the whole linear system by a multigrid approach and not by a lu factorization. Therefore I would like to use -ksp_type richardson -pc_type mg. In this case, the preconditioned problem P^(-1) (f-A x^n) is solved exactly or it performs just a V-cycle iteration? In both cases, since I am using a one-grid multigrid (just for debugging), it should anyway provide the exact solution at the first iteration, but it is not so. > >> >>> 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. >>>> >>>> >>> Yes, give Richardson many more iterations, -ksp_max_it. >>> >>> Matt >>> >>> >> I tried, but unfortunately nothing changed. >> Another strange phenomen is that, even with the standard KSP solver (which >> provides the right solution), if I use the flag -ksp_monitor, nothing is >> displayed. Is that an indicator of some troubles? >> > > 1) You must call KSPSetFromOptions() to use command line options > > 2) Run with -ksp_monitor -ksp_view and send the results, or there is no way > we can know what is happening. > > Matt > Thank you. I plotted the residual and it decreases so much slowly. It seems like it is using the non-preconditioned Richardson iteration x^(n+1) = x^n + P^(-1) (f-A x^n) with P=I instead of P=A. Thank you. Armando > >> Thank you in advance. >> Armando >> >> >>> 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 >>>> >>>> >>>>

4 3

Re: [petsc-users] Multigrid as a preconditioner
by coco＠dmi.unict.it 17 Feb '12

17 Feb '12

> Date: Fri, 17 Feb 2012 15:38:55 -0500 > From: Jed Brown <jedbrown(a)mcs.anl.gov> > Subject: Re: [petsc-users] petsc-users Digest, Vol 38, Issue 41 > To: PETSc users list <petsc-users(a)mcs.anl.gov> > Message-ID: > <CAM9tzSn+CaFbGRdRyc_bVQ8Wa1ACRn8LTRTg1EWxg0az9=YWEQ(a)mail.gmail.com> > Content-Type: text/plain; charset="utf-8" > > On Fri, Feb 17, 2012 at 14:09, <coco(a)dmi.unict.it> wrote: > >> Indeed I would like to solve the whole linear system by a multigrid >> approach and not by a lu factorization. Therefore I would like to use >> -ksp_type richardson -pc_type mg. >> In this case, the preconditioned problem P^(-1) (f-A x^n) is solved >> exactly or it performs just a V-cycle iteration? In both cases, since I am >> using a one-grid multigrid (just for debugging), it should anyway provide >> the exact solution at the first iteration, but it is not so. >> > > -pc_type mg with one level just applies a normal smoother. I've sometimes > thought it should do a coarse-level solve instead, but I haven't messed > with it. Barry, why doesn't it do a direct solve? > This explains why the residual decreases so slowly: because it applies the smoother instead of the coarse solver. > In general -pc_type mg does one multigrid cycle (usually a V or W cycle). > If you want to use multiple iterations, you can > > -pc_type ksp -ksp_pc_type mg > > which would use the default KSP (GMRES) as an iteration, preconditioned by > multigrid. The "outer" problem will see the result of this converged > iterative solve. I've perfectly understood that. Thank you very much. Armando

1 0

Re: [petsc-users] Multigrid as a preconditioner
by coco＠dmi.unict.it 17 Feb '12

17 Feb '12

> Date: Fri, 17 Feb 2012 13:29:16 -0600 > From: Matthew Knepley <knepley(a)gmail.com> > Subject: Re: [petsc-users] petsc-users Digest, Vol 38, Issue 41 > To: PETSc users list <petsc-users(a)mcs.anl.gov> > Message-ID: > <CAMYG4G=pQPtg_yM=im0ZZv73hfcbanKvW3vuzXj3zCTE80mR7g(a)mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > On Fri, Feb 17, 2012 at 1:09 PM, <coco(a)dmi.unict.it> wrote: > >> >> Thank you for the answer. >> >> Date: Fri, 17 Feb 2012 10:11:51 -0600 >>> From: Matthew Knepley <knepley(a)gmail.com> >>> Subject: Re: [petsc-users] Multigrid as a preconditioner >>> To: PETSc users list <petsc-users(a)mcs.anl.gov> >>> Message-ID: >>> <**CAMYG4GkKE6doSQ1FgHCSr1deZRYS0**kpvLQF8daqv0au==tOgyw(a)mail.** >>> gmail.com <tOgyw(a)mail.gmail.com>> >>> Content-Type: text/plain; charset="iso-8859-1" >>> >>> On Fri, Feb 17, 2012 at 6:38 AM, <coco(a)dmi.unict.it> wrote: >>> >>> Thank you very much for the answer, but some other doubts remain. >>>> >>>> Date: Thu, 16 Feb 2012 13:39:11 -0600 >>>> >>>>> From: Matthew Knepley <knepley(a)gmail.com> >>>>> Subject: Re: [petsc-users] Multigrid as a preconditioner >>>>> To: PETSc users list <petsc-users(a)mcs.anl.gov> >>>>> Message-ID: >>>>> <CAMYG4Gmt_mn_****hPVKOrbC452geco3Zz26SF05pGy6t7**** >>>>> RPAq+0Zg(a)mail.gmail.com<CAMYG4**Gmt_mn_**hPVKOrbC452geco3Zz26SF05pGy6t7 >>>>> **RPAq%2B0Zg(a)mail.gmail.com<CAMYG4Gmt_mn_hPVKOrbC452geco3Zz26SF05pGy6t7RPAq%252B0Zg(a)mail.gmail.com> >>>>> > >>>>> >>>>>> >>>>>> Content-Type: text/plain; charset="iso-8859-1" >>>>> >>>>> On Thu, Feb 16, 2012 at 12:54 PM, <coco(a)dmi.unict.it> wrote: >>>>> >>>>> 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. >>>>>> >>>>>> >>>>>> Your understanding of the Richardson iteration is flawed. You can >>>>> consult >>>>> Yousef Saad's book for the standard definition and anaysis. >>>>> >>>>> >>>>> I think my explanation was not so clear. What I would like to do is to >>>> use >>>> a preconditioned Richardson iteration: >>>> >>>> x^(n+1) = x^n + P^(-1) (f-A x^n) >>>> >>>> Choosing P=A, I should expect to obtain the exact solution at the first >>>> iteration. Then, the whole linear system is solved by the preconditioner >>>> method that I chose. Is it what Petsc would do? >>>> >>>> >>> I am not sure what you mean by "Is it what Petsc would do?". PETSc does >>> what you tell it to do. If you want it >>> to solve in one iteration, tell it to use LU, -ksp_type richardson >>> -pc_type >>> lu. >>> >>> >> Indeed I would like to solve the whole linear system by a multigrid >> approach and not by a lu factorization. Therefore I would like to use >> -ksp_type richardson -pc_type mg. >> In this case, the preconditioned problem P^(-1) (f-A x^n) is solved >> exactly or it performs just a V-cycle iteration? In both cases, since I am >> using a one-grid multigrid (just for debugging), it should anyway provide >> the exact solution at the first iteration, but it is not so. >> > > Great, however I have no idea what you are actually using if you do not > send what I asked for in the last message. > > Matt > > This is the output of the options: -ksp_monitor -ksp_view 0 KSP Residual norm 1.393931159942e+00 [...] 10000 KSP Residual norm 1.000095856678e+00 KSP Object: type: richardson Richardson: damping factor=1 maximum iterations=10000, initial guess is zero tolerances: relative=3.75657e-16, absolute=1e-50, divergence=10000 left preconditioning using PRECONDITIONED norm type for convergence test PC Object: type: mg MG: type is MULTIPLICATIVE, levels=1 cycles=v Cycles per PCApply=1 Coarse grid solver -- level 0 presmooths=1 postsmooths=1 ----- KSP Object:(mg_levels_0_) type: gmres GMRES: restart=30, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement GMRES: happy breakdown tolerance 1e-30 maximum iterations=1, initial guess is zero tolerances: relative=3.75657e-16, absolute=1e-50, divergence=10000 left preconditioning using PRECONDITIONED norm type for convergence test PC Object:(mg_levels_0_) type: none linear system matrix = precond matrix: Matrix Object: type=shell, rows=2662, cols=2662 linear system matrix = precond matrix: Matrix Object: type=mpiaij, rows=2662, cols=2662 total: nonzeros=12160, allocated nonzeros=31324 not using I-node (on process 0) routines Thank you. Armando >> >>> >>>> 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. >>>>>> >>>>>> >>>>>> Yes, give Richardson many more iterations, -ksp_max_it. >>>>> >>>>> Matt >>>>> >>>>> >>>>> I tried, but unfortunately nothing changed. >>>> Another strange phenomen is that, even with the standard KSP solver >>>> (which >>>> provides the right solution), if I use the flag -ksp_monitor, nothing is >>>> displayed. Is that an indicator of some troubles? >>>> >>>> >>> 1) You must call KSPSetFromOptions() to use command line options >>> >>> 2) Run with -ksp_monitor -ksp_view and send the results, or there is no >>> way >>> we can know what is happening. >>> >>> Matt >>> >>> >> Thank you. I plotted the residual and it decreases so much slowly. It >> seems like it is using the non-preconditioned Richardson iteration x^(n+1) >> = x^n + P^(-1) (f-A x^n) with P=I instead of P=A. >> >> Thank you. >> Armando >> >> >>> Thank you in advance. >>>> Armando >>>> >>>> >>>> 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 >>>>>> >>>>>> >>>>>> >>>>>>

1 0

ISAllGather withoutduplicates
by Blaise Bourdin 17 Feb '12

17 Feb '12

Hi, Is there an easy way to gather all values of an IS across all processes in a communicator while removing duplicates? Basically, I want to go from [0] Number of indices in set 2 [0] 0 1 [0] 1 2 [1] Number of indices in set 2 [1] 0 2 [1] 1 3 to [0] Number of indices in set 3 [0] 0 1 [0] 1 2 [0] 2 3 [1] Number of indices in set 3 [1] 0 1 [1] 1 2 [1] 2 3 The way I do it right now is ierr = ISGetTotalIndices(csIS,&labels);CHKERRQ(ierr); ierr = ISGetSize(csIS,&num_cs_global);CHKERRQ(ierr); ierr = PetscMalloc(num_cs_global * sizeof(PetscInt),&labels2); for (i = 0; i < num_cs_global; i++) { labels2[i] = labels[i]; } ierr = PetscSortRemoveDupsInt(&num_cs_global,labels2);CHKERRQ(ierr); ierr = ISCreateGeneral(comm,num_cs_global,labels2,PETSC_COPY_VALUES,&csIS_global);CHKERRQ(ierr); ierr = ISRestoreTotalIndices(csIS,&labels);CHKERRQ(ierr); ierr = PetscFree(labels2);CHKERRQ(ierr); but there has to be a better way (or at least one that does not involve copying from const PetscInt *labels to PetscInt *labels2, and then uses again PETSC_COPY_VALUES). Blaise -- Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

2 4

Re: [petsc-users] Multigrid as a preconditioner
by coco＠dmi.unict.it 17 Feb '12

17 Feb '12

Thank you very much for the answer, but some other doubts remain. > Date: Thu, 16 Feb 2012 13:39:11 -0600 > From: Matthew Knepley <knepley(a)gmail.com> > Subject: Re: [petsc-users] Multigrid as a preconditioner > To: PETSc users list <petsc-users(a)mcs.anl.gov> > Message-ID: > <CAMYG4Gmt_mn_hPVKOrbC452geco3Zz26SF05pGy6t7RPAq+0Zg(a)mail.gmail.com> > Content-Type: text/plain; charset="iso-8859-1" > > On Thu, Feb 16, 2012 at 12:54 PM, <coco(a)dmi.unict.it> wrote: > >> 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. >> > > Your understanding of the Richardson iteration is flawed. You can consult > Yousef Saad's book for the standard definition and anaysis. > I think my explanation was not so clear. What I would like to do is to use a preconditioned Richardson iteration: x^(n+1) = x^n + P^(-1) (f-A x^n) Choosing P=A, I should expect to obtain the exact solution at the first iteration. Then, the whole linear system is solved by the preconditioner method that I chose. Is it what Petsc would do? > >> 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. >> > > Yes, give Richardson many more iterations, -ksp_max_it. > > Matt > I tried, but unfortunately nothing changed. Another strange phenomen is that, even with the standard KSP solver (which provides the right solution), if I use the flag -ksp_monitor, nothing is displayed. Is that an indicator of some troubles? Thank you in advance. Armando > >> 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

What do with singular blocks in block matrix preconditioning?
by Thomas Witkowski 17 Feb '12

17 Feb '12

I consider a 2x2 block matrix (saddle point) with the left upper block being singular due to Neumann boundary conditions. The whole block matrix is still non-singular. I worked on some ideas for block preconditioning, but there is always some problem with the singular block. All publications I know assume the block to be definite. There is also some work on highly singular blocks, but this is here not the case. Does some of you know papers about block preconditioners for some class of 2x2 saddle point problems, where the left upper block is assumed to be positive semi-definite? From a more practical point of view, I have the problem that, independently of a special kind of block preconditioner, one has always to solve (or to approximate the solution) a system with the singular block with an arbitrary right hand side. But in general the right hand side does not fulfill the compatibility condition of having zero mean. Is there a way out of this problem? Thomas

3 4