Actual source code: ex10.c

petsc-3.4.4 2014-03-13
  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:   -trans  : solve transpose system instead\n\n";
 12: /*
 13:   This code can be used to test PETSc interface to other packages.\n\
 14:   Examples of command line options:       \n\
 15:    ./ex10 -f0 <datafile> -ksp_type preonly  \n\
 16:         -help -ksp_view                  \n\
 17:         -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
 18:         -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
 19:         -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
 20:    mpiexec -n <np> ./ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
 21:  \n\n";
 22: */
 23: /*T
 24:    Concepts: KSP^solving a linear system
 25:    Processors: n
 26: T*/

 28: /*
 29:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
 30:   automatically includes:
 31:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 32:      petscmat.h - matrices
 33:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 34:      petscviewer.h - viewers               petscpc.h  - preconditioners
 35: */
 36: #include <petscksp.h>
 37: #include <petsctime.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;
 52:   PetscReal      norm;
 53:   PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
 54:   PetscBool      preload=PETSC_TRUE,isSymmetric,cknorm=PETSC_FALSE,initialguessfile = PETSC_FALSE;
 55:   PetscMPIInt    rank;
 56:   char           initialguessfilename[PETSC_MAX_PATH_LEN];

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

 66:   /*
 67:      Determine files from which we read the two linear systems
 68:      (matrix and right-hand-side vector).
 69:   */
 70:   PetscOptionsGetString(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,"-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,"-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);

112:   flg  = PETSC_FALSE;
113:   PetscOptionsGetString(NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN,&flg);
114:   VecCreate(PETSC_COMM_WORLD,&b);
115:   if (flg) {   /* rhs is stored in a separate file */
116:     if (file[2][0] == '0') {
117:       PetscInt    m;
118:       PetscScalar one = 1.0;
119:       PetscInfo(0,"Using vector of ones for RHS\n");
120:       MatGetLocalSize(A,&m,NULL);
121:       VecSetSizes(b,m,PETSC_DECIDE);
122:       VecSetFromOptions(b);
123:       VecSet(b,one);
124:     } else {
125:       PetscViewerDestroy(&fd);
126:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
127:       VecSetFromOptions(b);
128:       VecLoad(b,fd);
129:     }
130:   } else {   /* rhs is stored in the same file as matrix */
131:     VecSetFromOptions(b);
132:     VecLoad(b,fd);
133:   }
134:   PetscViewerDestroy(&fd);

136:   /* Make A singular for testing zero-pivot of ilu factorization        */
137:   /* Example: ./ex10 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
138:   flg  = PETSC_FALSE;
139:   PetscOptionsGetBool(NULL, "-test_zeropivot", &flg,NULL);
140:   if (flg) {
141:     PetscInt          row,ncols;
142:     const PetscInt    *cols;
143:     const PetscScalar *vals;
144:     PetscBool         flg1=PETSC_FALSE;
145:     PetscScalar       *zeros;
146:     row  = 0;
147:     MatGetRow(A,row,&ncols,&cols,&vals);
148:     PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
149:     PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
150:     PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);
151:     if (flg1) {   /* set entire row as zero */
152:       MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
153:     } else {   /* only set (row,row) entry as zero */
154:       MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
155:     }
156:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
157:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
158:   }

160:   /* Check whether A is symmetric */
161:   flg  = PETSC_FALSE;
162:   PetscOptionsGetBool(NULL, "-check_symmetry", &flg,NULL);
163:   if (flg) {
164:     Mat Atrans;
165:     MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
166:     MatEqual(A, Atrans, &isSymmetric);
167:     if (isSymmetric) {
168:       MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);
169:     } else {
170:       PetscPrintf(PETSC_COMM_WORLD,"Warning: A is non-symmetric \n");
171:     }
172:     MatDestroy(&Atrans);
173:   }

175:   /*
176:      If the loaded matrix is larger than the vector (due to being padded
177:      to match the block size of the system), then create a new padded vector.
178:   */

180:   MatGetLocalSize(A,&m,&n);
181:   /*  if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);*/
182:   MatGetSize(A,&M,NULL);
183:   VecGetSize(b,&m);
184:   if (M != m) {   /* Create a new vector b by padding the old one */
185:     PetscInt    j,mvec,start,end,indx;
186:     Vec         tmp;
187:     PetscScalar *bold;

189:     VecCreate(PETSC_COMM_WORLD,&tmp);
190:     VecSetSizes(tmp,n,PETSC_DECIDE);
191:     VecSetFromOptions(tmp);
192:     VecGetOwnershipRange(b,&start,&end);
193:     VecGetLocalSize(b,&mvec);
194:     VecGetArray(b,&bold);
195:     for (j=0; j<mvec; j++) {
196:       indx = start+j;
197:       VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
198:     }
199:     VecRestoreArray(b,&bold);
200:     VecDestroy(&b);
201:     VecAssemblyBegin(tmp);
202:     VecAssemblyEnd(tmp);
203:     b    = tmp;
204:   }

206:   MatGetVecs(A,&x,NULL);
207:   VecDuplicate(b,&u);
208:   if (initialguessfile) {
209:     PetscViewer viewer2;
210:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,initialguessfilename,FILE_MODE_READ,&viewer2);
211:     VecLoad(x,viewer2);
212:     PetscViewerDestroy(&viewer2);
213:     initialguess = PETSC_TRUE;
214:   } else if (initialguess) {
215:     VecSet(x,1.0);
216:   } else {
217:     VecSet(x,0.0);
218:   }


221:   /* Check scaling in A */
222:   flg  = PETSC_FALSE;
223:   PetscOptionsGetBool(NULL, "-check_scaling", &flg,NULL);
224:   if (flg) {
225:     Vec       max, min;
226:     PetscInt  idx;
227:     PetscReal val;

229:     VecDuplicate(x, &max);
230:     VecDuplicate(x, &min);
231:     MatGetRowMaxAbs(A, max, NULL);
232:     MatGetRowMinAbs(A, min, NULL);
233:     {
234:       PetscViewer viewer;

236:       PetscViewerASCIIOpen(PETSC_COMM_WORLD, "max.data", &viewer);
237:       VecView(max, viewer);
238:       PetscViewerDestroy(&viewer);
239:       PetscViewerASCIIOpen(PETSC_COMM_WORLD, "min.data", &viewer);
240:       VecView(min, viewer);
241:       PetscViewerDestroy(&viewer);
242:     }
243:     VecView(max, PETSC_VIEWER_DRAW_WORLD);
244:     VecMax(max, &idx, &val);
245:     PetscPrintf(PETSC_COMM_WORLD, "Largest max row element %G at row %d\n", val, idx);
246:     VecView(min, PETSC_VIEWER_DRAW_WORLD);
247:     VecMin(min, &idx, &val);
248:     PetscPrintf(PETSC_COMM_WORLD, "Smallest min row element %G at row %d\n", val, idx);
249:     VecMin(max, &idx, &val);
250:     PetscPrintf(PETSC_COMM_WORLD, "Smallest max row element %G at row %d\n", val, idx);
251:     VecPointwiseDivide(max, max, min);
252:     VecMax(max, &idx, &val);
253:     PetscPrintf(PETSC_COMM_WORLD, "Largest row ratio %G at row %d\n", val, idx);
254:     VecView(max, PETSC_VIEWER_DRAW_WORLD);
255:     VecDestroy(&max);
256:     VecDestroy(&min);
257:   }

259:   /*  MatView(A,PETSC_VIEWER_STDOUT_WORLD); */
260:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
261:                     Setup solve for system
262:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
263:   /*
264:      Conclude profiling last stage; begin profiling next stage.
265:   */
266:   PetscPreLoadStage("KSPSetUpSolve");

268:   /*
269:      We also explicitly time this stage via PetscTime()
270:   */
271:   PetscTime(&tsetup1);

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

297:     /*
298:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
299:      enable more precise profiling of setting up the preconditioner.
300:      These calls are optional, since both will be called within
301:      KSPSolve() if they haven't been called already.
302:     */
303:     KSPSetUp(ksp);
304:     KSPSetUpOnBlocks(ksp);
305:     PetscTime(&tsetup2);
306:     tsetup = tsetup2 - tsetup1;

308:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
309:                          Solve system
310:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

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

411:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
412:       VecView(x, viewer);
413:       PetscViewerDestroy(&viewer);
414:     }

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

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

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

438:   PetscFinalize();
439:   return 0;
440: }