Actual source code: ex10.c

  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 spooles or superlu or superlu_dist or mumps \n\
 19:         -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package spooles or 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>

 40: int main(int argc,char **args)
 41: {
 42:   KSP            ksp;             /* linear solver context */
 43:   Mat            A;               /* matrix */
 44:   Vec            x,b,u;           /* approx solution, RHS, exact solution */
 45:   PetscViewer    fd;              /* viewer */
 46:   char           file[4][PETSC_MAX_PATH_LEN];     /* input file name */
 47:   PetscBool      table=PETSC_FALSE,flg,trans=PETSC_FALSE,initialguess = PETSC_FALSE;
 48:   PetscBool      outputSoln=PETSC_FALSE;
 50:   PetscInt       its,num_numfac,m,n,M;
 51:   PetscReal      norm;
 52:   PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
 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(PETSC_NULL,"-table",&table,PETSC_NULL);
 60:   PetscOptionsGetBool(PETSC_NULL,"-trans",&trans,PETSC_NULL);
 61:   PetscOptionsGetBool(PETSC_NULL,"-initialguess",&initialguess,PETSC_NULL);
 62:   PetscOptionsGetBool(PETSC_NULL,"-output_solution",&outputSoln,PETSC_NULL);
 63:   PetscOptionsGetString(PETSC_NULL,"-initialguessfilename",initialguessfilename,PETSC_MAX_PATH_LEN,&initialguessfile);

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

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

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

 98:     /* 
 99:        Open binary file.  Note that we use FILE_MODE_READ to indicate
100:        reading from this file.
101:     */
102:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);
103: 
104:     /*
105:        Load the matrix and vector; then destroy the viewer.
106:     */
107:     MatCreate(PETSC_COMM_WORLD,&A);
108:     MatSetFromOptions(A);
109:     MatLoad(A,fd);

111:     if (!preload){
112:       flg = PETSC_FALSE;
113:       PetscOptionsGetString(PETSC_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,PETSC_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 {
131:         VecSetFromOptions(b);
132:         VecLoad(b,fd);
133:       }
134:     }
135:     PetscViewerDestroy(&fd);

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

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

176:     /* 
177:        If the loaded matrix is larger than the vector (due to being padded 
178:        to match the block size of the system), then create a new padded vector.
179:     */
180: 
181:     MatGetLocalSize(A,&m,&n);
182:     if (m != n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
183:     MatGetSize(A,&M,PETSC_NULL);
184:     VecGetSize(b,&m);
185:     if (M != m) { /* Create a new vector b by padding the old one */
186:       PetscInt    j,mvec,start,end,indx;
187:       Vec         tmp;
188:       PetscScalar *bold;

190:       VecCreate(PETSC_COMM_WORLD,&tmp);
191:       VecSetSizes(tmp,n,PETSC_DECIDE);
192:       VecSetFromOptions(tmp);
193:       VecGetOwnershipRange(b,&start,&end);
194:       VecGetLocalSize(b,&mvec);
195:       VecGetArray(b,&bold);
196:       for (j=0; j<mvec; j++) {
197:         indx = start+j;
198:         VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
199:       }
200:       VecRestoreArray(b,&bold);
201:       VecDestroy(&b);
202:       VecAssemblyBegin(tmp);
203:       VecAssemblyEnd(tmp);
204:       b = tmp;
205:     }
206:     VecDuplicate(b,&x);
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(PETSC_NULL, "-check_scaling", &flg,PETSC_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, PETSC_NULL);
232:       MatGetRowMinAbs(A, min, PETSC_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 PetscGetTime()
270:     */
271:     PetscGetTime(&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(PETSC_NULL,"-num_numfac",&num_numfac,PETSC_NULL);
280:     while ( num_numfac-- ){
281:       PetscBool  lsqr;
282:       char       str[32];
283:       PetscOptionsGetString(PETSC_NULL,"-ksp_type",str,32,&lsqr);
284:       if (lsqr) {
285:         PetscStrcmp("lsqr",str,&lsqr);
286:       }
287:       if (lsqr) {
288:         Mat BtB;
289:         MatMatMultTranspose(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:       PetscGetTime(&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:       PetscGetTime(&tsolve1);
316:       if (trans) {
317:         KSPSolveTranspose(ksp,b,x);
318:         KSPGetIterationNumber(ksp,&its);
319:       } else {
320:         PetscInt  num_rhs=1;
321:         PetscOptionsGetInt(PETSC_NULL,"-num_rhs",&num_rhs,PETSC_NULL);
322:         cknorm = PETSC_FALSE;
323:         PetscOptionsGetBool(PETSC_NULL,"-cknorm",&cknorm,PETSC_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:           PetscPrintf(PETSC_COMM_WORLD,"  Residual norm %A\n",norm);
339:         }
340:       } /* while ( num_rhs-- ) */
341:       PetscGetTime(&tsolve2);
342:       tsolve = tsolve2 - tsolve1;

344:       /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
345:             Check error, print output, free data structures.
346:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

377:         /*
378:           Destroy the viewer
379:         */
380:         PetscViewerDestroy(&viewer);
381:       } else {
382:         PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
383:         PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);
384:       }
385:       PetscOptionsGetString(PETSC_NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
386:       if (flg) {
387:         PetscViewer viewer;
388:         Vec         xstar;
389:         PetscReal   norm;

391:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
392:         VecCreate(PETSC_COMM_WORLD,&xstar);
393:         VecLoad(xstar,viewer);
394:         VecAXPY(xstar, -1.0, x);
395:         VecNorm(xstar, NORM_2, &norm);
396:         PetscPrintf(PETSC_COMM_WORLD, "Error norm %A\n", norm);
397:         VecDestroy(&xstar);
398:         PetscViewerDestroy(&viewer);
399:       }
400:       if (outputSoln) {
401:         PetscViewer viewer;

403:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
404:         VecView(x, viewer);
405:         PetscViewerDestroy(&viewer);
406:       }

408:       flg  = PETSC_FALSE;
409:       PetscOptionsGetBool(PETSC_NULL, "-ksp_reason", &flg,PETSC_NULL);
410:       if (flg){
411:         KSPConvergedReason reason;
412:         KSPGetConvergedReason(ksp,&reason);
413:         PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
414:       }
415: 
416:     } /* while ( num_numfac-- ) */

418:     /* 
419:        Free work space.  All PETSc objects should be destroyed when they
420:        are no longer needed.
421:     */
422:     MatDestroy(&A); VecDestroy(&b);
423:     VecDestroy(&u); VecDestroy(&x);
424:     KSPDestroy(&ksp);
425:   PetscPreLoadEnd();
426:   /* -----------------------------------------------------------
427:                       End of linear solver loop
428:      ----------------------------------------------------------- */

430:   PetscFinalize();
431:   return 0;
432: }