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 illustratesn
6: preloading of instructions with the smaller system so that more accuraten
7: performance monitoring can be done with the larger one (that actuallyn
8: is the system of interest). See the 'Performance Hints' chapter of then
9: users manual for a discussion of preloading. Input parameters includen
10: -f0 <input_file> : first file to load (small system)n
11: -f1 <input_file> : second file to load (larger system)nn
12: -trans : solve transpose system insteadnn";
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: -pc_type lu -mat_aij_spooles/superlu/superlu_dist n
20: -pc_type cholesky -mat_baij_dscpack -matload_type mpibaij n
21: -pc_type cholesky -mat_sbaij_spooles -matload_type mpisbaijnn";
22: */
23: /*T
24: Concepts: SLES^solving a linear system
25: Processors: n
26: T*/
28: /*
29: Include "petscsles.h" so that we can use SLES 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 petscsles.h
39: #undef __FUNCT__
41: int main(int argc,char **args)
42: {
43: SLES sles; /* linear solver context */
44: Mat A; /* matrix */
45: Vec x,b,u; /* approx solution, RHS, exact solution */
46: PetscViewer fd; /* viewer */
47: char file[2][128]; /* input file name */
48: PetscTruth table,flg,trans,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;
54: int num_numfac;
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: ierr = MatLoad(fd,MATMPIAIJ,&A);
104: ierr = PetscPushErrorHandler(PetscIgnoreErrorHandler,PETSC_NULL);
105: ierrp = VecLoad(fd,&b);
106: ierr = 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 RHSn");
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: /* Check whether A is symmetric */
120: PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
121: if (flg) {
122: Mat Atrans;
123: MatTranspose(A, &Atrans);
124: MatEqual(A, Atrans, &isSymmetric);
125: if (isSymmetric) {
126: PetscPrintf(PETSC_COMM_WORLD,"A is symmetric n");
127: } else {
128: PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric n");
129: }
130: MatDestroy(Atrans);
131: }
133: /*
134: If the loaded matrix is larger than the vector (due to being padded
135: to match the block size of the system), then create a new padded vector.
136: */
137: {
138: int m,n,j,mvec,start,end,index;
139: Vec tmp;
140: PetscScalar *bold;
142: /* Create a new vector b by padding the old one */
143: MatGetLocalSize(A,&m,&n);
144: VecCreate(PETSC_COMM_WORLD,&tmp);
145: VecSetSizes(tmp,m,PETSC_DECIDE);
146: VecSetFromOptions(tmp);
147: VecGetOwnershipRange(b,&start,&end);
148: VecGetLocalSize(b,&mvec);
149: VecGetArray(b,&bold);
150: for (j=0; j<mvec; j++) {
151: index = start+j;
152: ierr = VecSetValues(tmp,1,&index,bold+j,INSERT_VALUES);
153: }
154: VecRestoreArray(b,&bold);
155: VecDestroy(b);
156: VecAssemblyBegin(tmp);
157: VecAssemblyEnd(tmp);
158: b = tmp;
159: }
160: VecDuplicate(b,&x);
161: VecDuplicate(b,&u);
162: VecSet(&zero,x);
164: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
165: Setup solve for system
166: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169: if (partition) {
170: MatPartitioning mpart;
171: IS mis,nis,isn,is;
172: int *count,size,rank;
173: Mat B;
174: MPI_Comm_size(PETSC_COMM_WORLD,&size);
175: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
176: PetscMalloc(size*sizeof(int),&count);
177: MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
178: MatPartitioningSetAdjacency(mpart, A);
179: /* MatPartitioningSetVertexWeights(mpart, weight); */
180: MatPartitioningSetFromOptions(mpart);
181: MatPartitioningApply(mpart, &mis);
182: MatPartitioningDestroy(mpart);
183: ISPartitioningToNumbering(mis,&nis);
184: ISPartitioningCount(mis,count);
185: ISDestroy(mis);
186: ISInvertPermutation(nis, count[rank], &is);
187: PetscFree(count);
188: ISDestroy(nis);
189: ISSort(is);
190: ISAllGather(is,&isn);
191: MatGetSubMatrix(A,is,isn,PETSC_DECIDE,MAT_INITIAL_MATRIX,&B);
193: /* need to move the vector also */
194: ISDestroy(is);
195: ISDestroy(isn);
196: MatDestroy(A);
197: A = B;
198: }
199:
200: /*
201: Conclude profiling last stage; begin profiling next stage.
202: */
203: PreLoadStage("SLESSetUp");
205: /*
206: We also explicitly time this stage via PetscGetTime()
207: */
208: PetscGetTime(&tsetup1);
210: /*
211: Create linear solver; set operators; set runtime options.
212: */
213: SLESCreate(PETSC_COMM_WORLD,&sles);
215: num_numfac = 1;
216: PetscOptionsGetInt(PETSC_NULL,"-num_numfac",&num_numfac,PETSC_NULL);
217: while ( num_numfac-- ){
218: /* SLESSetOperators(sles,A,A,DIFFERENT_NONZERO_PATTERN); */
219: SLESSetOperators(sles,A,A,SAME_NONZERO_PATTERN);
220: SLESSetFromOptions(sles);
222: /*
223: Here we explicitly call SLESSetUp() and SLESSetUpOnBlocks() to
224: enable more precise profiling of setting up the preconditioner.
225: These calls are optional, since both will be called within
226: SLESSolve() if they haven't been called already.
227: */
228: SLESSetUp(sles,b,x);
229: SLESSetUpOnBlocks(sles);
230: PetscGetTime(&tsetup2);
231: tsetup = tsetup2 - tsetup1;
233: /*
234: Tests "diagonal-scaling of preconditioned residual norm" as used
235: by many ODE integrator codes including PVODE. Note this is different
236: than diagonally scaling the matrix before computing the preconditioner
237: */
238: PetscOptionsHasName(PETSC_NULL,"-diagonal_scale",&diagonalscale);
239: if (diagonalscale) {
240: PC pc;
241: int j,start,end,n;
242: Vec scale;
243:
244: SLESGetPC(sles,&pc);
245: VecGetSize(x,&n);
246: VecDuplicate(x,&scale);
247: VecGetOwnershipRange(scale,&start,&end);
248: for (j=start; j<end; j++) {
249: VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
250: }
251: VecAssemblyBegin(scale);
252: VecAssemblyEnd(scale);
253: PCDiagonalScaleSet(pc,scale);
254: VecDestroy(scale);
256: }
258: PetscOptionsHasName(PETSC_NULL, "-null_space", &hasNullSpace);
259: if (hasNullSpace == PETSC_TRUE) {
260: MatNullSpace nullSpace;
261: PC pc;
263: MatNullSpaceCreate(PETSC_COMM_WORLD, 1, 0, PETSC_NULL, &nullSpace);
264: MatNullSpaceTest(nullSpace, A);
265: SLESGetPC(sles,&pc);
266: PCNullSpaceAttach(pc, nullSpace);
267: MatNullSpaceDestroy(nullSpace);
268: }
270: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
271: Solve system
272: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274: /*
275: Begin profiling next stage
276: */
277: PreLoadStage("SLESSolve");
279: /*
280: Solve linear system; we also explicitly time this stage.
281: */
282: PetscGetTime(&tsolve1);
283: if (trans) {
284: SLESSolveTranspose(sles,b,x,&its);
285: } else {
286: int num_rhs=1;
287: PetscOptionsGetInt(PETSC_NULL,"-num_rhs",&num_rhs,PETSC_NULL);
288: while ( num_rhs-- ) {
289: SLESSolve(sles,b,x,&its);
290: }
291: }
292: PetscGetTime(&tsolve2);
293: tsolve = tsolve2 - tsolve1;
295: /*
296: Conclude profiling this stage
297: */
298: PreLoadStage("Cleanup");
300: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
301: Check error, print output, free data structures.
302: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
304: /*
305: Check error
306: */
307: if (trans) {
308: MatMultTranspose(A,x,u);
309: } else {
310: MatMult(A,x,u);
311: }
312: VecAXPY(&none,b,u);
313: VecNorm(u,NORM_2,&norm);
315: /*
316: Write output (optinally using table for solver details).
317: - PetscPrintf() handles output for multiprocessor jobs
318: by printing from only one processor in the communicator.
319: - SLESView() prints information about the linear solver.
320: */
321: if (table) {
322: char *matrixname,slesinfo[120];
323: PetscViewer viewer;
325: /*
326: Open a string viewer; then write info to it.
327: */
328: PetscViewerStringOpen(PETSC_COMM_WORLD,slesinfo,120,&viewer);
329: SLESView(sles,viewer);
330: PetscStrrchr(file[PreLoadIt],'/',&matrixname);
331: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3d %2.0e %2.1e %2.1e %2.1e %s n",
332: matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,slesinfo);
334: /*
335: Destroy the viewer
336: */
337: PetscViewerDestroy(viewer);
338: } else {
339: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3dn",its);
340: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %An",norm);
341: }
343: } /* while ( num_numfac-- ) */
345: /*
346: Free work space. All PETSc objects should be destroyed when they
347: are no longer needed.
348: */
349: MatDestroy(A); VecDestroy(b);
350: VecDestroy(u); VecDestroy(x);
351: SLESDestroy(sles);
352: PreLoadEnd();
353: /* -----------------------------------------------------------
354: End of linear solver loop
355: ----------------------------------------------------------- */
357: PetscFinalize();
358: return 0;
359: }