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: }