Actual source code: ex30.c
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves a linear system.\n\
3: It is copied and intended to move dirty codes from ksp/tutorials/ex10.c and simplify ex10.c.\n\
4: Input parameters include\n\
5: -f0 <input_file> : first file to load (small system)\n\
6: -f1 <input_file> : second file to load (larger system)\n\n\
7: -trans : solve transpose system instead\n\n";
8: /*
9: This code can be used to test PETSc interface to other packages.\n\
10: Examples of command line options: \n\
11: ex30 -f0 <datafile> -ksp_type preonly \n\
12: -help -ksp_view \n\
13: -num_numfac <num_numfac> -num_rhs <num_rhs> \n\
14: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type superlu or superlu_dist or mumps \n\
15: -ksp_type preonly -pc_type cholesky -pc_factor_mat_solver_type mumps \n\
16: mpiexec -n <np> ex30 -f0 <datafile> -ksp_type cg -pc_type asm -pc_asm_type basic -sub_pc_type icc -mat_type sbaij
18: ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_type superlu -mat_superlu_conditionnumber -ckerror -mat_superlu_diagpivotthresh 0
19: ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type hypre -pc_hypre_type boomeramg -ksp_type fgmres -ckError
20: ./ex30 -f0 $D/small -mat_sigma -3.999999999999999 -ksp_type fgmres -pc_type lu -pc_factor_mat_solver_type petsc -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5 -ckerror
21: \n\n";
22: */
23: /*T
24: Concepts: KSP solving a linear system
25: Processors: n
26: T*/
28: #include <petscksp.h>
30: int main(int argc,char **args)
31: {
32: KSP ksp;
33: Mat A,B;
34: Vec x,b,u,b2; /* approx solution, RHS, exact solution */
35: PetscViewer fd; /* viewer */
36: char file[4][PETSC_MAX_PATH_LEN]; /* input file name */
37: PetscBool table = PETSC_FALSE,flg,flgB=PETSC_FALSE,trans=PETSC_FALSE,partition=PETSC_FALSE,initialguess = PETSC_FALSE;
38: PetscBool outputSoln=PETSC_FALSE;
39: PetscInt its,num_numfac;
40: PetscReal rnorm,enorm;
41: PetscBool preload=PETSC_TRUE,diagonalscale,isSymmetric,ckrnorm=PETSC_TRUE,Test_MatDuplicate=PETSC_FALSE,ckerror=PETSC_FALSE;
42: PetscMPIInt rank;
43: PetscScalar sigma;
44: PetscInt m;
46: PetscInitialize(&argc,&args,(char*)0,help);
47: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
48: PetscOptionsGetBool(NULL,NULL,"-table",&table,NULL);
49: PetscOptionsGetBool(NULL,NULL,"-trans",&trans,NULL);
50: PetscOptionsGetBool(NULL,NULL,"-partition",&partition,NULL);
51: PetscOptionsGetBool(NULL,NULL,"-initialguess",&initialguess,NULL);
52: PetscOptionsGetBool(NULL,NULL,"-output_solution",&outputSoln,NULL);
53: PetscOptionsGetBool(NULL,NULL,"-ckrnorm",&ckrnorm,NULL);
54: PetscOptionsGetBool(NULL,NULL,"-ckerror",&ckerror,NULL);
56: /*
57: Determine files from which we read the two linear systems
58: (matrix and right-hand-side vector).
59: */
60: PetscOptionsGetString(NULL,NULL,"-f",file[0],sizeof(file[0]),&flg);
61: if (flg) {
62: PetscStrcpy(file[1],file[0]);
63: preload = PETSC_FALSE;
64: } else {
65: PetscOptionsGetString(NULL,NULL,"-f0",file[0],sizeof(file[0]),&flg);
67: PetscOptionsGetString(NULL,NULL,"-f1",file[1],sizeof(file[1]),&flg);
68: if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
69: }
71: /* -----------------------------------------------------------
72: Beginning of linear solver loop
73: ----------------------------------------------------------- */
74: /*
75: Loop through the linear solve 2 times.
76: - The intention here is to preload and solve a small system;
77: then load another (larger) system and solve it as well.
78: This process preloads the instructions with the smaller
79: system so that more accurate performance monitoring (via
80: -log_view) can be done with the larger one (that actually
81: is the system of interest).
82: */
83: PetscPreLoadBegin(preload,"Load system");
85: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
86: Load system
87: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
89: /*
90: Open binary file. Note that we use FILE_MODE_READ to indicate
91: reading from this file.
92: */
93: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[PetscPreLoadIt],FILE_MODE_READ,&fd);
95: /*
96: Load the matrix and vector; then destroy the viewer.
97: */
98: MatCreate(PETSC_COMM_WORLD,&A);
99: MatSetFromOptions(A);
100: MatLoad(A,fd);
102: flg = PETSC_FALSE;
103: PetscOptionsGetString(NULL,NULL,"-rhs",file[2],sizeof(file[2]),&flg);
104: if (flg) { /* rhs is stored in a separate file */
105: PetscViewerDestroy(&fd);
106: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
107: VecCreate(PETSC_COMM_WORLD,&b);
108: VecLoad(b,fd);
109: } else {
110: /* if file contains no RHS, then use a vector of all ones */
111: PetscInfo(0,"Using vector of ones for RHS\n");
112: MatGetLocalSize(A,&m,NULL);
113: VecCreate(PETSC_COMM_WORLD,&b);
114: VecSetSizes(b,m,PETSC_DECIDE);
115: VecSetFromOptions(b);
116: VecSet(b,1.0);
117: PetscObjectSetName((PetscObject)b, "Rhs vector");
118: }
119: PetscViewerDestroy(&fd);
121: /* Test MatDuplicate() */
122: if (Test_MatDuplicate) {
123: MatDuplicate(A,MAT_COPY_VALUES,&B);
124: MatEqual(A,B,&flg);
125: if (!flg) {
126: PetscPrintf(PETSC_COMM_WORLD," A != B \n");
127: }
128: MatDestroy(&B);
129: }
131: /* Add a shift to A */
132: PetscOptionsGetScalar(NULL,NULL,"-mat_sigma",&sigma,&flg);
133: if (flg) {
134: PetscOptionsGetString(NULL,NULL,"-fB",file[2],sizeof(file[2]),&flgB);
135: if (flgB) {
136: /* load B to get A = A + sigma*B */
137: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[2],FILE_MODE_READ,&fd);
138: MatCreate(PETSC_COMM_WORLD,&B);
139: MatSetOptionsPrefix(B,"B_");
140: MatLoad(B,fd);
141: PetscViewerDestroy(&fd);
142: MatAXPY(A,sigma,B,DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
143: } else {
144: MatShift(A,sigma);
145: }
146: }
148: /* Make A singular for testing zero-pivot of ilu factorization */
149: /* Example: ./ex30 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
150: flg = PETSC_FALSE;
151: PetscOptionsGetBool(NULL,NULL, "-test_zeropivot", &flg,NULL);
152: if (flg) {
153: PetscInt row,ncols;
154: const PetscInt *cols;
155: const PetscScalar *vals;
156: PetscBool flg1=PETSC_FALSE;
157: PetscScalar *zeros;
158: row = 0;
159: MatGetRow(A,row,&ncols,&cols,&vals);
160: PetscCalloc1(ncols+1,&zeros);
161: flg1 = PETSC_FALSE;
162: PetscOptionsGetBool(NULL,NULL, "-set_row_zero", &flg1,NULL);
163: if (flg1) { /* set entire row as zero */
164: MatSetValues(A,1,&row,ncols,cols,zeros,INSERT_VALUES);
165: } else { /* only set (row,row) entry as zero */
166: MatSetValues(A,1,&row,1,&row,zeros,INSERT_VALUES);
167: }
168: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
169: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
170: }
172: /* Check whether A is symmetric */
173: flg = PETSC_FALSE;
174: PetscOptionsGetBool(NULL,NULL, "-check_symmetry", &flg,NULL);
175: if (flg) {
176: Mat Atrans;
177: MatTranspose(A, MAT_INITIAL_MATRIX,&Atrans);
178: MatEqual(A, Atrans, &isSymmetric);
179: if (isSymmetric) {
180: PetscPrintf(PETSC_COMM_WORLD,"A is symmetric \n");
181: } else {
182: PetscPrintf(PETSC_COMM_WORLD,"A is non-symmetric \n");
183: }
184: MatDestroy(&Atrans);
185: }
187: VecDuplicate(b,&b2);
188: VecDuplicate(b,&x);
189: PetscObjectSetName((PetscObject)x, "Solution vector");
190: VecDuplicate(b,&u);
191: PetscObjectSetName((PetscObject)u, "True Solution vector");
192: VecSet(x,0.0);
194: if (ckerror) { /* Set true solution */
195: VecSet(u,1.0);
196: MatMult(A,u,b);
197: }
199: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
200: Setup solve for system
201: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
203: if (partition) {
204: MatPartitioning mpart;
205: IS mis,nis,is;
206: PetscInt *count;
207: PetscMPIInt size;
208: Mat BB;
209: MPI_Comm_size(PETSC_COMM_WORLD,&size);
210: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
211: PetscMalloc1(size,&count);
212: MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
213: MatPartitioningSetAdjacency(mpart, A);
214: /* MatPartitioningSetVertexWeights(mpart, weight); */
215: MatPartitioningSetFromOptions(mpart);
216: MatPartitioningApply(mpart, &mis);
217: MatPartitioningDestroy(&mpart);
218: ISPartitioningToNumbering(mis,&nis);
219: ISPartitioningCount(mis,size,count);
220: ISDestroy(&mis);
221: ISInvertPermutation(nis, count[rank], &is);
222: PetscFree(count);
223: ISDestroy(&nis);
224: ISSort(is);
225: MatCreateSubMatrix(A,is,is,MAT_INITIAL_MATRIX,&BB);
227: /* need to move the vector also */
228: ISDestroy(&is);
229: MatDestroy(&A);
230: A = BB;
231: }
233: /*
234: Create linear solver; set operators; set runtime options.
235: */
236: KSPCreate(PETSC_COMM_WORLD,&ksp);
237: KSPSetInitialGuessNonzero(ksp,initialguess);
238: num_numfac = 1;
239: PetscOptionsGetInt(NULL,NULL,"-num_numfac",&num_numfac,NULL);
240: while (num_numfac--) {
242: KSPSetOperators(ksp,A,A);
243: KSPSetFromOptions(ksp);
245: /*
246: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
247: enable more precise profiling of setting up the preconditioner.
248: These calls are optional, since both will be called within
249: KSPSolve() if they haven't been called already.
250: */
251: KSPSetUp(ksp);
252: KSPSetUpOnBlocks(ksp);
254: /*
255: Tests "diagonal-scaling of preconditioned residual norm" as used
256: by many ODE integrator codes including SUNDIALS. Note this is different
257: than diagonally scaling the matrix before computing the preconditioner
258: */
259: diagonalscale = PETSC_FALSE;
260: PetscOptionsGetBool(NULL,NULL,"-diagonal_scale",&diagonalscale,NULL);
261: if (diagonalscale) {
262: PC pc;
263: PetscInt j,start,end,n;
264: Vec scale;
266: KSPGetPC(ksp,&pc);
267: VecGetSize(x,&n);
268: VecDuplicate(x,&scale);
269: VecGetOwnershipRange(scale,&start,&end);
270: for (j=start; j<end; j++) {
271: VecSetValue(scale,j,((PetscReal)(j+1))/((PetscReal)n),INSERT_VALUES);
272: }
273: VecAssemblyBegin(scale);
274: VecAssemblyEnd(scale);
275: PCSetDiagonalScale(pc,scale);
276: VecDestroy(&scale);
277: }
279: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
280: Solve system
281: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
282: /*
283: Solve linear system;
284: */
285: if (trans) {
286: KSPSolveTranspose(ksp,b,x);
287: KSPGetIterationNumber(ksp,&its);
288: } else {
289: PetscInt num_rhs=1;
290: PetscOptionsGetInt(NULL,NULL,"-num_rhs",&num_rhs,NULL);
292: while (num_rhs--) {
293: KSPSolve(ksp,b,x);
294: }
295: KSPGetIterationNumber(ksp,&its);
296: if (ckrnorm) { /* Check residual for each rhs */
297: if (trans) {
298: MatMultTranspose(A,x,b2);
299: } else {
300: MatMult(A,x,b2);
301: }
302: VecAXPY(b2,-1.0,b);
303: VecNorm(b2,NORM_2,&rnorm);
304: PetscPrintf(PETSC_COMM_WORLD," Number of iterations = %3D\n",its);
305: PetscPrintf(PETSC_COMM_WORLD," Residual norm %g\n",(double)rnorm);
306: }
307: if (ckerror && !trans) { /* Check error for each rhs */
308: /* VecView(x,PETSC_VIEWER_STDOUT_WORLD); */
309: VecAXPY(u,-1.0,x);
310: VecNorm(u,NORM_2,&enorm);
311: PetscPrintf(PETSC_COMM_WORLD," Error norm %g\n",(double)enorm);
312: }
314: } /* while (num_rhs--) */
316: /*
317: Write output (optinally using table for solver details).
318: - PetscPrintf() handles output for multiprocessor jobs
319: by printing from only one processor in the communicator.
320: - KSPView() prints information about the linear solver.
321: */
322: if (table && ckrnorm) {
323: char *matrixname,kspinfo[120];
324: PetscViewer viewer;
326: /*
327: Open a string viewer; then write info to it.
328: */
329: PetscViewerStringOpen(PETSC_COMM_WORLD,kspinfo,sizeof(kspinfo),&viewer);
330: KSPView(ksp,viewer);
331: PetscStrrchr(file[PetscPreLoadIt],'/',&matrixname);
332: PetscPrintf(PETSC_COMM_WORLD,"%-8.8s %3D %2.0e %s \n", matrixname,its,rnorm,kspinfo);
334: /*
335: Destroy the viewer
336: */
337: PetscViewerDestroy(&viewer);
338: }
340: PetscOptionsGetString(NULL,NULL,"-solution",file[3],sizeof(file[3]),&flg);
341: if (flg) {
342: PetscViewer viewer;
343: Vec xstar;
345: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file[3],FILE_MODE_READ,&viewer);
346: VecCreate(PETSC_COMM_WORLD,&xstar);
347: VecLoad(xstar,viewer);
348: VecAXPY(xstar, -1.0, x);
349: VecNorm(xstar, NORM_2, &enorm);
350: PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)enorm);
351: VecDestroy(&xstar);
352: PetscViewerDestroy(&viewer);
353: }
354: if (outputSoln) {
355: PetscViewer viewer;
357: PetscViewerBinaryOpen(PETSC_COMM_WORLD,"solution.petsc",FILE_MODE_WRITE,&viewer);
358: VecView(x, viewer);
359: PetscViewerDestroy(&viewer);
360: }
362: flg = PETSC_FALSE;
363: PetscOptionsGetBool(NULL,NULL, "-ksp_reason", &flg,NULL);
364: if (flg) {
365: KSPConvergedReason reason;
366: KSPGetConvergedReason(ksp,&reason);
367: PetscPrintf(PETSC_COMM_WORLD,"KSPConvergedReason: %D\n", reason);
368: }
370: } /* while (num_numfac--) */
372: /*
373: Free work space. All PETSc objects should be destroyed when they
374: are no longer needed.
375: */
376: MatDestroy(&A)); PetscCall(VecDestroy(&b);
377: VecDestroy(&u)); PetscCall(VecDestroy(&x);
378: VecDestroy(&b2);
379: KSPDestroy(&ksp);
380: if (flgB) MatDestroy(&B);
381: PetscPreLoadEnd();
382: /* -----------------------------------------------------------
383: End of linear solver loop
384: ----------------------------------------------------------- */
386: PetscFinalize();
387: return 0;
388: }
390: /*TEST
392: test:
393: args: -f ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type ilu -pc_factor_mat_ordering_type natural -num_numfac 2 -pc_factor_reuse_fill
394: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
395: output_file: output/ex30.out
397: test:
398: TODO: Matrix row/column sizes are not compatible with block size
399: suffix: 2
400: args: -f ${DATAFILESPATH}/matrices/arco1 -mat_type baij -matload_block_size 3 -ksp_type preonly -pc_type ilu -pc_factor_mat_ordering_type natural -num_numfac 2
401: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
403: test:
404: suffix: shiftnz
405: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5
406: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
408: test:
409: suffix: shiftpd
410: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type POSITIVE_DEFINITE
411: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
413: test:
414: suffix: shift_cholesky_aij
415: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5
416: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
417: output_file: output/ex30_shiftnz.out
419: test:
420: suffix: shiftpd_2
421: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE
422: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
424: test:
425: suffix: shift_cholesky_sbaij
426: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type NONZERO -pc_factor_shift_amount 1.e-5 -mat_type sbaij
427: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
428: output_file: output/ex30_shiftnz.out
430: test:
431: suffix: shiftpd_2_sbaij
432: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE -mat_type sbaij
433: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
434: output_file: output/ex30_shiftpd_2.out
436: test:
437: suffix: shiftinblocks
438: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type INBLOCKS
439: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
441: test:
442: suffix: shiftinblocks2
443: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS
444: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
445: output_file: output/ex30_shiftinblocks.out
447: test:
448: suffix: shiftinblockssbaij
449: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS -mat_type sbaij
450: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
451: output_file: output/ex30_shiftinblocks.out
453: TEST*/