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 -sles_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\n";
 21: */
 22: /*T
 23:    Concepts: SLES^solving a linear system
 24:    Processors: n
 25: T*/

 27: /* 
 28:   Include "petscsles.h" so that we can use SLES solvers.  Note that this file
 29:   automatically includes:
 30:      petsc.h       - base PETSc routines   petscvec.h - vectors
 31:      petscsys.h    - system routines       petscmat.h - matrices
 32:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 33:      petscviewer.h - viewers               petscpc.h  - preconditioners
 34: */
 35:  #include petscsles.h


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

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

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

 62:   /* 
 63:      Determine files from which we read the two linear systems
 64:      (matrix and right-hand-side vector).
 65:   */
 66:   PetscOptionsGetString(PETSC_NULL,"-f",file[0],127,&flg);
 67:   if (flg) {
 68:     PetscStrcpy(file[1],file[0]);
 69:   } else {
 70:     PetscOptionsGetString(PETSC_NULL,"-f0",file[0],127,&flg);
 71:     if (!flg) SETERRQ(1,"Must indicate binary file with the -f0 or -f option");
 72:     PetscOptionsGetString(PETSC_NULL,"-f1",file[1],127,&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:   PreLoadBegin(preload,"Load system");

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

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

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

119:     /* Add a shift to A */
120:     PetscOptionsGetScalar(PETSC_NULL,"-mat_sigma",&sigma,&flg);
121:     if(flg) {MatShift(&sigma,A); }

123:     /* Check whether A is symmetric */
124:     PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
125:     if (flg) {
126:       Mat Atrans;
127:       MatTranspose(A, &Atrans);
128:       MatEqual(A, Atrans, &isSymmetric);
129:       if (isSymmetric) {
130:         PetscPrintf(PETSC_COMM_WORLD,"A is symmetric \n");
131:       } else {
132:         PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric \n");
133:       }
134:       MatDestroy(Atrans);
135:     }

137:     /* 
138:        If the loaded matrix is larger than the vector (due to being padded 
139:        to match the block size of the system), then create a new padded vector.
140:     */
141:     {
142:       int         m,n,j,mvec,start,end,indx;
143:       Vec         tmp;
144:       PetscScalar *bold;

146:       /* Create a new vector b by padding the old one */
147:       MatGetLocalSize(A,&m,&n);
148:       VecCreate(PETSC_COMM_WORLD,&tmp);
149:       VecSetSizes(tmp,m,PETSC_DECIDE);
150:       VecSetFromOptions(tmp);
151:       VecGetOwnershipRange(b,&start,&end);
152:       VecGetLocalSize(b,&mvec);
153:       VecGetArray(b,&bold);
154:       for (j=0; j<mvec; j++) {
155:         indx = start+j;
156:         VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
157:       }
158:       VecRestoreArray(b,&bold);
159:       VecDestroy(b);
160:       VecAssemblyBegin(tmp);
161:       VecAssemblyEnd(tmp);
162:       b = tmp;
163:     }
164:     VecDuplicate(b,&x);
165:     VecDuplicate(b,&u);
166:     VecSet(&zero,x);

168:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
169:                       Setup solve for system
170:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */


173:     if (partition) {
174:       MatPartitioning mpart;
175:       IS              mis,nis,isn,is;
176:       int             *count,size,rank;
177:       Mat             B;
178:       MPI_Comm_size(PETSC_COMM_WORLD,&size);
179:       MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
180:       PetscMalloc(size*sizeof(int),&count);
181:       MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
182:       MatPartitioningSetAdjacency(mpart, A);
183:       /* MatPartitioningSetVertexWeights(mpart, weight); */
184:       MatPartitioningSetFromOptions(mpart);
185:       MatPartitioningApply(mpart, &mis);
186:       MatPartitioningDestroy(mpart);
187:       ISPartitioningToNumbering(mis,&nis);
188:       ISPartitioningCount(mis,count);
189:       ISDestroy(mis);
190:       ISInvertPermutation(nis, count[rank], &is);
191:       PetscFree(count);
192:       ISDestroy(nis);
193:       ISSort(is);
194:       ISAllGather(is,&isn);
195:       MatGetSubMatrix(A,is,isn,PETSC_DECIDE,MAT_INITIAL_MATRIX,&B);

197:       /* need to move the vector also */
198:       ISDestroy(is);
199:       ISDestroy(isn);
200:       MatDestroy(A);
201:       A    = B;
202:     }
203: 
204:     /*
205:        Conclude profiling last stage; begin profiling next stage.
206:     */
207:     PreLoadStage("SLESSetUp");

209:     /*
210:        We also explicitly time this stage via PetscGetTime()
211:     */
212:     PetscGetTime(&tsetup1);

214:     /*
215:        Create linear solver; set operators; set runtime options.
216:     */
217:     SLESCreate(PETSC_COMM_WORLD,&sles);

219:     num_numfac = 1;
220:     PetscOptionsGetInt(PETSC_NULL,"-num_numfac",&num_numfac,PETSC_NULL);
221:     while ( num_numfac-- ){
222:       /* SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN); */
223:     SLESSetOperators(sles,A,A,SAME_NONZERO_PATTERN);
224:     SLESSetFromOptions(sles);

226:     /* 
227:        Here we explicitly call SLESSetUp() and SLESSetUpOnBlocks() to
228:        enable more precise profiling of setting up the preconditioner.
229:        These calls are optional, since both will be called within
230:        SLESSolve() if they haven't been called already.
231:     */
232:     SLESSetUp(sles,b,x);
233:     SLESSetUpOnBlocks(sles);
234:     PetscGetTime(&tsetup2);
235:     tsetup = tsetup2 - tsetup1;

237:     /*
238:        Tests "diagonal-scaling of preconditioned residual norm" as used 
239:        by many ODE integrator codes including PVODE. Note this is different
240:        than diagonally scaling the matrix before computing the preconditioner
241:     */
242:     PetscOptionsHasName(PETSC_NULL,"-diagonal_scale",&diagonalscale);
243:     if (diagonalscale) {
244:       PC     pc;
245:       int    j,start,end,n;
246:       Vec    scale;
247: 
248:       SLESGetPC(sles,&pc);
249:       VecGetSize(x,&n);
250:       VecDuplicate(x,&scale);
251:       VecGetOwnershipRange(scale,&start,&end);
252:       for (j=start; j<end; j++) {
253:         VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
254:       }
255:       VecAssemblyBegin(scale);
256:       VecAssemblyEnd(scale);
257:       PCDiagonalScaleSet(pc,scale);
258:       VecDestroy(scale);

260:     }

262:     PetscOptionsHasName(PETSC_NULL, "-null_space", &hasNullSpace);
263:     if (hasNullSpace == PETSC_TRUE) {
264:       MatNullSpace nullSpace;
265:       PC           pc;

267:       MatNullSpaceCreate(PETSC_COMM_WORLD, 1, 0, PETSC_NULL, &nullSpace);
268:       MatNullSpaceTest(nullSpace, A);
269:       SLESGetPC(sles,&pc);
270:       PCNullSpaceAttach(pc, nullSpace);
271:       MatNullSpaceDestroy(nullSpace);
272:     }

274:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
275:                            Solve system
276:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

278:     /*
279:        Begin profiling next stage
280:     */
281:     PreLoadStage("SLESSolve");

283:     /*
284:        Solve linear system; we also explicitly time this stage.
285:     */
286:     PetscGetTime(&tsolve1);
287:     if (trans) {
288:       SLESSolveTranspose(sles,b,x,&its);
289:     } else {
290:       int  num_rhs=1;
291:       PetscOptionsGetInt(PETSC_NULL,"-num_rhs",&num_rhs,PETSC_NULL);
292:       while ( num_rhs-- ) {
293:         SLESSolve(sles,b,x,&its);
294:       }
295:     }
296:     PetscGetTime(&tsolve2);
297:     tsolve = tsolve2 - tsolve1;

299:    /* 
300:        Conclude profiling this stage
301:     */
302:     PreLoadStage("Cleanup");

304:     /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
305:             Check error, print output, free data structures.
306:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

308:     /* 
309:        Check error
310:     */
311:     if (trans) {
312:       MatMultTranspose(A,x,u);
313:     } else {
314:       MatMult(A,x,u);
315:     }
316:     VecAXPY(&none,b,u);
317:     VecNorm(u,NORM_2,&norm);

319:     /*
320:        Write output (optinally using table for solver details).
321:         - PetscPrintf() handles output for multiprocessor jobs 
322:           by printing from only one processor in the communicator.
323:         - SLESView() prints information about the linear solver.
324:     */
325:     if (table) {
326:       char        *matrixname,slesinfo[120];
327:       PetscViewer viewer;

329:       /*
330:          Open a string viewer; then write info to it.
331:       */
332:       PetscViewerStringOpen(PETSC_COMM_WORLD,slesinfo,120,&viewer);
333:       SLESView(sles,viewer);
334:       PetscStrrchr(file[PreLoadIt],'/',&matrixname);
335:       PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3d %2.0e %2.1e %2.1e %2.1e %s \n",
336:                 matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,slesinfo);

338:       /*
339:          Destroy the viewer
340:       */
341:       PetscViewerDestroy(viewer);
342:     } else {
343:       PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3d\n",its);
344:       PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);
345:     }

347:     PetscOptionsHasName(PETSC_NULL, "-ksp_reason", &flg);
348:     if (flg){
349:       KSP ksp;
350:       KSPConvergedReason reason;
351:       SLESGetKSP(sles,&ksp);
352:       KSPGetConvergedReason(ksp,&reason);
353:       PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %d\n", reason);
354:     }
355: 
356:     } /* while ( num_numfac-- ) */

358:     /* 
359:        Free work space.  All PETSc objects should be destroyed when they
360:        are no longer needed.
361:     */
362:     MatDestroy(A); VecDestroy(b);
363:     VecDestroy(u); VecDestroy(x);
364:     SLESDestroy(sles);
365:   PreLoadEnd();
366:   /* -----------------------------------------------------------
367:                       End of linear solver loop
368:      ----------------------------------------------------------- */

370:   PetscFinalize();
371:   return 0;
372: }