Actual source code: ex30.c

petsc-3.4.2 2013-07-02
  2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
  3: It is copied and intended to move dirty codes from ksp/examples/tutorials/ex10.c and simplify ex10.c.\n\
  4:   Input parameters include\n\
  5:   -f0 <input_file> : first file to load (small system)\n\
  6:   -f1 <input_file> : second file to load (larger system)\n\n\
  7:   -trans  : solve transpose system instead\n\n";
  8: /*
  9:   This code  can be used to test PETSc interface to other packages.\n\
 10:   Examples of command line options:       \n\
 11:    ex30 -f0 <datafile> -ksp_type preonly  \n\
 12:         -help -ksp_view                  \n\
 13:         -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
 14:         -ksp_type preonly -pc_type lu -pc_factor_mat_solver_package superlu or superlu_dist or mumps \n\
 15:         -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_package mumps \n\
 16:    mpiexec -n <np> ex30 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij

 18:    ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package superlu -mat_superlu_conditionnumber -ckerror -mat_superlu_diagpivotthresh 0
 19:    ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type hypre -pc_hypre_type boomeramg -ksp_type fgmres -ckError
 20:    ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package petsc -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5 -ckerror
 21:  \n\n";
 22: */
 23: /*T
 24:    Concepts: KSP solving a linear system
 25:    Processors: n
 26: T*/

 28: #include <petscksp.h>
 29: #include <petsctime.h>

 33: int main(int argc,char **args)
 34: {
 35:   KSP            ksp;
 36:   Mat            A,B;
 37:   Vec            x,b,u,b2;        /* approx solution, RHS, exact solution */
 38:   PetscViewer    fd;              /* viewer */
 39:   char           file[4][PETSC_MAX_PATH_LEN];     /* input file name */
 40:   PetscBool      table     = PETSC_FALSE,flg,flgB=PETSC_FALSE,trans=PETSC_FALSE,partition=PETSC_FALSE,initialguess = PETSC_FALSE;
 41:   PetscBool      outputSoln=PETSC_FALSE;
 43:   PetscInt       its,num_numfac,n,M;
 44:   PetscReal      rnorm,enorm;
 45:   PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
 46:   PetscBool      preload=PETSC_TRUE,diagonalscale,isSymmetric,ckrnorm=PETSC_TRUE,Test_MatDuplicate=PETSC_FALSE,ckerror=PETSC_FALSE;
 47:   PetscMPIInt    rank;
 48:   PetscScalar    sigma;
 49:   PetscInt       m;

 51:   PetscInitialize(&argc,&args,(char*)0,help);
 52:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 53:   PetscOptionsGetBool(NULL,"-table",&table,NULL);
 54:   PetscOptionsGetBool(NULL,"-trans",&trans,NULL);
 55:   PetscOptionsGetBool(NULL,"-partition",&partition,NULL);
 56:   PetscOptionsGetBool(NULL,"-initialguess",&initialguess,NULL);
 57:   PetscOptionsGetBool(NULL,"-output_solution",&outputSoln,NULL);
 58:   PetscOptionsGetBool(NULL,"-ckrnorm",&ckrnorm,NULL);
 59:   PetscOptionsGetBool(NULL,"-ckerror",&ckerror,NULL);

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

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

 90:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 91:                          Load system
 92:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

100:   /*
101:      Load the matrix and vector; then destroy the viewer.
102:   */
103:   MatCreate(PETSC_COMM_WORLD,&A);
104:   MatLoad(A,fd);

106:   if (!preload) {
107:     flg  = PETSC_FALSE;
108:     PetscOptionsGetString(NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN,&flg);
109:     if (flg) {   /* rhs is stored in a separate file */
110:       PetscViewerDestroy(&fd);
111:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
112:     } else {
113:       /* if file contains no RHS, then use a vector of all ones */
114:       PetscInfo(0,"Using vector of ones for RHS\n");
115:       MatGetLocalSize(A,&m,NULL);
116:       VecCreate(PETSC_COMM_WORLD,&b);
117:       VecSetSizes(b,m,PETSC_DECIDE);
118:       VecSetFromOptions(b);
119:       VecSet(b,1.0);
120:       PetscObjectSetName((PetscObject)b, "Rhs vector");
121:     }
122:   }
123:   PetscViewerDestroy(&fd);

125:   /* Test MatDuplicate() */
126:   if (Test_MatDuplicate) {
127:     MatDuplicate(A,MAT_COPY_VALUES,&B);
128:     MatEqual(A,B,&flg);
129:     if (!flg) {
130:       PetscPrintf(PETSC_COMM_WORLD,"  A != B \n");
131:     }
132:     MatDestroy(&B);
133:   }

135:   /* Add a shift to A */
136:   PetscOptionsGetScalar(NULL,"-mat_sigma",&sigma,&flg);
137:   if (flg) {
138:     PetscOptionsGetString(NULL,"-fB",file[2],PETSC_MAX_PATH_LEN,&flgB);
139:     if (flgB) {
140:       /* load B to get A = A + sigma*B */
141:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
142:       MatCreate(PETSC_COMM_WORLD,&B);
143:       MatSetOptionsPrefix(B,"B_");
144:       MatLoad(B,fd);
145:       PetscViewerDestroy(&fd);
146:       MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN);   /* A <- sigma*B + A */
147:     } else {
148:       MatShift(A,sigma);
149:     }
150:   }

152:   /* Make A singular for testing zero-pivot of ilu factorization        */
153:   /* Example: ./ex30 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
154:   flg  = PETSC_FALSE;
155:   PetscOptionsGetBool(NULL, "-test_zeropivot", &flg,NULL);
156:   if (flg) {
157:     PetscInt          row,ncols;
158:     const PetscInt    *cols;
159:     const PetscScalar *vals;
160:     PetscBool         flg1=PETSC_FALSE;
161:     PetscScalar       *zeros;
162:     row  = 0;
163:     MatGetRow(A,row,&ncols,&cols,&vals);
164:     PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
165:     PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
166:     flg1 = PETSC_FALSE;
167:     PetscOptionsGetBool(NULL, "-set_row_zero", &flg1,NULL);
168:     if (flg1) {   /* set entire row as zero */
169:       MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
170:     } else {   /* only set (row,row) entry as zero */
171:       MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
172:     }
173:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
174:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
175:   }

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

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

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

206:     VecCreate(PETSC_COMM_WORLD,&tmp);
207:     VecSetSizes(tmp,n,PETSC_DECIDE);
208:     VecSetFromOptions(tmp);
209:     VecGetOwnershipRange(b,&start,&end);
210:     VecGetLocalSize(b,&mvec);
211:     VecGetArray(b,&bold);
212:     for (j=0; j<mvec; j++) {
213:       indx = start+j;
214:       VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
215:     }
216:     VecRestoreArray(b,&bold);
217:     VecDestroy(&b);
218:     VecAssemblyBegin(tmp);
219:     VecAssemblyEnd(tmp);
220:     b    = tmp;
221:   }
222:   VecDuplicate(b,&b2);
223:   VecDuplicate(b,&x);
224:   PetscObjectSetName((PetscObject)x, "Solution vector");
225:   VecDuplicate(b,&u);
226:   PetscObjectSetName((PetscObject)u, "True Solution vector");
227:   VecSet(x,0.0);

229:   if (ckerror) {   /* Set true solution */
230:     VecSet(u,1.0);
231:     MatMult(A,u,b);
232:   }

234:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
235:                     Setup solve for system
236:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

238:   if (partition) {
239:     MatPartitioning mpart;
240:     IS              mis,nis,is;
241:     PetscInt        *count;
242:     PetscMPIInt     size;
243:     Mat             BB;
244:     MPI_Comm_size(PETSC_COMM_WORLD,&size);
245:     MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
246:     PetscMalloc(size*sizeof(PetscInt),&count);
247:     MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
248:     MatPartitioningSetAdjacency(mpart, A);
249:     /* MatPartitioningSetVertexWeights(mpart, weight); */
250:     MatPartitioningSetFromOptions(mpart);
251:     MatPartitioningApply(mpart, &mis);
252:     MatPartitioningDestroy(&mpart);
253:     ISPartitioningToNumbering(mis,&nis);
254:     ISPartitioningCount(mis,size,count);
255:     ISDestroy(&mis);
256:     ISInvertPermutation(nis, count[rank], &is);
257:     PetscFree(count);
258:     ISDestroy(&nis);
259:     ISSort(is);
260:     MatGetSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&BB);

262:     /* need to move the vector also */
263:     ISDestroy(&is);
264:     MatDestroy(&A);
265:     A    = BB;
266:   }

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--) {

282:     KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
283:     KSPSetFromOptions(ksp);

285:     /*
286:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
287:      enable more precise profiling of setting up the preconditioner.
288:      These calls are optional, since both will be called within
289:      KSPSolve() if they haven't been called already.
290:     */
291:     KSPSetUp(ksp);
292:     KSPSetUpOnBlocks(ksp);
293:     PetscTime(&tsetup2);
294:     tsetup = tsetup2 - tsetup1;

296:     /*
297:      Tests "diagonal-scaling of preconditioned residual norm" as used
298:      by many ODE integrator codes including SUNDIALS. Note this is different
299:      than diagonally scaling the matrix before computing the preconditioner
300:     */
301:     diagonalscale = PETSC_FALSE;
302:     PetscOptionsGetBool(NULL,"-diagonal_scale",&diagonalscale,NULL);
303:     if (diagonalscale) {
304:       PC       pc;
305:       PetscInt j,start,end,n;
306:       Vec      scale;

308:       KSPGetPC(ksp,&pc);
309:       VecGetSize(x,&n);
310:       VecDuplicate(x,&scale);
311:       VecGetOwnershipRange(scale,&start,&end);
312:       for (j=start; j<end; j++) {
313:         VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
314:       }
315:       VecAssemblyBegin(scale);
316:       VecAssemblyEnd(scale);
317:       PCSetDiagonalScale(pc,scale);
318:       VecDestroy(&scale);
319:     }

321:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
322:                          Solve system
323:       - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
324:     /*
325:      Solve linear system; we also explicitly time this stage.
326:     */
327:     PetscTime(&tsolve1);
328:     if (trans) {
329:       KSPSolveTranspose(ksp,b,x);
330:       KSPGetIterationNumber(ksp,&its);
331:     } else {
332:       PetscInt num_rhs=1;
333:       PetscOptionsGetInt(NULL,"-num_rhs",&num_rhs,NULL);

335:       while (num_rhs--) {
336:         KSPSolve(ksp,b,x);
337:       }
338:       KSPGetIterationNumber(ksp,&its);
339:       if (ckrnorm) {     /* Check residual for each rhs */
340:         if (trans) {
341:           MatMultTranspose(A,x,b2);
342:         } else {
343:           MatMult(A,x,b2);
344:         }
345:         VecAXPY(b2,-1.0,b);
346:         VecNorm(b2,NORM_2,&rnorm);
347:         PetscPrintf(PETSC_COMM_WORLD,"  Number of iterations = %3D\n",its);
348:         PetscPrintf(PETSC_COMM_WORLD,"  Residual norm %G\n",rnorm);
349:       }
350:       if (ckerror && !trans) {    /* Check error for each rhs */
351:         /* VecView(x,PETSC_VIEWER_STDOUT_WORLD); */
352:         VecAXPY(u,-1.0,x);
353:         VecNorm(u,NORM_2,&enorm);
354:         PetscPrintf(PETSC_COMM_WORLD,"  Error norm %G\n",enorm);
355:       }

357:     }   /* while (num_rhs--) */
358:     PetscTime(&tsolve2);
359:     tsolve = tsolve2 - tsolve1;


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 && ckrnorm) {
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,rnorm,tsetup+tsolve,tsetup,tsolve,kspinfo);

381:       /*
382:         Destroy the viewer
383:       */
384:       PetscViewerDestroy(&viewer);
385:     }

387:     PetscOptionsGetString(NULL,"-solution",file[3],PETSC_MAX_PATH_LEN,&flg);
388:     if (flg) {
389:       PetscViewer viewer;
390:       Vec         xstar;

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

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

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

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

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

433:   PetscFinalize();
434:   return 0;
435: }