Actual source code: ex10.c
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: This version first preloads and solves a small system, then loads \n\
4: another (larger) system and solves it as well. This example illustrates\n\
5: preloading of instructions with the smaller system so that more accurate\n\
6: performance monitoring can be done with the larger one (that actually\n\
7: is the system of interest). See the 'Performance Hints' chapter of the\n\
8: users manual for a discussion of preloading. Input parameters include\n\
9: -f0 <input_file> : first file to load (small system)\n\
10: -f1 <input_file> : second file to load (larger system)\n\n\
11: -trans : solve transpose system instead\n\n";
12: /*
13: This code can be used to test PETSc interface to other packages.\n\
14: Examples of command line options: \n\
15: ex10 -f0 <datafile> -ksp_type preonly \n\
16: -help -ksp_view \n\
17: -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
18: -ksp_type preonly -pc_type lu -mat_type aijspooles/superlu/superlu_dist/aijmumps \n\
19: -ksp_type preonly -pc_type cholesky -mat_type sbaijspooles/dscpack/sbaijmumps \n\
20: -f0 <A> -fB <B> -mat_type sbaijmumps -ksp_type preonly -pc_type cholesky -test_inertia -mat_sigma <sigma> \n\
21: mpirun -np <np> ex10 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
22: \n\n";
23: */
24: /*T
25: Concepts: KSP^solving a linear system
26: Processors: n
27: T*/
29: /*
30: Include "petscksp.h" so that we can use KSP solvers. Note that this file
31: automatically includes:
32: petsc.h - base PETSc routines petscvec.h - vectors
33: petscsys.h - system routines petscmat.h - matrices
34: petscis.h - index sets petscksp.h - Krylov subspace methods
35: petscviewer.h - viewers petscpc.h - preconditioners
36: */
37: #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][PETSC_MAX_PATH_LEN]; /* input file name */
48: PetscTruth table,flg,flgB=PETSC_FALSE,trans=PETSC_FALSE,partition=PETSC_FALSE;
50: PetscInt its,num_numfac,m,n,M;
51: PetscReal norm;
52: PetscLogDouble tsetup,tsetup1,tsetup2,tsolve,tsolve1,tsolve2;
53: PetscTruth preload=PETSC_TRUE,diagonalscale,isSymmetric,cknorm=PETSC_FALSE,Test_MatDuplicate=PETSC_FALSE;
54: PetscMPIInt rank;
55: PetscScalar sigma;
57: PetscInitialize(&argc,&args,(char *)0,help);
58: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
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],PETSC_MAX_PATH_LEN-1,&flg);
68: if (flg) {
69: PetscStrcpy(file[1],file[0]);
70: preload = PETSC_FALSE;
71: } else {
72: PetscOptionsGetString(PETSC_NULL,"-f0",file[0],PETSC_MAX_PATH_LEN-1,&flg);
73: if (!flg) SETERRQ(1,"Must indicate binary file with the -f0 or -f option");
74: PetscOptionsGetString(PETSC_NULL,"-f1",file[1],PETSC_MAX_PATH_LEN-1,&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 FILE_MODE_READ to indicate
98: reading from this file.
99: */
100: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PreLoadIt],FILE_MODE_READ,&fd);
101:
102: /*
103: Load the matrix and vector; then destroy the viewer.
104: */
105: MatLoad(fd,MATAIJ,&A);
106:
107: if (!preload){
108: flg = PETSC_FALSE;
109: PetscOptionsGetString(PETSC_NULL,"-rhs",file[2],PETSC_MAX_PATH_LEN-1,&flg);
110: if (flg){ /* rhs is stored in a separate file */
111: PetscViewerDestroy(fd);
112: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
113: }
114: }
115: if (rank){
116: PetscExceptionTry1(VecLoad(fd,PETSC_NULL,&b),PETSC_ERR_FILE_UNEXPECTED);
117: } else {
118: PetscExceptionTry1(VecLoad(fd,PETSC_NULL,&b),PETSC_ERR_FILE_READ);
119: }
120: if (PetscExceptionCaught(ierr,PETSC_ERR_FILE_UNEXPECTED) || PetscExceptionCaught(ierr,PETSC_ERR_FILE_READ)) { /* if file contains no RHS, then use a vector of all ones */
121: PetscInt m;
122: PetscScalar one = 1.0;
123: PetscInfo(0,"Using vector of ones for RHS\n");
124: MatGetLocalSize(A,&m,PETSC_NULL);
125: VecCreate(PETSC_COMM_WORLD,&b);
126: VecSetSizes(b,m,PETSC_DECIDE);
127: VecSetFromOptions(b);
128: VecSet(b,one);
129: } else
130: PetscViewerDestroy(fd);
132: /* Test MatDuplicate() */
133: if (Test_MatDuplicate){
134: MatDuplicate(A,MAT_COPY_VALUES,&B);
135: MatEqual(A,B,&flg);
136: if (!flg){
137: PetscPrintf(PETSC_COMM_WORLD," A != B \n");
138: }
139: MatDestroy(B);
140: }
142: /* Add a shift to A */
143: PetscOptionsGetScalar(PETSC_NULL,"-mat_sigma",&sigma,&flg);
144: if (flg) {
145: PetscOptionsGetString(PETSC_NULL,"-fB",file[2],PETSC_MAX_PATH_LEN-1,&flgB);
146: if (flgB){
147: /* load B to get A = A + sigma*B */
148: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
149: MatLoad(fd,MATAIJ,&B);
150: PetscViewerDestroy(fd);
151: MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
152: } else {
153: MatShift(A,sigma);
154: }
155: }
157: /* Make A singular for testing zero-pivot of ilu factorization */
158: /* Example: ./ex10 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
159: PetscOptionsHasName(PETSC_NULL, "-test_zeropivot", &flg);
160: if (flg) {
161: PetscInt row,ncols;
162: const PetscInt *cols;
163: const PetscScalar *vals;
164: PetscTruth flg1=PETSC_FALSE;
165: PetscScalar *zeros;
166: row = 0;
167: MatGetRow(A,row,&ncols,&cols,&vals);
168: PetscMalloc(sizeof(PetscScalar)*(ncols+1),&zeros);
169: PetscMemzero(zeros,(ncols+1)*sizeof(PetscScalar));
170: PetscOptionsHasName(PETSC_NULL, "-set_row_zero", &flg1);
171: if (flg1){ /* set entire row as zero */
172: MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
173: } else { /* only set (row,row) entry as zero */
174: MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
175: }
176: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
177: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
178: }
180: /* Check whether A is symmetric */
181: PetscOptionsHasName(PETSC_NULL, "-check_symmetry", &flg);
182: if (flg) {
183: Mat Atrans;
184: MatTranspose(A, &Atrans);
185: MatEqual(A, Atrans, &isSymmetric);
186: if (isSymmetric) {
187: PetscPrintf(PETSC_COMM_WORLD,"A is symmetric \n");
188: } else {
189: PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric \n");
190: }
191: MatDestroy(Atrans);
192: }
194: /*
195: If the loaded matrix is larger than the vector (due to being padded
196: to match the block size of the system), then create a new padded vector.
197: */
198:
199: MatGetLocalSize(A,&m,&n);
200: if (m != n) {
201: SETERRQ2(PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%d, %d)", m, n);
202: }
203: MatGetSize(A,&M,PETSC_NULL);
204: VecGetSize(b,&m);
205: if (M != m) { /* Create a new vector b by padding the old one */
206: PetscInt j,mvec,start,end,indx;
207: Vec tmp;
208: PetscScalar *bold;
210: VecCreate(PETSC_COMM_WORLD,&tmp);
211: VecSetSizes(tmp,n,PETSC_DECIDE);
212: VecSetFromOptions(tmp);
213: VecGetOwnershipRange(b,&start,&end);
214: VecGetLocalSize(b,&mvec);
215: VecGetArray(b,&bold);
216: for (j=0; j<mvec; j++) {
217: indx = start+j;
218: VecSetValues(tmp,1,&indx,bold+j,INSERT_VALUES);
219: }
220: VecRestoreArray(b,&bold);
221: VecDestroy(b);
222: VecAssemblyBegin(tmp);
223: VecAssemblyEnd(tmp);
224: b = tmp;
225: }
226: VecDuplicate(b,&x);
227: VecDuplicate(b,&u);
228: VecSet(x,0.0);
230: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
231: Setup solve for system
232: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
235: if (partition) {
236: MatPartitioning mpart;
237: IS mis,nis,isn,is;
238: PetscInt *count;
239: PetscMPIInt size;
240: Mat BB;
241: MPI_Comm_size(PETSC_COMM_WORLD,&size);
242: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
243: PetscMalloc(size*sizeof(PetscInt),&count);
244: MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
245: MatPartitioningSetAdjacency(mpart, A);
246: /* MatPartitioningSetVertexWeights(mpart, weight); */
247: MatPartitioningSetFromOptions(mpart);
248: MatPartitioningApply(mpart, &mis);
249: MatPartitioningDestroy(mpart);
250: ISPartitioningToNumbering(mis,&nis);
251: ISPartitioningCount(mis,count);
252: ISDestroy(mis);
253: ISInvertPermutation(nis, count[rank], &is);
254: PetscFree(count);
255: ISDestroy(nis);
256: ISSort(is);
257: ISAllGather(is,&isn);
258: MatGetSubMatrix(A,is,isn,PETSC_DECIDE,MAT_INITIAL_MATRIX,&BB);
260: /* need to move the vector also */
261: ISDestroy(is);
262: ISDestroy(isn);
263: MatDestroy(A);
264: A = BB;
265: }
266:
267: /*
268: Conclude profiling last stage; begin profiling next stage.
269: */
270: PreLoadStage("KSPSetUp");
272: /*
273: We also explicitly time this stage via PetscGetTime()
274: */
275: PetscGetTime(&tsetup1);
277: /*
278: Create linear solver; set operators; set runtime options.
279: */
280: KSPCreate(PETSC_COMM_WORLD,&ksp);
281: KSPSetType(ksp,KSPLGMRES);
282: num_numfac = 1;
283: PetscOptionsGetInt(PETSC_NULL,"-num_numfac",&num_numfac,PETSC_NULL);
284: while ( num_numfac-- ){
285:
286: KSPSetOperators(ksp,A,A,SAME_NONZERO_PATTERN);
287: KSPSetFromOptions(ksp);
289: /*
290: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
291: enable more precise profiling of setting up the preconditioner.
292: These calls are optional, since both will be called within
293: KSPSolve() if they haven't been called already.
294: */
295: KSPSetUp(ksp);
296: KSPSetUpOnBlocks(ksp);
297: PetscGetTime(&tsetup2);
298: tsetup = tsetup2 - tsetup1;
300: /*
301: Test MatGetInertia()
302: Usage:
303: ex10 -f0 <mat_binaryfile> -ksp_type preonly -pc_type cholesky -mat_type seqsbaij -test_inertia -mat_sigma <sigma>
304: */
305: PetscOptionsHasName(PETSC_NULL,"-test_inertia",&flg);
306: if (flg){
307: PC pc;
308: PetscInt nneg, nzero, npos;
309: Mat F;
310:
311: KSPGetPC(ksp,&pc);
312: PCGetFactoredMatrix(pc,&F);
313: MatGetInertia(F,&nneg,&nzero,&npos);
314: PetscPrintf(PETSC_COMM_SELF," MatInertia: nneg: %D, nzero: %D, npos: %D\n",nneg,nzero,npos);
315: }
317: /*
318: Tests "diagonal-scaling of preconditioned residual norm" as used
319: by many ODE integrator codes including SUNDIALS. Note this is different
320: than diagonally scaling the matrix before computing the preconditioner
321: */
322: PetscOptionsHasName(PETSC_NULL,"-diagonal_scale",&diagonalscale);
323: if (diagonalscale) {
324: PC pc;
325: PetscInt j,start,end,n;
326: Vec scale;
327:
328: KSPGetPC(ksp,&pc);
329: VecGetSize(x,&n);
330: VecDuplicate(x,&scale);
331: VecGetOwnershipRange(scale,&start,&end);
332: for (j=start; j<end; j++) {
333: VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
334: }
335: VecAssemblyBegin(scale);
336: VecAssemblyEnd(scale);
337: PCDiagonalScaleSet(pc,scale);
338: VecDestroy(scale);
339: }
341: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
342: Solve system
343: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
345: /*
346: Begin profiling next stage
347: */
348: PreLoadStage("KSPSolve");
350: /*
351: Solve linear system; we also explicitly time this stage.
352: */
353: PetscGetTime(&tsolve1);
354: if (trans) {
355: KSPSolveTranspose(ksp,b,x);
356: KSPGetIterationNumber(ksp,&its);
357: } else {
358: PetscInt num_rhs=1;
359: PetscOptionsGetInt(PETSC_NULL,"-num_rhs",&num_rhs,PETSC_NULL);
360: PetscOptionsHasName(PETSC_NULL,"-cknorm",&cknorm);
361: while ( num_rhs-- ) {
362: KSPSolve(ksp,b,x);
363: }
364: KSPGetIterationNumber(ksp,&its);
365: if (cknorm){ /* Check error for each rhs */
366: if (trans) {
367: MatMultTranspose(A,x,u);
368: } else {
369: MatMult(A,x,u);
370: }
371: VecAXPY(u,-1.0,b);
372: VecNorm(u,NORM_2,&norm);
373: PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);
374: PetscPrintf(PETSC_COMM_WORLD," Residual norm %A\n",norm);
375: }
376: } /* while ( num_rhs-- ) */
377: PetscGetTime(&tsolve2);
378: tsolve = tsolve2 - tsolve1;
380: /*
381: Conclude profiling this stage
382: */
383: PreLoadStage("Cleanup");
385: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
386: Check error, print output, free data structures.
387: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
389: /*
390: Check error
391: */
392: if (trans) {
393: MatMultTranspose(A,x,u);
394: } else {
395: MatMult(A,x,u);
396: }
397: VecAXPY(u,-1.0,b);
398: VecNorm(u,NORM_2,&norm);
400: /*
401: Write output (optinally using table for solver details).
402: - PetscPrintf() handles output for multiprocessor jobs
403: by printing from only one processor in the communicator.
404: - KSPView() prints information about the linear solver.
405: */
406: if (table) {
407: char *matrixname,kspinfo[120];
408: PetscViewer viewer;
410: /*
411: Open a string viewer; then write info to it.
412: */
413: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,120,&viewer);
414: KSPView(ksp,viewer);
415: PetscStrrchr(file[PreLoadIt],'/',&matrixname);
416: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %2.1e %2.1e %2.1e %s \n",
417: matrixname,its,norm,tsetup+tsolve,tsetup,tsolve,kspinfo);
419: /*
420: Destroy the viewer
421: */
422: PetscViewerDestroy(viewer);
423: } else {
424: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
425: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %A\n",norm);
426: }
428: PetscOptionsHasName(PETSC_NULL, "-ksp_reason", &flg);
429: if (flg){
430: KSPConvergedReason reason;
431: KSPGetConvergedReason(ksp,&reason);
432: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
433: }
434:
435: } /* while ( num_numfac-- ) */
437: /*
438: Free work space. All PETSc objects should be destroyed when they
439: are no longer needed.
440: */
441: MatDestroy(A); VecDestroy(b);
442: VecDestroy(u); VecDestroy(x);
443: KSPDestroy(ksp);
444: if (flgB) { MatDestroy(B); }
445: PreLoadEnd();
446: /* -----------------------------------------------------------
447: End of linear solver loop
448: ----------------------------------------------------------- */
450: PetscFinalize();
451: return 0;
452: }