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