Actual source code: ex7.c

  1: static char help[] = "Block Jacobi preconditioner for solving a linear system in parallel with KSP.\n\
  2: The code indicates the\n\
  3: procedures for setting the particular block sizes and for using different\n\
  4: linear solvers on the individual blocks.\n\n";

  6: /*
  7:    Note:  This example focuses on ways to customize the block Jacobi
  8:    preconditioner. See ex1.c and ex2.c for more detailed comments on
  9:    the basic usage of KSP (including working with matrices and vectors).

 11:    Recall: The block Jacobi method is equivalent to the ASM preconditioner
 12:    with zero overlap.
 13:  */

 15: /*T
 16:    Concepts: KSP^customizing the block Jacobi preconditioner
 17:    Processors: n
 18: T*/

 20: /*
 21:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 22:   automatically includes:
 23:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 24:      petscmat.h - matrices
 25:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 26:      petscviewer.h - viewers               petscpc.h  - preconditioners
 27: */
 28: #include <petscksp.h>

 30: int main(int argc,char **args)
 31: {
 32:   Vec            x,b,u;      /* approx solution, RHS, exact solution */
 33:   Mat            A;            /* linear system matrix */
 34:   KSP            ksp;         /* KSP context */
 35:   KSP            *subksp;     /* array of local KSP contexts on this processor */
 36:   PC             pc;           /* PC context */
 37:   PC             subpc;        /* PC context for subdomain */
 38:   PetscReal      norm;         /* norm of solution error */
 39:   PetscInt       i,j,Ii,J,*blks,m = 4,n;
 40:   PetscMPIInt    rank,size;
 41:   PetscInt       its,nlocal,first,Istart,Iend;
 42:   PetscScalar    v,one = 1.0,none = -1.0;
 43:   PetscBool      isbjacobi;

 45:   PetscInitialize(&argc,&args,(char*)0,help);
 46:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 47:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 48:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 49:   n    = m+2;

 51:   /* -------------------------------------------------------------------
 52:          Compute the matrix and right-hand-side vector that define
 53:          the linear system, Ax = b.
 54:      ------------------------------------------------------------------- */

 56:   /*
 57:      Create and assemble parallel matrix
 58:   */
 59:   MatCreate(PETSC_COMM_WORLD,&A);
 60:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 61:   MatSetFromOptions(A);
 62:   MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
 63:   MatSeqAIJSetPreallocation(A,5,NULL);
 64:   MatGetOwnershipRange(A,&Istart,&Iend);
 65:   for (Ii=Istart; Ii<Iend; Ii++) {
 66:     v = -1.0; i = Ii/n; j = Ii - i*n;
 67:     if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 68:     if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 69:     if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 70:     if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 71:     v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
 72:   }
 73:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 74:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 75:   MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);

 77:   /*
 78:      Create parallel vectors
 79:   */
 80:   MatCreateVecs(A,&u,&b);
 81:   VecDuplicate(u,&x);

 83:   /*
 84:      Set exact solution; then compute right-hand-side vector.
 85:   */
 86:   VecSet(u,one);
 87:   MatMult(A,u,b);

 89:   /*
 90:      Create linear solver context
 91:   */
 92:   KSPCreate(PETSC_COMM_WORLD,&ksp);

 94:   /*
 95:      Set operators. Here the matrix that defines the linear system
 96:      also serves as the preconditioning matrix.
 97:   */
 98:   KSPSetOperators(ksp,A,A);

100:   /*
101:      Set default preconditioner for this program to be block Jacobi.
102:      This choice can be overridden at runtime with the option
103:         -pc_type <type>

105:      IMPORTANT NOTE: Since the inners solves below are constructed to use
106:      iterative methods (such as KSPGMRES) the outer Krylov method should
107:      be set to use KSPFGMRES since it is the only Krylov method (plus KSPFCG)
108:      that allows the preconditioners to be nonlinear (that is have iterative methods
109:      inside them). The reason these examples work is because the number of
110:      iterations on the inner solves is left at the default (which is 10,000)
111:      and the tolerance on the inner solves is set to be a tight value of around 10^-6.
112:   */
113:   KSPGetPC(ksp,&pc);
114:   PCSetType(pc,PCBJACOBI);

116:   /* -------------------------------------------------------------------
117:                    Define the problem decomposition
118:      ------------------------------------------------------------------- */

120:   /*
121:      Call PCBJacobiSetTotalBlocks() to set individually the size of
122:      each block in the preconditioner.  This could also be done with
123:      the runtime option
124:          -pc_bjacobi_blocks <blocks>
125:      Also, see the command PCBJacobiSetLocalBlocks() to set the
126:      local blocks.

128:       Note: The default decomposition is 1 block per processor.
129:   */
130:   PetscMalloc1(m,&blks);
131:   for (i=0; i<m; i++) blks[i] = n;
132:   PCBJacobiSetTotalBlocks(pc,m,blks);
133:   PetscFree(blks);

135:   /*
136:     Set runtime options
137:   */
138:   KSPSetFromOptions(ksp);

140:   /* -------------------------------------------------------------------
141:                Set the linear solvers for the subblocks
142:      ------------------------------------------------------------------- */

144:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145:        Basic method, should be sufficient for the needs of most users.
146:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

148:      By default, the block Jacobi method uses the same solver on each
149:      block of the problem.  To set the same solver options on all blocks,
150:      use the prefix -sub before the usual PC and KSP options, e.g.,
151:           -sub_pc_type <pc> -sub_ksp_type <ksp> -sub_ksp_rtol 1.e-4
152:   */

154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:         Advanced method, setting different solvers for various blocks.
156:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

158:      Note that each block's KSP context is completely independent of
159:      the others, and the full range of uniprocessor KSP options is
160:      available for each block. The following section of code is intended
161:      to be a simple illustration of setting different linear solvers for
162:      the individual blocks.  These choices are obviously not recommended
163:      for solving this particular problem.
164:   */
165:   PetscObjectTypeCompare((PetscObject)pc,PCBJACOBI,&isbjacobi);
166:   if (isbjacobi) {
167:     /*
168:        Call KSPSetUp() to set the block Jacobi data structures (including
169:        creation of an internal KSP context for each block).

171:        Note: KSPSetUp() MUST be called before PCBJacobiGetSubKSP().
172:     */
173:     KSPSetUp(ksp);

175:     /*
176:        Extract the array of KSP contexts for the local blocks
177:     */
178:     PCBJacobiGetSubKSP(pc,&nlocal,&first,&subksp);

180:     /*
181:        Loop over the local blocks, setting various KSP options
182:        for each block.
183:     */
184:     for (i=0; i<nlocal; i++) {
185:       KSPGetPC(subksp[i],&subpc);
186:       if (rank == 0) {
187:         if (i%2) {
188:           PCSetType(subpc,PCILU);
189:         } else {
190:           PCSetType(subpc,PCNONE);
191:           KSPSetType(subksp[i],KSPBCGS);
192:           KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
193:         }
194:       } else {
195:         PCSetType(subpc,PCJACOBI);
196:         KSPSetType(subksp[i],KSPGMRES);
197:         KSPSetTolerances(subksp[i],1.e-6,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
198:       }
199:     }
200:   }

202:   /* -------------------------------------------------------------------
203:                       Solve the linear system
204:      ------------------------------------------------------------------- */

206:   /*
207:      Solve the linear system
208:   */
209:   KSPSolve(ksp,b,x);

211:   /* -------------------------------------------------------------------
212:                       Check solution and clean up
213:      ------------------------------------------------------------------- */

215:   /*
216:      Check the error
217:   */
218:   VecAXPY(x,none,u);
219:   VecNorm(x,NORM_2,&norm);
220:   KSPGetIterationNumber(ksp,&its);
221:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);

223:   /*
224:      Free work space.  All PETSc objects should be destroyed when they
225:      are no longer needed.
226:   */
227:   KSPDestroy(&ksp);
228:   VecDestroy(&u));  PetscCall(VecDestroy(&x);
229:   VecDestroy(&b));  PetscCall(MatDestroy(&A);
230:   PetscFinalize();
231:   return 0;
232: }

234: /*TEST

236:    test:
237:       suffix: 1
238:       nsize: 2
239:       args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

241:    test:
242:       suffix: 2
243:       nsize: 2
244:       args: -ksp_view ::ascii_info_detail

246:    test:
247:       suffix: viennacl
248:       requires: viennacl
249:       args: -ksp_monitor_short -mat_type aijviennacl
250:       output_file: output/ex7_mpiaijcusparse.out

252:    test:
253:       suffix: viennacl_2
254:       nsize: 2
255:       requires: viennacl
256:       args: -ksp_monitor_short -mat_type aijviennacl
257:       output_file: output/ex7_mpiaijcusparse_2.out

259:    test:
260:       suffix: mpiaijcusparse
261:       requires: cuda
262:       args: -ksp_monitor_short -mat_type aijcusparse

264:    test:
265:       suffix: mpiaijcusparse_2
266:       nsize: 2
267:       requires: cuda
268:       args: -ksp_monitor_short -mat_type aijcusparse

270:    test:
271:       suffix: mpiaijcusparse_simple
272:       requires: cuda
273:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse -sub_ksp_type preonly -sub_pc_type ilu

275:    test:
276:       suffix: mpiaijcusparse_simple_2
277:       nsize: 2
278:       requires: cuda
279:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse -sub_ksp_type preonly -sub_pc_type ilu

281:    test:
282:       suffix: mpiaijcusparse_3
283:       requires: cuda
284:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse

286:    test:
287:       suffix: mpiaijcusparse_4
288:       nsize: 2
289:       requires: cuda
290:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse

292:    testset:
293:       args: -ksp_monitor_short -pc_type gamg -ksp_view -pc_gamg_esteig_ksp_type cg -pc_gamg_esteig_ksp_max_it 10
294:       test:
295:         suffix: gamg_cuda
296:         nsize: {{1 2}separate output}
297:         requires: cuda
298:         args: -mat_type aijcusparse
299:         # triggers cusparse MatTransposeMat operation when squaring the graph
300:         args: -pc_gamg_sym_graph 0 -pc_gamg_threshold -1 -pc_gamg_square_graph 1
301:       test:
302:         suffix: gamg_kokkos
303:         nsize: {{1 2}separate output}
304:         requires: !sycl kokkos_kernels
305:         args: -mat_type aijkokkos

307: TEST*/