Actual source code: ex7.c

petsc-master 2019-12-10
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);
 51:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 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);

 79:   /*
 80:      Create parallel vectors
 81:   */
 82:   VecCreate(PETSC_COMM_WORLD,&u);
 83:   VecSetSizes(u,PETSC_DECIDE,m*n);
 84:   VecSetFromOptions(u);
 85:   VecDuplicate(u,&b);
 86:   VecDuplicate(b,&x);

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

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

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

105:   /*
106:      Set default preconditioner for this program to be block Jacobi.
107:      This choice can be overridden at runtime with the option
108:         -pc_type <type>
109:         
110:      IMPORTANT NOTE: Since the inners solves below are constructed to use 
111:      iterative methods (such as KSPGMRES) the outer Krylov method should
112:      be set to use KSPFGMRES since it is the only Krylov method (plus KSPFCG)
113:      that allows the preconditioners to be nonlinear (that is have iterative methods
114:      inside them). The reason these examples work is because the number of 
115:      iterations on the inner solves is left at the default (which is 10,000)
116:      and the tolerance on the inner solves is set to be a tight value of around 10^-6.
117:   */
118:   KSPGetPC(ksp,&pc);
119:   PCSetType(pc,PCBJACOBI);


122:   /* -------------------------------------------------------------------
123:                    Define the problem decomposition
124:      ------------------------------------------------------------------- */

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

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


142:   /* -------------------------------------------------------------------
143:                Set the linear solvers for the subblocks
144:      ------------------------------------------------------------------- */

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

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

156:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157:         Advanced method, setting different solvers for various blocks.
158:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

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

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

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

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

204:   /* -------------------------------------------------------------------
205:                       Solve the linear system
206:      ------------------------------------------------------------------- */

208:   /*
209:     Set runtime options
210:   */
211:   KSPSetFromOptions(ksp);

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

218:   /* -------------------------------------------------------------------
219:                       Check solution and clean up
220:      ------------------------------------------------------------------- */

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

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


242: /*TEST

244:    test:
245:       nsize: 2
246:       args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always> ex7_1.tmp 2>&1

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

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

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

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

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

277:    test:
278:       suffix: mpiaijcusparse_simple
279:       requires: cuda
280:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse -vec_type cuda -sub_ksp_type preonly -sub_pc_type ilu

282:    test:
283:       suffix: mpiaijcusparse_simple_2
284:       nsize: 2
285:       requires: cuda
286:       args: -ksp_monitor_short -mat_type aijcusparse -sub_pc_factor_mat_solver_type cusparse -vec_type cuda -sub_ksp_type preonly -sub_pc_type ilu

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

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

299: TEST*/