Actual source code: ex10.c

petsc-3.7.3 2016-07-24
Report Typos and Errors
  2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
  3: This version first preloads and solves a small system, then loads \n\
  4: another (larger) system and solves it as well.  This example illustrates\n\
  5: preloading of instructions with the smaller system so that more accurate\n\
  6: performance monitoring can be done with the larger one (that actually\n\
  7: is the system of interest).  See the 'Performance Hints' chapter of the\n\
  8: users manual for a discussion of preloading.  Input parameters include\n\
  9:   -f0 <input_file> : first file to load (small system)\n\
 10:   -f1 <input_file> : second file to load (larger system)\n\n\
 11:   -nearnulldim <0> : number of vectors in the near-null space immediately following matrix\n\n\
 12:   -trans  : solve transpose system instead\n\n";
 13: /*
 14:   This code can be used to test PETSc interface to other packages.\n\
 15:   Examples of command line options:       \n\
 16:    ./ex10 -f0 <datafile> -ksp_type preonly  \n\
 17:         -help -ksp_view                  \n\
 18:         -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
 19:         -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
 20:         -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
 21:    mpiexec -n <np> ./ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
 22:  \n\n";
 23: */
 24: /*T
 25:    Concepts: KSP^solving a linear system
 26:    Processors: n
 27: T*/

 29: /*
 30:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 31:   automatically includes:
 32:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 33:      petscmat.h - matrices
 34:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 35:      petscviewer.h - viewers               petscpc.h  - preconditioners
 36: */
 37: #include <petscksp.h>

 41: int main(int argc,char **args)
 42: {
 43:   KSP            ksp;             /* linear solver context */
 44:   Mat            A;               /* matrix */
 45:   Vec            x,b,u;           /* approx solution, RHS, exact solution */
 46:   PetscViewer    fd;              /* viewer */
 47:   char           file[4][PETSC_MAX_PATH_LEN];     /* input file name */
 48:   PetscBool      table     =PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE;
 49:   PetscBool      outputSoln=PETSC_FALSE;
 51:   PetscInt       its,num_numfac,m,n,M,nearnulldim = 0;
 52:   PetscReal      norm;
 53:   PetscBool      preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE;
 54:   PetscMPIInt    rank;
 55:   char           initialguessfilename[PETSC_MAX_PATH_LEN];

 57:   PetscInitialize(&argc,&args,(char*)0,help);
 58:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 59:   PetscOptionsGetBool(NULL,NULL,"-table",&table,NULL);
 60:   PetscOptionsGetBool(NULL,NULL,"-trans",&trans,NULL);
 61:   PetscOptionsGetBool(NULL,NULL,"-initialguess",&initialguess,NULL);
 62:   PetscOptionsGetBool(NULL,NULL,"-output_solution",&outputSoln,NULL);
 63:   PetscOptionsGetString(NULL,NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);
 64:   PetscOptionsGetInt(NULL,NULL,"-nearnulldim",&nearnulldim,NULL);

 66:   /*
 67:      Determine files from which we read the two linear systems
 68:      (matrix and right-hand-side vector).
 69:   */
 70:   PetscOptionsGetString(NULL,NULL,"-f",file[0],PETSC_MAX_PATH_LEN,&flg);
 71:   if (flg) {
 72:     PetscStrcpy(file[1],file[0]);
 73:     preload = PETSC_FALSE;
 74:   } else {
 75:     PetscOptionsGetString(NULL,NULL,"-f0",file[0],PETSC_MAX_PATH_LEN,&flg);
 76:     if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 or -f option");
 77:     PetscOptionsGetString(NULL,NULL,"-f1",file[1],PETSC_MAX_PATH_LEN,&flg);
 78:     if (!flg) preload = PETSC_FALSE;   /* don't bother with second system */
 79:   }

 81:   /* -----------------------------------------------------------
 82:                   Beginning of linear solver loop
 83:      ----------------------------------------------------------- */
 84:   /*
 85:      Loop through the linear solve 2 times.
 86:       - The intention here is to preload and solve a small system;
 87:         then load another (larger) system and solve it as well.
 88:         This process preloads the instructions with the smaller
 89:         system so that more accurate performance monitoring (via
 90:         -log_summary) can be done with the larger one (that actually
 91:         is the system of interest).
 92:   */
 93:   PetscPreLoadBegin(preload,"Load system");

 95:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 96:                          Load system
 97:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 99:   /*
100:      Open binary file.  Note that we use FILE_MODE_READ to indicate
101:      reading from this file.
102:   */
103:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);

105:   /*
106:      Load the matrix and vector; then destroy the viewer.
107:   */
108:   MatCreate(PETSC_COMM_WORLD,&A);
109:   MatSetFromOptions(A);
110:   MatLoad(A,fd);
111:   if (nearnulldim) {
112:     MatNullSpace nullsp;
113:     Vec          *nullvecs;
114:     PetscInt     i;
115:     PetscMalloc1(nearnulldim,&nullvecs);
116:     for (i=0; i<nearnulldim; i++) {
117:       VecCreate(PETSC_COMM_WORLD,&nullvecs[i]);
118:       VecLoad(nullvecs[i],fd);
119:     }
120:     MatNullSpaceCreate(PETSC_COMM_WORLD,PETSC_FALSE,nearnulldim,nullvecs,&nullsp);
121:     MatSetNearNullSpace(A,nullsp);
122:     for (i=0; i<nearnulldim; i++) {VecDestroy(&nullvecs[i]);}
123:     PetscFree(nullvecs);
124:     MatNullSpaceDestroy(&nullsp);
125:   }

127:   flg  = PETSC_FALSE;
128:   PetscOptionsGetString(NULL,NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN,&flg);
129:   VecCreate(PETSC_COMM_WORLD,&b);
130:   if (flg) {   /* rhs is stored in a separate file */
131:     if (file[2][0] == '0' || file[2][0] == 0) {
132:       PetscInt    m;
133:       PetscScalar one = 1.0;
134:       PetscInfo(0,"Using vector of ones for RHS\n");
135:       MatGetLocalSize(A,&m,NULL);
136:       VecSetSizes(b,m,PETSC_DECIDE);
137:       VecSetFromOptions(b);
138:       VecSet(b,one);
139:     } else {
140:       PetscViewerDestroy(&fd);
141:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
142:       VecSetFromOptions(b);
143:       VecLoad(b,fd);
144:     }
145:   } else {   /* rhs is stored in the same file as matrix */
146:     VecSetFromOptions(b);
147:     VecLoad(b,fd);
148:   }
149:   PetscViewerDestroy(&fd);

151:   /* Make A singular for testing zero-pivot of ilu factorization        */
152:   /* Example: ./ex10 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_type <shift_type> */
153:   flg  = PETSC_FALSE;
154:   PetscOptionsGetBool(NULL,NULL, "-test_zeropivot", &flg,NULL);
155:   if (flg) {
156:     PetscInt          row=0;
157:     PetscBool         flg1=PETSC_FALSE;
158:     PetscScalar       zero=0.0;

160:     PetscOptionsGetBool(NULL,NULL, "-set_row_zero", &flg1,NULL);
161:     if (flg1) {   /* set a row as zeros */
162:       MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,PETSC_TRUE);
163:       MatZeroRows(A,1,&row,0.0,NULL,NULL);
164:     } else if (!rank) {
165:       /* set (row,row) entry as zero */
166:       MatSetValues(A,1,&row,1,&row,&zero,INSERT_VALUES);
167:     }
168:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
169:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
170:   }

172:   /* Check whether A is symmetric */
173:   flg  = PETSC_FALSE;
174:   PetscOptionsGetBool(NULL,NULL, "-check_symmetry", &flg,NULL);
175:   if (flg) {
176:     Mat Atrans;
177:     MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
178:     MatEqual(A, Atrans, &isSymmetric);
179:     if (isSymmetric) {
180:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
181:     } else {
182:       PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");
183:     }
184:     MatDestroy(&Atrans);
185:   }

187:   /*
188:      If the loaded matrix is larger than the vector (due to being padded
189:      to match the block size of the system), then create a new padded vector.
190:   */

192:   MatGetLocalSize(A,&m,&n);
193:   /*  if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);*/
194:   MatGetSize(A,&M,NULL);
195:   VecGetSize(b,&m);
196:   if (M != m) {   /* Create a new vector b by padding the old one */
197:     PetscInt    j,mvec,start,end,indx;
198:     Vec         tmp;
199:     PetscScalar *bold;

201:     VecCreate(PETSC_COMM_WORLD,&tmp);
202:     VecSetSizes(tmp,n,PETSC_DECIDE);
203:     VecSetFromOptions(tmp);
204:     VecGetOwnershipRange(b,&start,&end);
205:     VecGetLocalSize(b,&mvec);
206:     VecGetArray(b,&bold);
207:     for (j=0; j<mvec; j++) {
208:       indx = start+j;
209:       VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
210:     }
211:     VecRestoreArray(b,&bold);
212:     VecDestroy(&b);
213:     VecAssemblyBegin(tmp);
214:     VecAssemblyEnd(tmp);
215:     b    = tmp;
216:   }

218:   MatCreateVecs(A,&x,NULL);
219:   VecDuplicate(b,&u);
220:   if (initialguessfile) {
221:     PetscViewer viewer2;
222:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer2);
223:     VecLoad(x,viewer2);
224:     PetscViewerDestroy(&viewer2);
225:     initialguess = PETSC_TRUE;
226:   } else if (initialguess) {
227:     VecSet(x,1.0);
228:   } else {
229:     VecSet(x,0.0);
230:   }


233:   /* Check scaling in A */
234:   flg  = PETSC_FALSE;
235:   PetscOptionsGetBool(NULL,NULL, "-check_scaling", &flg,NULL);
236:   if (flg) {
237:     Vec       max, min;
238:     PetscInt  idx;
239:     PetscReal val;

241:     VecDuplicate(x, &max);
242:     VecDuplicate(x, &min);
243:     MatGetRowMaxAbs(A, max, NULL);
244:     MatGetRowMinAbs(A, min, NULL);
245:     {
246:       PetscViewer viewer;

248:       PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);
249:       VecView(max, viewer);
250:       PetscViewerDestroy(&viewer);
251:       PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);
252:       VecView(min, viewer);
253:       PetscViewerDestroy(&viewer);
254:     }
255:     VecView(max, PETSC_VIEWER_DRAW_WORLD);
256:     VecMax(max, &idx, &val);
257:     PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %g at row %D\n", (double)val, idx);
258:     VecView(min, PETSC_VIEWER_DRAW_WORLD);
259:     VecMin(min, &idx, &val);
260:     PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %g at row %D\n", (double)val, idx);
261:     VecMin(max, &idx, &val);
262:     PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %g at row %D\n", (double)val, idx);
263:     VecPointwiseDivide(max, max, min);
264:     VecMax(max, &idx, &val);
265:     PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %g at row %D\n", (double)val, idx);
266:     VecView(max, PETSC_VIEWER_DRAW_WORLD);
267:     VecDestroy(&max);
268:     VecDestroy(&min);
269:   }

271:   /*  MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
272:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
273:                     Setup solve for system
274:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
275:   /*
276:      Conclude profiling last stage; begin profiling next stage.
277:   */
278:   PetscPreLoadStage("KSPSetUpSolve");

280:   /*
281:      Create linear solver; set operators; set runtime options.
282:   */
283:   KSPCreate(PETSC_COMM_WORLD,&ksp);
284:   KSPSetInitialGuessNonzero(ksp,initialguess);
285:   num_numfac = 1;
286:   PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL);
287:   while (num_numfac--) {
288:     PetscBool lsqr;
289:     char      str[32];
290:     PetscOptionsGetString(NULL,NULL,"-ksp_type",str,32,&lsqr);
291:     if (lsqr) {
292:       PetscStrcmp("lsqr",str,&lsqr);
293:     }
294:     if (lsqr) {
295:       Mat BtB;
296:       MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,4,&BtB);
297:       KSPSetOperators(ksp,A,BtB);
298:       MatDestroy(&BtB);
299:     } else {
300:       KSPSetOperators(ksp,A,A);
301:     }
302:     KSPSetFromOptions(ksp);

304:     /*
305:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
306:      enable more precise profiling of setting up the preconditioner.
307:      These calls are optional, since both will be called within
308:      KSPSolve() if they haven't been called already.
309:     */
310:     KSPSetUp(ksp);
311:     KSPSetUpOnBlocks(ksp);

313:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
314:                          Solve system
315:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

317:     /*
318:      Solve linear system; 
319:     */
320:     if (trans) {
321:       KSPSolveTranspose(ksp,b,x);
322:       KSPGetIterationNumber(ksp,&its);
323:     } else {
324:       PetscInt num_rhs=1;
325:       PetscOptionsGetInt(NULL,NULL,"-num_rhs",&num_rhs,NULL);
326:       cknorm = PETSC_FALSE;
327:       PetscOptionsGetBool(NULL,NULL,"-cknorm",&cknorm,NULL);
328:       while (num_rhs--) {
329:         if (num_rhs == 1) VecSet(x,0.0);
330:         KSPSolve(ksp,b,x);
331:       }
332:       KSPGetIterationNumber(ksp,&its);
333:       if (cknorm) {     /* Check error for each rhs */
334:         if (trans) {
335:           MatMultTranspose(A,x,u);
336:         } else {
337:           MatMult(A,x,u);
338:         }
339:         VecAXPY(u,-1.0,b);
340:         VecNorm(u,NORM_2,&norm);
341:         PetscPrintf(PETSC_COMM_WORLD,"  Number of iterations = %3D\n",its);
342:         if (!PetscIsNanScalar(norm)) {
343:           if (norm < 1.e-12) {
344:             PetscPrintf(PETSC_COMM_WORLD,"  Residual norm < 1.e-12\n");
345:           } else {
346:             PetscPrintf(PETSC_COMM_WORLD,"  Residual norm %g\n",(double)norm);
347:           }
348:         }
349:       }
350:     }   /* while (num_rhs--) */

352:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
353:           Check error, print output, free data structures.
354:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

356:     /*
357:        Check error
358:     */
359:     if (trans) {
360:       MatMultTranspose(A,x,u);
361:     } else {
362:       MatMult(A,x,u);
363:     }
364:     VecAXPY(u,-1.0,b);
365:     VecNorm(u,NORM_2,&norm);
366:     /*
367:      Write output (optinally using table for solver details).
368:       - PetscPrintf() handles output for multiprocessor jobs
369:         by printing from only one processor in the communicator.
370:       - KSPView() prints information about the linear solver.
371:     */
372:     if (table) {
373:       char        *matrixname,kspinfo[120];
374:       PetscViewer viewer;

376:       /*
377:        Open a string viewer; then write info to it.
378:       */
379:       PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
380:       KSPView(ksp,viewer);
381:       PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
382:       PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n",matrixname,its,norm,kspinfo);

384:       /*
385:         Destroy the viewer
386:       */
387:       PetscViewerDestroy(&viewer);
388:     } else {
389:       PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
390:       if (!PetscIsNanScalar(norm)) {
391:         if (norm < 1.e-12 && !PetscIsNanScalar((PetscScalar)norm)) {
392:           PetscPrintf(PETSC_COMM_WORLD,"  Residual norm < 1.e-12\n");
393:         } else {
394:           PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
395:         }
396:       }
397:     }
398:     PetscOptionsGetString(NULL,NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
399:     if (flg) {
400:       PetscViewer viewer;
401:       Vec         xstar;
402:       PetscReal   norm;

404:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
405:       VecCreate(PETSC_COMM_WORLD,&xstar);
406:       VecLoad(xstar,viewer);
407:       VecAXPY(xstar, -1.0, x);
408:       VecNorm(xstar, NORM_2, &norm);
409:       PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)norm);
410:       VecDestroy(&xstar);
411:       PetscViewerDestroy(&viewer);
412:     }
413:     if (outputSoln) {
414:       PetscViewer viewer;

416:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
417:       VecView(x, viewer);
418:       PetscViewerDestroy(&viewer);
419:     }

421:     flg  = PETSC_FALSE;
422:     PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);
423:     if (flg) {
424:       KSPConvergedReason reason;
425:       KSPGetConvergedReason(ksp,&reason);
426:       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
427:     }

429:   }   /* while (num_numfac--) */

431:   /*
432:      Free work space.  All PETSc objects should be destroyed when they
433:      are no longer needed.
434:   */
435:   MatDestroy(&A); VecDestroy(&b);
436:   VecDestroy(&u); VecDestroy(&x);
437:   KSPDestroy(&ksp);
438:   PetscPreLoadEnd();
439:   /* -----------------------------------------------------------
440:                       End of linear solver loop
441:      ----------------------------------------------------------- */

443:   PetscFinalize();
444:   return 0;
445: }