Actual source code: ex10.c

petsc-master 2014-12-21
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:     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,"-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:     PetscCalloc1(ncols+1,&zeros);
164:     PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);
165:     if (flg1) {   /* set entire row as zero */
166:       MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
167:     } else {   /* only set (row,row) entry as zero */
168:       MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
169:     }
170:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
171:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
172:   }

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

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

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

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

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


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

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

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

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

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

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

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

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