Actual source code: ex10.c

petsc-3.5.1 2014-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,"-table",&table,NULL);
 60:   PetscOptionsGetBool(NULL,"-trans",&trans,NULL);
 61:   PetscOptionsGetBool(NULL,"-initialguess",&initialguess,NULL);
 62:   PetscOptionsGetBool(NULL,"-output_solution",&outputSoln,NULL);
 63:   PetscOptionsGetString(NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);
 64:   PetscOptionsGetInt(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,"-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);
111:   if (nearnulldim) {
112:     MatNullSpace nullsp;
113:     Vec *nullvecs;
114:     PetscInt i;
115:     PetscMalloc(nearnulldim*sizeof(nullvecs[0]),&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,"-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_nonzero */
153:   flg  = PETSC_FALSE;
154:   PetscOptionsGetBool(NULL, "-test_zeropivot", &flg,NULL);
155:   if (flg) {
156:     PetscInt          row,ncols;
157:     const PetscInt    *cols;
158:     const PetscScalar *vals;
159:     PetscBool         flg1=PETSC_FALSE;
160:     PetscScalar       *zeros;
161:     row  = 0;
162:     MatGetRow(A,row,&ncols,&cols,&vals);
163:     PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
164:     PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
165:     PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);
166:     if (flg1) {   /* set entire row as zero */
167:       MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
168:     } else {   /* only set (row,row) entry as zero */
169:       MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
170:     }
171:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
172:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
173:   }

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

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

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

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

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


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

244:     VecDuplicate(x, &max);
245:     VecDuplicate(x, &min);
246:     MatGetRowMaxAbs(A, max, NULL);
247:     MatGetRowMinAbs(A, min, NULL);
248:     {
249:       PetscViewer viewer;

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

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

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

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

316:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
317:                          Solve system
318:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

385:       /*
386:         Destroy the viewer
387:       */
388:       PetscViewerDestroy(&viewer);
389:     } else {
390:       PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
391:       if (norm < 1.e-12) {
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:     PetscOptionsGetString(NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
398:     if (flg) {
399:       PetscViewer viewer;
400:       Vec         xstar;
401:       PetscReal   norm;

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

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

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

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

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

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