Actual source code: ex9.c

petsc-master 2017-01-20
Report Typos and Errors

  2: static char help[] = "The solution of 2 different linear systems with different linear solvers.\n\
  3: Also, this example illustrates the repeated\n\
  4: solution of linear systems, while reusing matrix, vector, and solver data\n\
  5: structures throughout the process.  Note the various stages of event logging.\n\n";

  7: /*T
  8:    Concepts: KSP^repeatedly solving linear systems;
  9:    Concepts: PetscLog^profiling multiple stages of code;
 10:    Concepts: PetscLog^user-defined event profiling;
 11:    Processors: n
 12: T*/

 14: /*
 15:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 16:   automatically includes:
 17:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 18:      petscmat.h - matrices
 19:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 20:      petscviewer.h - viewers               petscpc.h  - preconditioners
 21: */
 22:  #include <petscksp.h>

 24: /*
 25:    Declare user-defined routines
 26: */
 27: extern PetscErrorCode CheckError(Vec,Vec,Vec,PetscInt,PetscReal,PetscLogEvent);
 28: extern PetscErrorCode MyKSPMonitor(KSP,PetscInt,PetscReal,void*);

 30: int main(int argc,char **args)
 31: {
 32:   Vec            x1,b1,x2,b2; /* solution and RHS vectors for systems #1 and #2 */
 33:   Vec            u;              /* exact solution vector */
 34:   Mat            C1,C2;         /* matrices for systems #1 and #2 */
 35:   KSP            ksp1,ksp2;   /* KSP contexts for systems #1 and #2 */
 36:   PetscInt       ntimes = 3;     /* number of times to solve the linear systems */
 37:   PetscLogEvent  CHECK_ERROR;    /* event number for error checking */
 38:   PetscInt       ldim,low,high,iglobal,Istart,Iend,Istart2,Iend2;
 39:   PetscInt       Ii,J,i,j,m = 3,n = 2,its,t;
 41:   PetscBool      flg = PETSC_FALSE;
 42:   PetscScalar    v;
 43:   PetscMPIInt    rank,size;
 44: #if defined(PETSC_USE_LOG)
 45:   PetscLogStage stages[3];
 46: #endif

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

 55:   /*
 56:      Register various stages for profiling
 57:   */
 58:   PetscLogStageRegister("Prelim setup",&stages[0]);
 59:   PetscLogStageRegister("Linear System 1",&stages[1]);
 60:   PetscLogStageRegister("Linear System 2",&stages[2]);

 62:   /*
 63:      Register a user-defined event for profiling (error checking).
 64:   */
 65:   CHECK_ERROR = 0;
 66:   PetscLogEventRegister("Check Error",KSP_CLASSID,&CHECK_ERROR);

 68:   /* - - - - - - - - - - - - Stage 0: - - - - - - - - - - - - - -
 69:                         Preliminary Setup
 70:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 72:   PetscLogStagePush(stages[0]);

 74:   /*
 75:      Create data structures for first linear system.
 76:       - Create parallel matrix, specifying only its global dimensions.
 77:         When using MatCreate(), the matrix format can be specified at
 78:         runtime. Also, the parallel partitioning of the matrix is
 79:         determined by PETSc at runtime.
 80:       - Create parallel vectors.
 81:         - When using VecSetSizes(), we specify only the vector's global
 82:           dimension; the parallel partitioning is determined at runtime.
 83:         - Note: We form 1 vector from scratch and then duplicate as needed.
 84:   */
 85:   MatCreate(PETSC_COMM_WORLD,&C1);
 86:   MatSetSizes(C1,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 87:   MatSetFromOptions(C1);
 88:   MatSetUp(C1);
 89:   MatGetOwnershipRange(C1,&Istart,&Iend);
 90:   VecCreate(PETSC_COMM_WORLD,&u);
 91:   VecSetSizes(u,PETSC_DECIDE,m*n);
 92:   VecSetFromOptions(u);
 93:   VecDuplicate(u,&b1);
 94:   VecDuplicate(u,&x1);

 96:   /*
 97:      Create first linear solver context.
 98:      Set runtime options (e.g., -pc_type <type>).
 99:      Note that the first linear system uses the default option
100:      names, while the second linear systme uses a different
101:      options prefix.
102:   */
103:   KSPCreate(PETSC_COMM_WORLD,&ksp1);
104:   KSPSetFromOptions(ksp1);

106:   /*
107:      Set user-defined monitoring routine for first linear system.
108:   */
109:   PetscOptionsGetBool(NULL,NULL,"-my_ksp_monitor",&flg,NULL);
110:   if (flg) {KSPMonitorSet(ksp1,MyKSPMonitor,NULL,0);}

112:   /*
113:      Create data structures for second linear system.
114:   */
115:   MatCreate(PETSC_COMM_WORLD,&C2);
116:   MatSetSizes(C2,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
117:   MatSetFromOptions(C2);
118:   MatSetUp(C2);
119:   MatGetOwnershipRange(C2,&Istart2,&Iend2);
120:   VecDuplicate(u,&b2);
121:   VecDuplicate(u,&x2);

123:   /*
124:      Create second linear solver context
125:   */
126:   KSPCreate(PETSC_COMM_WORLD,&ksp2);

128:   /*
129:      Set different options prefix for second linear system.
130:      Set runtime options (e.g., -s2_pc_type <type>)
131:   */
132:   KSPAppendOptionsPrefix(ksp2,"s2_");
133:   KSPSetFromOptions(ksp2);

135:   /*
136:      Assemble exact solution vector in parallel.  Note that each
137:      processor needs to set only its local part of the vector.
138:   */
139:   VecGetLocalSize(u,&ldim);
140:   VecGetOwnershipRange(u,&low,&high);
141:   for (i=0; i<ldim; i++) {
142:     iglobal = i + low;
143:     v       = (PetscScalar)(i + 100*rank);
144:     VecSetValues(u,1,&iglobal,&v,ADD_VALUES);
145:   }
146:   VecAssemblyBegin(u);
147:   VecAssemblyEnd(u);

149:   /*
150:      Log the number of flops for computing vector entries
151:   */
152:   PetscLogFlops(2.0*ldim);

154:   /*
155:      End curent profiling stage
156:   */
157:   PetscLogStagePop();

159:   /* --------------------------------------------------------------
160:                         Linear solver loop:
161:       Solve 2 different linear systems several times in succession
162:      -------------------------------------------------------------- */

164:   for (t=0; t<ntimes; t++) {

166:     /* - - - - - - - - - - - - Stage 1: - - - - - - - - - - - - - -
167:                  Assemble and solve first linear system
168:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

170:     /*
171:        Begin profiling stage #1
172:     */
173:     PetscLogStagePush(stages[1]);

175:     /*
176:        Initialize all matrix entries to zero.  MatZeroEntries() retains
177:        the nonzero structure of the matrix for sparse formats.
178:     */
179:     if (t > 0) {MatZeroEntries(C1);}

181:     /*
182:        Set matrix entries in parallel.  Also, log the number of flops
183:        for computing matrix entries.
184:         - Each processor needs to insert only elements that it owns
185:           locally (but any non-local elements will be sent to the
186:           appropriate processor during matrix assembly).
187:         - Always specify global row and columns of matrix entries.
188:     */
189:     for (Ii=Istart; Ii<Iend; Ii++) {
190:       v = -1.0; i = Ii/n; j = Ii - i*n;
191:       if (i>0)   {J = Ii - n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
192:       if (i<m-1) {J = Ii + n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
193:       if (j>0)   {J = Ii - 1; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
194:       if (j<n-1) {J = Ii + 1; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
195:       v = 4.0; MatSetValues(C1,1,&Ii,1,&Ii,&v,ADD_VALUES);
196:     }
197:     for (Ii=Istart; Ii<Iend; Ii++) { /* Make matrix nonsymmetric */
198:       v = -1.0*(t+0.5); i = Ii/n;
199:       if (i>0)   {J = Ii - n; MatSetValues(C1,1,&Ii,1,&J,&v,ADD_VALUES);}
200:     }
201:     PetscLogFlops(2.0*(Iend-Istart));

203:     /*
204:        Assemble matrix, using the 2-step process:
205:          MatAssemblyBegin(), MatAssemblyEnd()
206:        Computations can be done while messages are in transition
207:        by placing code between these two statements.
208:     */
209:     MatAssemblyBegin(C1,MAT_FINAL_ASSEMBLY);
210:     MatAssemblyEnd(C1,MAT_FINAL_ASSEMBLY);

212:     /*
213:        Indicate same nonzero structure of successive linear system matrices
214:     */
215:     MatSetOption(C1,MAT_NEW_NONZERO_LOCATIONS,PETSC_TRUE);

217:     /*
218:        Compute right-hand-side vector
219:     */
220:     MatMult(C1,u,b1);

222:     /*
223:        Set operators. Here the matrix that defines the linear system
224:        also serves as the preconditioning matrix.
225:     */
226:     KSPSetOperators(ksp1,C1,C1);

228:     /*
229:        Use the previous solution of linear system #1 as the initial
230:        guess for the next solve of linear system #1.  The user MUST
231:        call KSPSetInitialGuessNonzero() in indicate use of an initial
232:        guess vector; otherwise, an initial guess of zero is used.
233:     */
234:     if (t>0) {
235:       KSPSetInitialGuessNonzero(ksp1,PETSC_TRUE);
236:     }

238:     /*
239:        Solve the first linear system.  Here we explicitly call
240:        KSPSetUp() for more detailed performance monitoring of
241:        certain preconditioners, such as ICC and ILU.  This call
242:        is optional, ase KSPSetUp() will automatically be called
243:        within KSPSolve() if it hasn't been called already.
244:     */
245:     KSPSetUp(ksp1);
246:     KSPSolve(ksp1,b1,x1);
247:     KSPGetIterationNumber(ksp1,&its);

249:     /*
250:        Check error of solution to first linear system
251:     */
252:     CheckError(u,x1,b1,its,1.e-4,CHECK_ERROR);

254:     /* - - - - - - - - - - - - Stage 2: - - - - - - - - - - - - - -
255:                  Assemble and solve second linear system
256:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

258:     /*
259:        Conclude profiling stage #1; begin profiling stage #2
260:     */
261:     PetscLogStagePop();
262:     PetscLogStagePush(stages[2]);

264:     /*
265:        Initialize all matrix entries to zero
266:     */
267:     if (t > 0) {MatZeroEntries(C2);}

269:     /*
270:        Assemble matrix in parallel. Also, log the number of flops
271:        for computing matrix entries.
272:         - To illustrate the features of parallel matrix assembly, we
273:           intentionally set the values differently from the way in
274:           which the matrix is distributed across the processors.  Each
275:           entry that is not owned locally will be sent to the appropriate
276:           processor during MatAssemblyBegin() and MatAssemblyEnd().
277:         - For best efficiency the user should strive to set as many
278:           entries locally as possible.
279:      */
280:     for (i=0; i<m; i++) {
281:       for (j=2*rank; j<2*rank+2; j++) {
282:         v = -1.0;  Ii = j + n*i;
283:         if (i>0)   {J = Ii - n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
284:         if (i<m-1) {J = Ii + n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
285:         if (j>0)   {J = Ii - 1; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
286:         if (j<n-1) {J = Ii + 1; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
287:         v = 6.0 + t*0.5; MatSetValues(C2,1,&Ii,1,&Ii,&v,ADD_VALUES);
288:       }
289:     }
290:     for (Ii=Istart2; Ii<Iend2; Ii++) { /* Make matrix nonsymmetric */
291:       v = -1.0*(t+0.5); i = Ii/n;
292:       if (i>0)   {J = Ii - n; MatSetValues(C2,1,&Ii,1,&J,&v,ADD_VALUES);}
293:     }
294:     MatAssemblyBegin(C2,MAT_FINAL_ASSEMBLY);
295:     MatAssemblyEnd(C2,MAT_FINAL_ASSEMBLY);
296:     PetscLogFlops(2.0*(Iend-Istart));

298:     /*
299:        Indicate same nonzero structure of successive linear system matrices
300:     */
301:     MatSetOption(C2,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);

303:     /*
304:        Compute right-hand-side vector
305:     */
306:     MatMult(C2,u,b2);

308:     /*
309:        Set operators. Here the matrix that defines the linear system
310:        also serves as the preconditioning matrix.  Indicate same nonzero
311:        structure of successive preconditioner matrices by setting flag
312:        SAME_NONZERO_PATTERN.
313:     */
314:     KSPSetOperators(ksp2,C2,C2);

316:     /*
317:        Solve the second linear system
318:     */
319:     KSPSetUp(ksp2);
320:     KSPSolve(ksp2,b2,x2);
321:     KSPGetIterationNumber(ksp2,&its);

323:     /*
324:        Check error of solution to second linear system
325:     */
326:     CheckError(u,x2,b2,its,1.e-4,CHECK_ERROR);

328:     /*
329:        Conclude profiling stage #2
330:     */
331:     PetscLogStagePop();
332:   }
333:   /* --------------------------------------------------------------
334:                        End of linear solver loop
335:      -------------------------------------------------------------- */

337:   /*
338:      Free work space.  All PETSc objects should be destroyed when they
339:      are no longer needed.
340:   */
341:   KSPDestroy(&ksp1); KSPDestroy(&ksp2);
342:   VecDestroy(&x1);   VecDestroy(&x2);
343:   VecDestroy(&b1);   VecDestroy(&b2);
344:   MatDestroy(&C1);   MatDestroy(&C2);
345:   VecDestroy(&u);

347:   PetscFinalize();
348:   return ierr;
349: }
350: /* ------------------------------------------------------------- */
351: /*
352:     CheckError - Checks the error of the solution.

354:     Input Parameters:
355:     u - exact solution
356:     x - approximate solution
357:     b - work vector
358:     its - number of iterations for convergence
359:     tol - tolerance
360:     CHECK_ERROR - the event number for error checking
361:                   (for use with profiling)

363:     Notes:
364:     In order to profile this section of code separately from the
365:     rest of the program, we register it as an "event" with
366:     PetscLogEventRegister() in the main program.  Then, we indicate
367:     the start and end of this event by respectively calling
368:         PetscLogEventBegin(CHECK_ERROR,u,x,b,0);
369:         PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
370:     Here, we specify the objects most closely associated with
371:     the event (the vectors u,x,b).  Such information is optional;
372:     we could instead just use 0 instead for all objects.
373: */
374: PetscErrorCode CheckError(Vec u,Vec x,Vec b,PetscInt its,PetscReal tol,PetscLogEvent CHECK_ERROR)
375: {
376:   PetscScalar    none = -1.0;
377:   PetscReal      norm;

380:   PetscLogEventBegin(CHECK_ERROR,u,x,b,0);

382:   /*
383:      Compute error of the solution, using b as a work vector.
384:   */
385:   VecCopy(x,b);
386:   VecAXPY(b,none,u);
387:   VecNorm(b,NORM_2,&norm);
388:   if (norm > tol) {
389:     PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);
390:   }
391:   PetscLogEventEnd(CHECK_ERROR,u,x,b,0);
392:   return 0;
393: }
394: /* ------------------------------------------------------------- */
395: /*
396:    MyKSPMonitor - This is a user-defined routine for monitoring
397:    the KSP iterative solvers.

399:    Input Parameters:
400:      ksp   - iterative context
401:      n     - iteration number
402:      rnorm - 2-norm (preconditioned) residual value (may be estimated)
403:      dummy - optional user-defined monitor context (unused here)
404: */
405: PetscErrorCode MyKSPMonitor(KSP ksp,PetscInt n,PetscReal rnorm,void *dummy)
406: {
407:   Vec            x;

410:   /*
411:      Build the solution vector
412:   */
413:   KSPBuildSolution(ksp,NULL,&x);

415:   /*
416:      Write the solution vector and residual norm to stdout.
417:       - PetscPrintf() handles output for multiprocessor jobs
418:         by printing from only one processor in the communicator.
419:       - The parallel viewer PETSC_VIEWER_STDOUT_WORLD handles
420:         data from multiple processors so that the output
421:         is not jumbled.
422:   */
423:   PetscPrintf(PETSC_COMM_WORLD,"iteration %D solution vector:\n",n);
424:   VecView(x,PETSC_VIEWER_STDOUT_WORLD);
425:   PetscPrintf(PETSC_COMM_WORLD,"iteration %D KSP Residual norm %14.12e \n",n,rnorm);
426:   return 0;
427: }