Actual source code: ex7.c

petsc-master 2021-01-19
Report Typos and Errors
  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*/



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

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

 48:   PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
 49:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 50:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr);
 51:   MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr);
 52:   n    = m+2;

 54:   /* -------------------------------------------------------------------
 55:          Compute the matrix and right-hand-side vector that define
 56:          the linear system, Ax = b.
 57:      ------------------------------------------------------------------- */

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

 80:   /*
 81:      Create parallel vectors
 82:   */
 83:   MatCreateVecs(A,&u,&b);
 84:   VecDuplicate(u,&x);

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

 92:   /*
 93:      Create linear solver context
 94:   */
 95:   KSPCreate(PETSC_COMM_WORLD,&ksp);

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

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

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


120:   /* -------------------------------------------------------------------
121:                    Define the problem decomposition
122:      ------------------------------------------------------------------- */

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

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

139:   /*
140:     Set runtime options
141:   */
142:   KSPSetFromOptions(ksp);


145:   /* -------------------------------------------------------------------
146:                Set the linear solvers for the subblocks
147:      ------------------------------------------------------------------- */

149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:        Basic method, should be sufficient for the needs of most users.
151:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

159:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
160:         Advanced method, setting different solvers for various blocks.
161:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

176:        Note: KSPSetUp() MUST be called before PCBJacobiGetSubKSP().
177:     */
178:     KSPSetUp(ksp);

180:     /*
181:        Extract the array of KSP contexts for the local blocks
182:     */
183:     PCBJacobiGetSubKSP(pc,&nlocal,&first,&subksp);

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

207:   /* -------------------------------------------------------------------
208:                       Solve the linear system
209:      ------------------------------------------------------------------- */

211:   /*
212:      Solve the linear system
213:   */
214:   KSPSolve(ksp,b,x);

216:   /* -------------------------------------------------------------------
217:                       Check solution and clean up
218:      ------------------------------------------------------------------- */

220:   /*
221:      Check the error
222:   */
223:   VecAXPY(x,none,u);
224:   VecNorm(x,NORM_2,&norm);
225:   KSPGetIterationNumber(ksp,&its);
226:   PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);

228:   /*
229:      Free work space.  All PETSc objects should be destroyed when they
230:      are no longer needed.
231:   */
232:   KSPDestroy(&ksp);
233:   VecDestroy(&u);  VecDestroy(&x);
234:   VecDestroy(&b);  MatDestroy(&A);
235:   PetscFinalize();
236:   return ierr;
237: }


240: /*TEST

242:    test:
243:       suffix: 1
244:       nsize: 2
245:       args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always

247:    test:
248:       suffix: 2
249:       nsize: 2
250:       args: -ksp_view

252:    test:
253:       suffix: viennacl
254:       requires: viennacl
255:       args: -ksp_monitor_short -mat_type aijviennacl
256:       output_file: output/ex7_mpiaijcusparse.out

258:    test:
259:       suffix: viennacl_2
260:       nsize: 2
261:       requires: viennacl
262:       args: -ksp_monitor_short -mat_type aijviennacl
263:       output_file: output/ex7_mpiaijcusparse_2.out

265:    test:
266:       suffix: mpiaijcusparse
267:       requires: cuda
268:       args: -ksp_monitor_short -mat_type aijcusparse

270:    test:
271:       suffix: mpiaijcusparse_2
272:       nsize: 2
273:       requires: cuda
274:       args: -ksp_monitor_short -mat_type aijcusparse

276:    test:
277:       suffix: mpiaijcusparse_simple
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_simple_2
283:       nsize: 2
284:       requires: cuda
285:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse -sub_ksp_type preonly -sub_pc_type ilu

287:    test:
288:       suffix: mpiaijcusparse_3
289:       requires: cuda
290:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse

292:    test:
293:       suffix: mpiaijcusparse_4
294:       nsize: 2
295:       requires: cuda
296:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse

298:    testset:
299:       args: -ksp_monitor_short -pc_type gamg -ksp_view
300:       test:
301:         suffix: gamg_cuda
302:         nsize: {{1 2}separate output}
303:         requires: cuda
304:         args: -mat_type aijcusparse
305:         # triggers cusparse MatTransposeMat operation when squaring the graph
306:         args: -pc_gamg_sym_graph 0 -pc_gamg_threshold -1 -pc_gamg_square_graph 1
307:       test:
308:         suffix: gamg_kokkos
309:         nsize: {{1 2}separate output}
310:         requires: kokkos_kernels
311:         args: -mat_type aijkokkos

313: TEST*/