Actual source code: ex10.c

petsc-master 2017-04-20
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>

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

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

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

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

 93:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 94:                          Load system
 95:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

311:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
312:                          Solve system
313:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

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

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

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

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

441:   PetscFinalize();
442:   return ierr;
443: }