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