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;
42: PetscFunctionBeginUser;
43: PetscCall(PetscInitialize(&argc, &args, (char *)0, help));
44: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
45: PetscCall(PetscOptionsGetBool(NULL, NULL, "-table", &table, NULL));
46: PetscCall(PetscOptionsGetBool(NULL, NULL, "-trans", &trans, NULL));
47: PetscCall(PetscOptionsGetBool(NULL, NULL, "-partition", &partition, NULL));
48: PetscCall(PetscOptionsGetBool(NULL, NULL, "-initialguess", &initialguess, NULL));
49: PetscCall(PetscOptionsGetBool(NULL, NULL, "-output_solution", &outputSoln, NULL));
50: PetscCall(PetscOptionsGetBool(NULL, NULL, "-ckrnorm", &ckrnorm, NULL));
51: PetscCall(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: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file[0], sizeof(file[0]), &flg));
58: if (flg) {
59: PetscCall(PetscStrncpy(file[1], file[0], sizeof(file[1])));
60: preload = PETSC_FALSE;
61: } else {
62: PetscCall(PetscOptionsGetString(NULL, NULL, "-f0", file[0], sizeof(file[0]), &flg));
63: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER_INPUT, "Must indicate binary file with the -f0 or -f option");
64: PetscCall(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: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[PetscPreLoadIt], FILE_MODE_READ, &fd));
92: /*
93: Load the matrix and vector; then destroy the viewer.
94: */
95: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
96: PetscCall(MatSetFromOptions(A));
97: PetscCall(MatLoad(A, fd));
99: flg = PETSC_FALSE;
100: PetscCall(PetscOptionsGetString(NULL, NULL, "-rhs", file[2], sizeof(file[2]), &flg));
101: if (flg) { /* rhs is stored in a separate file */
102: PetscCall(PetscViewerDestroy(&fd));
103: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2], FILE_MODE_READ, &fd));
104: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
105: PetscCall(VecLoad(b, fd));
106: } else {
107: /* if file contains no RHS, then use a vector of all ones */
108: PetscCall(PetscInfo(0, "Using vector of ones for RHS\n"));
109: PetscCall(MatGetLocalSize(A, &m, NULL));
110: PetscCall(VecCreate(PETSC_COMM_WORLD, &b));
111: PetscCall(VecSetSizes(b, m, PETSC_DECIDE));
112: PetscCall(VecSetFromOptions(b));
113: PetscCall(VecSet(b, 1.0));
114: PetscCall(PetscObjectSetName((PetscObject)b, "Rhs vector"));
115: }
116: PetscCall(PetscViewerDestroy(&fd));
118: /* Test MatDuplicate() */
119: if (Test_MatDuplicate) {
120: PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &B));
121: PetscCall(MatEqual(A, B, &flg));
122: if (!flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD, " A != B \n"));
123: PetscCall(MatDestroy(&B));
124: }
126: /* Add a shift to A */
127: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-mat_sigma", &sigma, &flg));
128: if (flg) {
129: PetscCall(PetscOptionsGetString(NULL, NULL, "-fB", file[2], sizeof(file[2]), &flgB));
130: if (flgB) {
131: /* load B to get A = A + sigma*B */
132: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[2], FILE_MODE_READ, &fd));
133: PetscCall(MatCreate(PETSC_COMM_WORLD, &B));
134: PetscCall(MatSetOptionsPrefix(B, "B_"));
135: PetscCall(MatLoad(B, fd));
136: PetscCall(PetscViewerDestroy(&fd));
137: PetscCall(MatAXPY(A, sigma, B, DIFFERENT_NONZERO_PATTERN)); /* A <- sigma*B + A */
138: } else {
139: PetscCall(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: PetscCall(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: PetscCall(MatGetRow(A, row, &ncols, &cols, &vals));
155: PetscCall(PetscCalloc1(ncols + 1, &zeros));
156: flg1 = PETSC_FALSE;
157: PetscCall(PetscOptionsGetBool(NULL, NULL, "-set_row_zero", &flg1, NULL));
158: if (flg1) { /* set entire row as zero */
159: PetscCall(MatSetValues(A, 1, &row, ncols, cols, zeros, INSERT_VALUES));
160: } else { /* only set (row,row) entry as zero */
161: PetscCall(MatSetValues(A, 1, &row, 1, &row, zeros, INSERT_VALUES));
162: }
163: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
164: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
165: }
167: /* Check whether A is symmetric */
168: flg = PETSC_FALSE;
169: PetscCall(PetscOptionsGetBool(NULL, NULL, "-check_symmetry", &flg, NULL));
170: if (flg) {
171: Mat Atrans;
172: PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atrans));
173: PetscCall(MatEqual(A, Atrans, &isSymmetric));
174: if (isSymmetric) {
175: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "A is symmetric \n"));
176: } else {
177: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "A is non-symmetric \n"));
178: }
179: PetscCall(MatDestroy(&Atrans));
180: }
182: PetscCall(VecDuplicate(b, &b2));
183: PetscCall(VecDuplicate(b, &x));
184: PetscCall(PetscObjectSetName((PetscObject)x, "Solution vector"));
185: PetscCall(VecDuplicate(b, &u));
186: PetscCall(PetscObjectSetName((PetscObject)u, "True Solution vector"));
187: PetscCall(VecSet(x, 0.0));
189: if (ckerror) { /* Set true solution */
190: PetscCall(VecSet(u, 1.0));
191: PetscCall(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: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
205: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
206: PetscCall(PetscMalloc1(size, &count));
207: PetscCall(MatPartitioningCreate(PETSC_COMM_WORLD, &mpart));
208: PetscCall(MatPartitioningSetAdjacency(mpart, A));
209: /* PetscCall(MatPartitioningSetVertexWeights(mpart, weight)); */
210: PetscCall(MatPartitioningSetFromOptions(mpart));
211: PetscCall(MatPartitioningApply(mpart, &mis));
212: PetscCall(MatPartitioningDestroy(&mpart));
213: PetscCall(ISPartitioningToNumbering(mis, &nis));
214: PetscCall(ISPartitioningCount(mis, size, count));
215: PetscCall(ISDestroy(&mis));
216: PetscCall(ISInvertPermutation(nis, count[rank], &is));
217: PetscCall(PetscFree(count));
218: PetscCall(ISDestroy(&nis));
219: PetscCall(ISSort(is));
220: PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &BB));
222: /* need to move the vector also */
223: PetscCall(ISDestroy(&is));
224: PetscCall(MatDestroy(&A));
225: A = BB;
226: }
228: /*
229: Create linear solver; set operators; set runtime options.
230: */
231: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
232: PetscCall(KSPSetInitialGuessNonzero(ksp, initialguess));
233: num_numfac = 1;
234: PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_numfac", &num_numfac, NULL));
235: while (num_numfac--) {
236: PetscCall(KSPSetOperators(ksp, A, A));
237: PetscCall(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: PetscCall(KSPSetUp(ksp));
246: PetscCall(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: PetscCall(PetscOptionsGetBool(NULL, NULL, "-diagonal_scale", &diagonalscale, NULL));
255: if (diagonalscale) {
256: PC pc;
257: PetscInt j, start, end, n;
258: Vec scale;
260: PetscCall(KSPGetPC(ksp, &pc));
261: PetscCall(VecGetSize(x, &n));
262: PetscCall(VecDuplicate(x, &scale));
263: PetscCall(VecGetOwnershipRange(scale, &start, &end));
264: for (j = start; j < end; j++) PetscCall(VecSetValue(scale, j, ((PetscReal)(j + 1)) / ((PetscReal)n), INSERT_VALUES));
265: PetscCall(VecAssemblyBegin(scale));
266: PetscCall(VecAssemblyEnd(scale));
267: PetscCall(PCSetDiagonalScale(pc, scale));
268: PetscCall(VecDestroy(&scale));
269: }
271: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
272: Solve system
273: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274: /*
275: Solve linear system;
276: */
277: if (trans) {
278: PetscCall(KSPSolveTranspose(ksp, b, x));
279: PetscCall(KSPGetIterationNumber(ksp, &its));
280: } else {
281: PetscInt num_rhs = 1;
282: PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_rhs", &num_rhs, NULL));
284: while (num_rhs--) PetscCall(KSPSolve(ksp, b, x));
285: PetscCall(KSPGetIterationNumber(ksp, &its));
286: if (ckrnorm) { /* Check residual for each rhs */
287: if (trans) {
288: PetscCall(MatMultTranspose(A, x, b2));
289: } else {
290: PetscCall(MatMult(A, x, b2));
291: }
292: PetscCall(VecAXPY(b2, -1.0, b));
293: PetscCall(VecNorm(b2, NORM_2, &rnorm));
294: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Number of iterations = %3" PetscInt_FMT "\n", its));
295: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Residual norm %g\n", (double)rnorm));
296: }
297: if (ckerror && !trans) { /* Check error for each rhs */
298: /* PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); */
299: PetscCall(VecAXPY(u, -1.0, x));
300: PetscCall(VecNorm(u, NORM_2, &enorm));
301: PetscCall(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 = NULL, kspinfo[120];
314: PetscViewer viewer;
316: /*
317: Open a string viewer; then write info to it.
318: */
319: PetscCall(PetscViewerStringOpen(PETSC_COMM_WORLD, kspinfo, sizeof(kspinfo), &viewer));
320: PetscCall(KSPView(ksp, viewer));
321: PetscCall(PetscStrrchr(file[PetscPreLoadIt], '/', &matrixname));
322: PetscCall(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: PetscCall(PetscViewerDestroy(&viewer));
328: }
330: PetscCall(PetscOptionsGetString(NULL, NULL, "-solution", file[3], sizeof(file[3]), &flg));
331: if (flg) {
332: PetscViewer viewer;
333: Vec xstar;
335: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[3], FILE_MODE_READ, &viewer));
336: PetscCall(VecCreate(PETSC_COMM_WORLD, &xstar));
337: PetscCall(VecLoad(xstar, viewer));
338: PetscCall(VecAXPY(xstar, -1.0, x));
339: PetscCall(VecNorm(xstar, NORM_2, &enorm));
340: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)enorm));
341: PetscCall(VecDestroy(&xstar));
342: PetscCall(PetscViewerDestroy(&viewer));
343: }
344: if (outputSoln) {
345: PetscViewer viewer;
347: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "solution.petsc", FILE_MODE_WRITE, &viewer));
348: PetscCall(VecView(x, viewer));
349: PetscCall(PetscViewerDestroy(&viewer));
350: }
352: flg = PETSC_FALSE;
353: PetscCall(PetscOptionsGetBool(NULL, NULL, "-ksp_reason", &flg, NULL));
354: if (flg) {
355: KSPConvergedReason reason;
356: PetscCall(KSPGetConvergedReason(ksp, &reason));
357: PetscCall(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: PetscCall(MatDestroy(&A));
367: PetscCall(VecDestroy(&b));
368: PetscCall(VecDestroy(&u));
369: PetscCall(VecDestroy(&x));
370: PetscCall(VecDestroy(&b2));
371: PetscCall(KSPDestroy(&ksp));
372: if (flgB) PetscCall(MatDestroy(&B));
373: PetscPreLoadEnd();
374: /* -----------------------------------------------------------
375: End of linear solver loop
376: ----------------------------------------------------------- */
378: PetscCall(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*/