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: */
24: #include <petscksp.h>
26: int main(int argc, char **args)
27: {
28: KSP ksp;
29: Mat A, B;
30: Vec x, b, u, b2; /* approx solution, RHS, exact solution */
31: PetscViewer fd; /* viewer */
32: char file[4][PETSC_MAX_PATH_LEN]; /* input file name */
33: PetscBool table = PETSC_FALSE, flg, flgB = PETSC_FALSE, trans = PETSC_FALSE, partition = PETSC_FALSE, initialguess = PETSC_FALSE;
34: PetscBool outputSoln = PETSC_FALSE;
35: PetscInt its, num_numfac;
36: PetscReal rnorm, enorm;
37: PetscBool preload = PETSC_TRUE, diagonalscale, isSymmetric, ckrnorm = PETSC_TRUE, Test_MatDuplicate = PETSC_FALSE, ckerror = PETSC_FALSE;
38: PetscMPIInt rank;
39: PetscScalar sigma;
40: PetscInt m;
43: PetscInitialize(&argc, &args, (char *)0, help);
44: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
45: PetscOptionsGetBool(NULL, NULL, "-table", &table, NULL);
46: PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL);
47: PetscOptionsGetBool(NULL, NULL, "-partition", &partition, NULL);
48: PetscOptionsGetBool(NULL, NULL, "-initialguess", &initialguess, NULL);
49: PetscOptionsGetBool(NULL, NULL, "-output_solution", &outputSoln, NULL);
50: PetscOptionsGetBool(NULL, NULL, "-ckrnorm", &ckrnorm, NULL);
51: PetscOptionsGetBool(NULL, NULL, "-ckerror", &ckerror, NULL);
53: /*
54: Determine files from which we read the two linear systems
55: (matrix and right-hand-side vector).
56: */
57: PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg);
58: if (flg) {
59: PetscStrcpy(file[1], file[0]);
60: preload = PETSC_FALSE;
61: } else {
62: PetscOptionsGetString(NULL, NULL, "-f0", file[0], sizeof(file[0]), &flg);
64: PetscOptionsGetString(NULL, NULL, "-f1", file[1], sizeof(file[1]), &flg);
65: if (!flg) preload = PETSC_FALSE; /* don't bother with second system */
66: }
68: /* -----------------------------------------------------------
69: Beginning of linear solver loop
70: ----------------------------------------------------------- */
71: /*
72: Loop through the linear solve 2 times.
73: - The intention here is to preload and solve a small system;
74: then load another (larger) system and solve it as well.
75: This process preloads the instructions with the smaller
76: system so that more accurate performance monitoring (via
77: -log_view) can be done with the larger one (that actually
78: is the system of interest).
79: */
80: PetscPreLoadBegin(preload, "Load system");
82: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
83: Load system
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: /*
87: Open binary file. Note that we use FILE_MODE_READ to indicate
88: reading from this file.
89: */
90: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[PetscPreLoadIt], FILE_MODE_READ, &fd);
92: /*
93: Load the matrix and vector; then destroy the viewer.
94: */
95: MatCreate(PETSC_COMM_WORLD, &A);
96: MatSetFromOptions(A);
97: MatLoad(A, fd);
99: flg = PETSC_FALSE;
100: PetscOptionsGetString(NULL, NULL, "-rhs", file[2], sizeof(file[2]), &flg);
101: if (flg) { /* rhs is stored in a separate file */
102: PetscViewerDestroy(&fd);
103: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2], FILE_MODE_READ, &fd);
104: VecCreate(PETSC_COMM_WORLD, &b);
105: VecLoad(b, fd);
106: } else {
107: /* if file contains no RHS, then use a vector of all ones */
108: PetscInfo(0, "Using vector of ones for RHS\n");
109: MatGetLocalSize(A, &m, NULL);
110: VecCreate(PETSC_COMM_WORLD, &b);
111: VecSetSizes(b, m, PETSC_DECIDE);
112: VecSetFromOptions(b);
113: VecSet(b, 1.0);
114: PetscObjectSetName((PetscObject)b, "Rhs vector");
115: }
116: PetscViewerDestroy(&fd);
118: /* Test MatDuplicate() */
119: if (Test_MatDuplicate) {
120: MatDuplicate(A, MAT_COPY_VALUES, &B);
121: MatEqual(A, B, &flg);
122: if (!flg) PetscPrintf(PETSC_COMM_WORLD, " A != B \n");
123: MatDestroy(&B);
124: }
126: /* Add a shift to A */
127: PetscOptionsGetScalar(NULL, NULL, "-mat_sigma", &sigma, &flg);
128: if (flg) {
129: PetscOptionsGetString(NULL, NULL, "-fB", file[2], sizeof(file[2]), &flgB);
130: if (flgB) {
131: /* load B to get A = A + sigma*B */
132: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2], FILE_MODE_READ, &fd);
133: MatCreate(PETSC_COMM_WORLD, &B);
134: MatSetOptionsPrefix(B, "B_");
135: MatLoad(B, fd);
136: PetscViewerDestroy(&fd);
137: MatAXPY(A, sigma, B, DIFFERENT_NONZERO_PATTERN); /* A <- sigma*B + A */
138: } else {
139: MatShift(A, sigma);
140: }
141: }
143: /* Make A singular for testing zero-pivot of ilu factorization */
144: /* Example: ./ex30 -f0 <datafile> -test_zeropivot -set_row_zero -pc_factor_shift_nonzero */
145: flg = PETSC_FALSE;
146: PetscOptionsGetBool(NULL, NULL, "-test_zeropivot", &flg, NULL);
147: if (flg) {
148: PetscInt row, ncols;
149: const PetscInt *cols;
150: const PetscScalar *vals;
151: PetscBool flg1 = PETSC_FALSE;
152: PetscScalar *zeros;
153: row = 0;
154: MatGetRow(A, row, &ncols, &cols, &vals);
155: PetscCalloc1(ncols + 1, &zeros);
156: flg1 = PETSC_FALSE;
157: PetscOptionsGetBool(NULL, NULL, "-set_row_zero", &flg1, NULL);
158: if (flg1) { /* set entire row as zero */
159: MatSetValues(A, 1, &row, ncols, cols, zeros, INSERT_VALUES);
160: } else { /* only set (row,row) entry as zero */
161: MatSetValues(A, 1, &row, 1, &row, zeros, INSERT_VALUES);
162: }
163: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
164: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
165: }
167: /* Check whether A is symmetric */
168: flg = PETSC_FALSE;
169: PetscOptionsGetBool(NULL, NULL, "-check_symmetry", &flg, NULL);
170: if (flg) {
171: Mat Atrans;
172: MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans);
173: MatEqual(A, Atrans, &isSymmetric);
174: if (isSymmetric) {
175: PetscPrintf(PETSC_COMM_WORLD, "A is symmetric \n");
176: } else {
177: PetscPrintf(PETSC_COMM_WORLD, "A is non-symmetric \n");
178: }
179: MatDestroy(&Atrans);
180: }
182: VecDuplicate(b, &b2);
183: VecDuplicate(b, &x);
184: PetscObjectSetName((PetscObject)x, "Solution vector");
185: VecDuplicate(b, &u);
186: PetscObjectSetName((PetscObject)u, "True Solution vector");
187: VecSet(x, 0.0);
189: if (ckerror) { /* Set true solution */
190: VecSet(u, 1.0);
191: MatMult(A, u, b);
192: }
194: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
195: Setup solve for system
196: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: if (partition) {
199: MatPartitioning mpart;
200: IS mis, nis, is;
201: PetscInt *count;
202: PetscMPIInt size;
203: Mat BB;
204: MPI_Comm_size(PETSC_COMM_WORLD, &size);
205: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
206: PetscMalloc1(size, &count);
207: MatPartitioningCreate(PETSC_COMM_WORLD, &mpart);
208: MatPartitioningSetAdjacency(mpart, A);
209: /* MatPartitioningSetVertexWeights(mpart, weight); */
210: MatPartitioningSetFromOptions(mpart);
211: MatPartitioningApply(mpart, &mis);
212: MatPartitioningDestroy(&mpart);
213: ISPartitioningToNumbering(mis, &nis);
214: ISPartitioningCount(mis, size, count);
215: ISDestroy(&mis);
216: ISInvertPermutation(nis, count[rank], &is);
217: PetscFree(count);
218: ISDestroy(&nis);
219: ISSort(is);
220: MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &BB);
222: /* need to move the vector also */
223: ISDestroy(&is);
224: MatDestroy(&A);
225: A = BB;
226: }
228: /*
229: Create linear solver; set operators; set runtime options.
230: */
231: KSPCreate(PETSC_COMM_WORLD, &ksp);
232: KSPSetInitialGuessNonzero(ksp, initialguess);
233: num_numfac = 1;
234: PetscOptionsGetInt(NULL, NULL, "-num_numfac", &num_numfac, NULL);
235: while (num_numfac--) {
236: KSPSetOperators(ksp, A, A);
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: KSPSetUp(ksp);
246: KSPSetUpOnBlocks(ksp);
248: /*
249: Tests "diagonal-scaling of preconditioned residual norm" as used
250: by many ODE integrator codes including SUNDIALS. Note this is different
251: than diagonally scaling the matrix before computing the preconditioner
252: */
253: diagonalscale = PETSC_FALSE;
254: PetscOptionsGetBool(NULL, NULL, "-diagonal_scale", &diagonalscale, NULL);
255: if (diagonalscale) {
256: PC pc;
257: PetscInt j, start, end, n;
258: Vec scale;
260: KSPGetPC(ksp, &pc);
261: VecGetSize(x, &n);
262: VecDuplicate(x, &scale);
263: VecGetOwnershipRange(scale, &start, &end);
264: for (j = start; j < end; j++) VecSetValue(scale, j, ((PetscReal)(j + 1)) / ((PetscReal)n), INSERT_VALUES);
265: VecAssemblyBegin(scale);
266: VecAssemblyEnd(scale);
267: PCSetDiagonalScale(pc, scale);
268: VecDestroy(&scale);
269: }
271: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
272: Solve system
273: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274: /*
275: Solve linear system;
276: */
277: if (trans) {
278: KSPSolveTranspose(ksp, b, x);
279: KSPGetIterationNumber(ksp, &its);
280: } else {
281: PetscInt num_rhs = 1;
282: PetscOptionsGetInt(NULL, NULL, "-num_rhs", &num_rhs, NULL);
284: while (num_rhs--) KSPSolve(ksp, b, x);
285: KSPGetIterationNumber(ksp, &its);
286: if (ckrnorm) { /* Check residual for each rhs */
287: if (trans) {
288: MatMultTranspose(A, x, b2);
289: } else {
290: MatMult(A, x, b2);
291: }
292: VecAXPY(b2, -1.0, b);
293: VecNorm(b2, NORM_2, &rnorm);
294: PetscPrintf(PETSC_COMM_WORLD, " Number of iterations = %3" PetscInt_FMT "\n", its);
295: PetscPrintf(PETSC_COMM_WORLD, " Residual norm %g\n", (double)rnorm);
296: }
297: if (ckerror && !trans) { /* Check error for each rhs */
298: /* VecView(x,PETSC_VIEWER_STDOUT_WORLD); */
299: VecAXPY(u, -1.0, x);
300: VecNorm(u, NORM_2, &enorm);
301: PetscPrintf(PETSC_COMM_WORLD, " Error norm %g\n", (double)enorm);
302: }
304: } /* while (num_rhs--) */
306: /*
307: Write output (optionally using table for solver details).
308: - PetscPrintf() handles output for multiprocessor jobs
309: by printing from only one processor in the communicator.
310: - KSPView() prints information about the linear solver.
311: */
312: if (table && ckrnorm) {
313: char *matrixname, kspinfo[120];
314: PetscViewer viewer;
316: /*
317: Open a string viewer; then write info to it.
318: */
319: PetscViewerStringOpen(PETSC_COMM_WORLD, kspinfo, sizeof(kspinfo), &viewer);
320: KSPView(ksp, viewer);
321: PetscStrrchr(file[PetscPreLoadIt], '/', &matrixname);
322: PetscPrintf(PETSC_COMM_WORLD, "%-8.8s %3" PetscInt_FMT " %2.0e %s \n", matrixname, its, (double)rnorm, kspinfo);
324: /*
325: Destroy the viewer
326: */
327: PetscViewerDestroy(&viewer);
328: }
330: PetscOptionsGetString(NULL, NULL, "-solution", file[3], sizeof(file[3]), &flg);
331: if (flg) {
332: PetscViewer viewer;
333: Vec xstar;
335: PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[3], FILE_MODE_READ, &viewer);
336: VecCreate(PETSC_COMM_WORLD, &xstar);
337: VecLoad(xstar, viewer);
338: VecAXPY(xstar, -1.0, x);
339: VecNorm(xstar, NORM_2, &enorm);
340: PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)enorm);
341: VecDestroy(&xstar);
342: PetscViewerDestroy(&viewer);
343: }
344: if (outputSoln) {
345: PetscViewer viewer;
347: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "solution.petsc", FILE_MODE_WRITE, &viewer);
348: VecView(x, viewer);
349: PetscViewerDestroy(&viewer);
350: }
352: flg = PETSC_FALSE;
353: PetscOptionsGetBool(NULL, NULL, "-ksp_reason", &flg, NULL);
354: if (flg) {
355: KSPConvergedReason reason;
356: KSPGetConvergedReason(ksp, &reason);
357: PetscPrintf(PETSC_COMM_WORLD, "KSPConvergedReason: %s\n", KSPConvergedReasons[reason]);
358: }
360: } /* while (num_numfac--) */
362: /*
363: Free work space. All PETSc objects should be destroyed when they
364: are no longer needed.
365: */
366: MatDestroy(&A);
367: VecDestroy(&b);
368: VecDestroy(&u);
369: VecDestroy(&x);
370: VecDestroy(&b2);
371: KSPDestroy(&ksp);
372: if (flgB) MatDestroy(&B);
373: PetscPreLoadEnd();
374: /* -----------------------------------------------------------
375: End of linear solver loop
376: ----------------------------------------------------------- */
378: PetscFinalize();
379: return 0;
380: }
382: /*TEST
384: test:
385: args: -f ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type ilu -pc_factor_mat_ordering_type natural -num_numfac 2 -pc_factor_reuse_fill
386: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
387: output_file: output/ex30.out
389: test:
390: TODO: Matrix row/column sizes are not compatible with block size
391: suffix: 2
392: 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
393: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
395: test:
396: suffix: shiftnz
397: 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
398: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
400: test:
401: suffix: shiftpd
402: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type POSITIVE_DEFINITE
403: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
405: test:
406: suffix: shift_cholesky_aij
407: 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
408: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
409: output_file: output/ex30_shiftnz.out
411: test:
412: suffix: shiftpd_2
413: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE
414: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
416: test:
417: suffix: shift_cholesky_sbaij
418: 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
419: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
420: output_file: output/ex30_shiftnz.out
422: test:
423: suffix: shiftpd_2_sbaij
424: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE -mat_type sbaij
425: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
426: output_file: output/ex30_shiftpd_2.out
428: test:
429: suffix: shiftinblocks
430: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type INBLOCKS
431: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
433: test:
434: suffix: shiftinblocks2
435: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS
436: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
437: output_file: output/ex30_shiftinblocks.out
439: test:
440: suffix: shiftinblockssbaij
441: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS -mat_type sbaij
442: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
443: output_file: output/ex30_shiftinblocks.out
445: TEST*/