Actual source code: ex10.c

  1: /*$Id: ex10.c,v 1.58 2001/09/11 16:33:29 bsmith Exp $*/

  3: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
  4: This version first preloads and solves a small system, then loads \n\
  5: another (larger) system and solves it as well.  This example illustrates\n\
  6: preloading of instructions with the smaller system so that more accurate\n\
  7: performance monitoring can be done with the larger one (that actually\n\
  8: is the system of interest).  See the 'Performance Hints' chapter of the\n\
  9: users manual for a discussion of preloading.  Input parameters include\n\
 10:   -f0 <input_file> : first file to load (small system)\n\
 11:   -f1 <input_file> : second file to load (larger system)\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 -matload_type seqaijspooles/superlu/superlu_dist/aijmumps \n\
 20:         -ksp_type preonly -pc_type cholesky -matload_type seqsbaijspooles/dscpack/sbaijmumps    \n\
 21:         -f0 <A> -fB <B> -matload_type sbaijmumps -ksp_type preonly -pc_type cholesky -test_inertia -mat_sigma <sigma> \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:      petsc.h       - base PETSc routines   petscvec.h - vectors
 32:      petscsys.h    - system routines       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


 41: int main(int argc,char **args)
 42: {
 43:   KSP            ksp;             /* linear solver context */
 44:   Mat            A,B;            /* matrix */
 45:   Vec            x,b,u;          /* approx solution, RHS, exact solution */
 46:   PetscViewer    fd;               /* viewer */
 47:   char           file[3][128];     /* input file name */
 48:   PetscTruth     table,flg,flgB=PETSC_FALSE,trans=PETSC_FALSE,partition=PETSC_FALSE;
 49:   int            ierr,its,ierrp;
 50:   PetscReal      norm;
 51:   PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
 52:   PetscScalar    zero = 0.0,none = -1.0;
 53:   PetscTruth     preload=PETSC_TRUE,diagonalscale,hasNullSpace,isSymmetric,cknorm=PETSC_FALSE;
 54:   int            num_numfac;
 55:   PetscScalar    sigma;

 57:   PetscInitialize(&argc,&args,(char *)0,help);

 59:   PetscOptionsHasName(PETSC_NULL,"-table",&table);
 60:   PetscOptionsHasName(PETSC_NULL,"-trans",&trans);
 61:   PetscOptionsHasName(PETSC_NULL,"-partition",&partition);

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

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

 92:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 93:                            Load system
 94:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 96:     /* 
 97:        Open binary file.  Note that we use PETSC_FILE_RDONLY to indicate
 98:        reading from this file.
 99:     */
100:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PreLoadIt],PETSC_FILE_RDONLY,&fd);

102:     /*
103:        Load the matrix and vector; then destroy the viewer.
104:     */
105:     MatLoad(fd,MATAIJ,&A);
106:     PetscPushErrorHandler(PetscIgnoreErrorHandler,PETSC_NULL);
107:     ierrp = VecLoad(fd,PETSC_NULL,&b);
108:     PetscPopErrorHandler();
109:     if (ierrp) { /* if file contains no RHS, then use a vector of all ones */
110:       int         m;
111:       PetscScalar one = 1.0;
112:       PetscLogInfo(0,"Using vector of ones for RHS\n");
113:       MatGetLocalSize(A,&m,PETSC_NULL);
114:       VecCreate(PETSC_COMM_WORLD,&b);
115:       VecSetSizes(b,m,PETSC_DECIDE);
116:       VecSetFromOptions(b);
117:       VecSet(&one,b);
118:     }
119:     PetscViewerDestroy(fd);

121:     /* Add a shift to A */
122:     PetscOptionsGetScalar(PETSC_NULL,"-mat_sigma",&sigma,&flg);
123:     if(flg) {
124:       PetscOptionsGetString(PETSC_NULL,"-fB",file[2],127,&flgB);
125:       if (flgB){
126:         /* load B to get A = A + sigma*B */
127:         PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],PETSC_FILE_RDONLY,&fd);
128:         MatLoad(fd,MATAIJ,&B);
129:         PetscViewerDestroy(fd);
130:         MatAXPY(&sigma,B,A,DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
131:       } else {
132:         MatShift(&sigma,A);
133:       }
134:     }

136:     /* Check whether A is symmetric */
137:     PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
138:     if (flg) {
139:       Mat Atrans;
140:       MatTranspose(A, &Atrans);
141:       MatEqual(A, Atrans, &isSymmetric);
142:       if (isSymmetric) {
143:         PetscPrintf(PETSC_COMM_WORLD,"A is symmetric \n");
144:       } else {
145:         PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric \n");
146:       }
147:       MatDestroy(Atrans);
148:     }

150:     /* 
151:        If the loaded matrix is larger than the vector (due to being padded 
152:        to match the block size of the system), then create a new padded vector.
153:     */
154:     {
155:       int         m,n,j,mvec,start,end,indx;
156:       Vec         tmp;
157:       PetscScalar *bold;

159:       /* Create a new vector b by padding the old one */
160:       MatGetLocalSize(A,&m,&n);
161:       VecCreate(PETSC_COMM_WORLD,&tmp);
162:       VecSetSizes(tmp,m,PETSC_DECIDE);
163:       VecSetFromOptions(tmp);
164:       VecGetOwnershipRange(b,&start,&end);
165:       VecGetLocalSize(b,&mvec);
166:       VecGetArray(b,&bold);
167:       for (j=0; j<mvec; j++) {
168:         indx = start+j;
169:         VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
170:       }
171:       VecRestoreArray(b,&bold);
172:       VecDestroy(b);
173:       VecAssemblyBegin(tmp);
174:       VecAssemblyEnd(tmp);
175:       b = tmp;
176:     }
177:     VecDuplicate(b,&x);
178:     VecDuplicate(b,&u);
179:     VecSet(&zero,x);

181:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
182:                       Setup solve for system
183:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


186:     if (partition) {
187:       MatPartitioning mpart;
188:       IS              mis,nis,isn,is;
189:       int             *count,size,rank;
190:       Mat             BB;
191:       MPI_Comm_size(PETSC_COMM_WORLD,&size);
192:       MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
193:       PetscMalloc(size*sizeof(int),&count);
194:       MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
195:       MatPartitioningSetAdjacency(mpart, A);
196:       /* MatPartitioningSetVertexWeights(mpart, weight); */
197:       MatPartitioningSetFromOptions(mpart);
198:       MatPartitioningApply(mpart, &mis);
199:       MatPartitioningDestroy(mpart);
200:       ISPartitioningToNumbering(mis,&nis);
201:       ISPartitioningCount(mis,count);
202:       ISDestroy(mis);
203:       ISInvertPermutation(nis, count[rank], &is);
204:       PetscFree(count);
205:       ISDestroy(nis);
206:       ISSort(is);
207:       ISAllGather(is,&isn);
208:       MatGetSubMatrix(A,is,isn,PETSC_DECIDE,MAT_INITIAL_MATRIX,&BB);

210:       /* need to move the vector also */
211:       ISDestroy(is);
212:       ISDestroy(isn);
213:       MatDestroy(A);
214:       A    = BB;
215:     }
216: 
217:     /*
218:        Conclude profiling last stage; begin profiling next stage.
219:     */
220:     PreLoadStage("KSPSetUp");

222:     /*
223:        We also explicitly time this stage via PetscGetTime()
224:     */
225:     PetscGetTime(&tsetup1);

227:     /*
228:        Create linear solver; set operators; set runtime options.
229:     */
230:     KSPCreate(PETSC_COMM_WORLD,&ksp);

232:     num_numfac = 1;
233:     PetscOptionsGetInt(PETSC_NULL,"-num_numfac",&num_numfac,PETSC_NULL);
234:     while ( num_numfac-- ){
235:       /* KSPSetOperators(ksp,A,A,DIFFERENT_NONZERO_PATTERN); */
236:     KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
237:     KSPSetFromOptions(ksp);

239:     /* 
240:        Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
241:        enable more precise profiling of setting up the preconditioner.
242:        These calls are optional, since both will be called within
243:        KSPSolve() if they haven't been called already.
244:     */
245:     KSPSetRhs(ksp,b);
246:     KSPSetSolution(ksp,x);
247:     KSPSetUp(ksp);
248:     KSPSetUpOnBlocks(ksp);
249:     PetscGetTime(&tsetup2);
250:     tsetup = tsetup2 - tsetup1;

252:     /*
253:       Test MatGetInertia()
254:       Usage:
255:       ex10 -f0 <mat_binaryfile> -ksp_type preonly -pc_type cholesky -mat_type seqsbaij -test_inertia -mat_sigma <sigma>
256:      */
257:     PetscOptionsHasName(PETSC_NULL,"-test_inertia",&flg);
258:     if (flg){
259:       PC      pc;
260:       int     nneg, nzero, npos;
261:       Mat     F;
262: 
263:       KSPGetPC(ksp,&pc);
264:       PCGetFactoredMatrix(pc,&F);
265:       MatGetInertia(F,&nneg,&nzero,&npos);
266:       PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %d, nzero: %d, npos: %d\n",nneg,nzero,npos);
267:     }

269:     /*
270:        Tests "diagonal-scaling of preconditioned residual norm" as used 
271:        by many ODE integrator codes including PVODE. Note this is different
272:        than diagonally scaling the matrix before computing the preconditioner
273:     */
274:     PetscOptionsHasName(PETSC_NULL,"-diagonal_scale",&diagonalscale);
275:     if (diagonalscale) {
276:       PC     pc;
277:       int    j,start,end,n;
278:       Vec    scale;
279: 
280:       KSPGetPC(ksp,&pc);
281:       VecGetSize(x,&n);
282:       VecDuplicate(x,&scale);
283:       VecGetOwnershipRange(scale,&start,&end);
284:       for (j=start; j<end; j++) {
285:         VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
286:       }
287:       VecAssemblyBegin(scale);
288:       VecAssemblyEnd(scale);
289:       PCDiagonalScaleSet(pc,scale);
290:       VecDestroy(scale);

292:     }

294:     PetscOptionsHasName(PETSC_NULL, "-null_space", &hasNullSpace);
295:     if (hasNullSpace == PETSC_TRUE) {
296:       MatNullSpace nullSpace;
297:       PC           pc;

299:       MatNullSpaceCreate(PETSC_COMM_WORLD, 1, 0, PETSC_NULL, &nullSpace);
300:       MatNullSpaceTest(nullSpace, A);
301:       KSPGetPC(ksp,&pc);
302:       PCNullSpaceAttach(pc, nullSpace);
303:       MatNullSpaceDestroy(nullSpace);
304:     }

306:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
307:                            Solve system
308:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

310:     /*
311:        Begin profiling next stage
312:     */
313:     PreLoadStage("KSPSolve");

315:     /*
316:        Solve linear system; we also explicitly time this stage.
317:     */
318:     PetscGetTime(&tsolve1);
319:     if (trans) {
320:       KSPSolveTranspose(ksp);
321:     } else {
322:       int  num_rhs=1;
323:       PetscOptionsGetInt(PETSC_NULL,"-num_rhs",&num_rhs,PETSC_NULL);
324:       PetscOptionsHasName(PETSC_NULL,"-cknorm",&cknorm);
325:       while ( num_rhs-- ) {
326:         KSPSolve(ksp);
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(&none,b,u);
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:    /* 
345:        Conclude profiling this stage
346:     */
347:     PreLoadStage("Cleanup");

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

353:     /* 
354:        Check error
355:     */
356:     if (trans) {
357:       MatMultTranspose(A,x,u);
358:     } else {
359:       MatMult(A,x,u);
360:     }
361:     VecAXPY(&none,b,u);
362:     VecNorm(u,NORM_2,&norm);

364:     /*
365:        Write output (optinally using table for solver details).
366:         - PetscPrintf() handles output for multiprocessor jobs 
367:           by printing from only one processor in the communicator.
368:         - KSPView() prints information about the linear solver.
369:     */
370:     if (table) {
371:       char        *matrixname,kspinfo[120];
372:       PetscViewer viewer;

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

383:       /*
384:          Destroy the viewer
385:       */
386:       PetscViewerDestroy(viewer);
387:     } else {
388:       PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
389:       PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);
390:     }

392:     PetscOptionsHasName(PETSC_NULL, "-ksp_reason", &flg);
393:     if (flg){
394:       KSPConvergedReason reason;
395:       KSPGetConvergedReason(ksp,&reason);
396:       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %d\n", reason);
397:     }
398: 
399:     } /* while ( num_numfac-- ) */

401:     /* 
402:        Free work space.  All PETSc objects should be destroyed when they
403:        are no longer needed.
404:     */
405:     MatDestroy(A); VecDestroy(b);
406:     VecDestroy(u); VecDestroy(x);
407:     KSPDestroy(ksp);
408:     if (flgB) { MatDestroy(B); }
409:   PreLoadEnd();
410:   /* -----------------------------------------------------------
411:                       End of linear solver loop
412:      ----------------------------------------------------------- */

414:   PetscFinalize();
415:   return 0;
416: }