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, diagonalscale, 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"));
186: PetscCall(VecSet(x, 0.0));
188: if (ckerror) { /* Set true solution */
189: PetscCall(VecSet(u, 1.0));
190: PetscCall(MatMult(A, u, b));
191: }
193: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
194: Setup solve for system
195: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197: if (partition) {
198: MatPartitioning mpart;
199: IS mis, nis, is;
200: PetscInt *count;
201: PetscMPIInt size;
202: Mat BB;
203: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
204: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
205: PetscCall(PetscMalloc1(size, &count));
206: PetscCall(MatPartitioningCreate(PETSC_COMM_WORLD, &mpart));
207: PetscCall(MatPartitioningSetAdjacency(mpart, A));
208: /* PetscCall(MatPartitioningSetVertexWeights(mpart, weight)); */
209: PetscCall(MatPartitioningSetFromOptions(mpart));
210: PetscCall(MatPartitioningApply(mpart, &mis));
211: PetscCall(MatPartitioningDestroy(&mpart));
212: PetscCall(ISPartitioningToNumbering(mis, &nis));
213: PetscCall(ISPartitioningCount(mis, size, count));
214: PetscCall(ISDestroy(&mis));
215: PetscCall(ISInvertPermutation(nis, count[rank], &is));
216: PetscCall(PetscFree(count));
217: PetscCall(ISDestroy(&nis));
218: PetscCall(ISSort(is));
219: PetscCall(MatCreateSubMatrix(A, is, is, MAT_INITIAL_MATRIX, &BB));
221: /* need to move the vector also */
222: PetscCall(ISDestroy(&is));
223: PetscCall(MatDestroy(&A));
224: A = BB;
225: }
227: /*
228: Create linear solver; set operators; set runtime options.
229: */
230: PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp));
231: PetscCall(KSPSetInitialGuessNonzero(ksp, initialguess));
232: num_numfac = 1;
233: PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_numfac", &num_numfac, NULL));
234: while (num_numfac--) {
235: PetscCall(KSPSetOperators(ksp, A, A));
236: PetscCall(KSPSetFromOptions(ksp));
238: /*
239: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
240: enable more precise profiling of setting up the preconditioner.
241: These calls are optional, since both will be called within
242: KSPSolve() if they haven't been called already.
243: */
244: PetscCall(KSPSetUp(ksp));
245: PetscCall(KSPSetUpOnBlocks(ksp));
247: /*
248: Tests "diagonal-scaling of preconditioned residual norm" as used
249: by many ODE integrator codes including SUNDIALS. Note this is different
250: than diagonally scaling the matrix before computing the preconditioner
251: */
252: diagonalscale = PETSC_FALSE;
253: PetscCall(PetscOptionsGetBool(NULL, NULL, "-diagonal_scale", &diagonalscale, NULL));
254: if (diagonalscale) {
255: PC pc;
256: PetscInt j, start, end, n;
257: Vec scale;
259: PetscCall(KSPGetPC(ksp, &pc));
260: PetscCall(VecGetSize(x, &n));
261: PetscCall(VecDuplicate(x, &scale));
262: PetscCall(VecGetOwnershipRange(scale, &start, &end));
263: for (j = start; j < end; j++) PetscCall(VecSetValue(scale, j, ((PetscReal)(j + 1)) / ((PetscReal)n), INSERT_VALUES));
264: PetscCall(VecAssemblyBegin(scale));
265: PetscCall(VecAssemblyEnd(scale));
266: PetscCall(PCSetDiagonalScale(pc, scale));
267: PetscCall(VecDestroy(&scale));
268: }
270: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
271: Solve system
272: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
273: /*
274: Solve linear system;
275: */
276: if (trans) {
277: PetscCall(KSPSolveTranspose(ksp, b, x));
278: PetscCall(KSPGetIterationNumber(ksp, &its));
279: } else {
280: PetscInt num_rhs = 1;
281: PetscCall(PetscOptionsGetInt(NULL, NULL, "-num_rhs", &num_rhs, NULL));
283: while (num_rhs--) PetscCall(KSPSolve(ksp, b, x));
284: PetscCall(KSPGetIterationNumber(ksp, &its));
285: if (ckrnorm) { /* Check residual for each rhs */
286: if (trans) {
287: PetscCall(MatMultTranspose(A, x, b2));
288: } else {
289: PetscCall(MatMult(A, x, b2));
290: }
291: PetscCall(VecAXPY(b2, -1.0, b));
292: PetscCall(VecNorm(b2, NORM_2, &rnorm));
293: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Number of iterations = %3" PetscInt_FMT "\n", its));
294: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Residual norm %g\n", (double)rnorm));
295: }
296: if (ckerror && !trans) { /* Check error for each rhs */
297: /* PetscCall(VecView(x,PETSC_VIEWER_STDOUT_WORLD)); */
298: PetscCall(VecAXPY(u, -1.0, x));
299: PetscCall(VecNorm(u, NORM_2, &enorm));
300: PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Error norm %g\n", (double)enorm));
301: }
303: } /* while (num_rhs--) */
305: /*
306: Write output (optionally using table for solver details).
307: - PetscPrintf() handles output for multiprocessor jobs
308: by printing from only one processor in the communicator.
309: - KSPView() prints information about the linear solver.
310: */
311: if (table && ckrnorm) {
312: char *matrixname = NULL, kspinfo[120];
313: PetscViewer viewer;
315: /*
316: Open a string viewer; then write info to it.
317: */
318: PetscCall(PetscViewerStringOpen(PETSC_COMM_WORLD, kspinfo, sizeof(kspinfo), &viewer));
319: PetscCall(KSPView(ksp, viewer));
320: PetscCall(PetscStrrchr(file[PetscPreLoadIt], '/', &matrixname));
321: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%-8.8s %3" PetscInt_FMT " %2.0e %s \n", matrixname, its, (double)rnorm, kspinfo));
323: /*
324: Destroy the viewer
325: */
326: PetscCall(PetscViewerDestroy(&viewer));
327: }
329: PetscCall(PetscOptionsGetString(NULL, NULL, "-solution", file[3], sizeof(file[3]), &flg));
330: if (flg) {
331: PetscViewer viewer;
332: Vec xstar;
334: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file[3], FILE_MODE_READ, &viewer));
335: PetscCall(VecCreate(PETSC_COMM_WORLD, &xstar));
336: PetscCall(VecLoad(xstar, viewer));
337: PetscCall(VecAXPY(xstar, -1.0, x));
338: PetscCall(VecNorm(xstar, NORM_2, &enorm));
339: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Error norm %g\n", (double)enorm));
340: PetscCall(VecDestroy(&xstar));
341: PetscCall(PetscViewerDestroy(&viewer));
342: }
343: if (outputSoln) {
344: PetscViewer viewer;
346: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "solution.petsc", FILE_MODE_WRITE, &viewer));
347: PetscCall(VecView(x, viewer));
348: PetscCall(PetscViewerDestroy(&viewer));
349: }
351: flg = PETSC_FALSE;
352: PetscCall(PetscOptionsGetBool(NULL, NULL, "-ksp_reason", &flg, NULL));
353: if (flg) {
354: KSPConvergedReason reason;
355: PetscCall(KSPGetConvergedReason(ksp, &reason));
356: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "KSPConvergedReason: %s\n", KSPConvergedReasons[reason]));
357: }
359: } /* while (num_numfac--) */
361: /*
362: Free work space. All PETSc objects should be destroyed when they
363: are no longer needed.
364: */
365: PetscCall(MatDestroy(&A));
366: PetscCall(VecDestroy(&b));
367: PetscCall(VecDestroy(&u));
368: PetscCall(VecDestroy(&x));
369: PetscCall(VecDestroy(&b2));
370: PetscCall(KSPDestroy(&ksp));
371: if (flgB) PetscCall(MatDestroy(&B));
372: PetscPreLoadEnd();
373: /* -----------------------------------------------------------
374: End of linear solver loop
375: ----------------------------------------------------------- */
377: PetscCall(PetscFinalize());
378: return 0;
379: }
381: /*TEST
383: test:
384: args: -f ${DATAFILESPATH}/matrices/small -ksp_type preonly -pc_type ilu -pc_factor_mat_ordering_type natural -num_numfac 2 -pc_factor_reuse_fill
385: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
386: output_file: output/ex30.out
388: test:
389: TODO: Matrix row/column sizes are not compatible with block size
390: suffix: 2
391: 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
392: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
394: test:
395: suffix: shiftnz
396: 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
397: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
399: test:
400: suffix: shiftpd
401: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type POSITIVE_DEFINITE
402: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
404: test:
405: suffix: shift_cholesky_aij
406: 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
407: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
408: output_file: output/ex30_shiftnz.out
410: test:
411: suffix: shiftpd_2
412: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE
413: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
415: test:
416: suffix: shift_cholesky_sbaij
417: 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
418: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
419: output_file: output/ex30_shiftnz.out
421: test:
422: suffix: shiftpd_2_sbaij
423: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type POSITIVE_DEFINITE -mat_type sbaij
424: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
425: output_file: output/ex30_shiftpd_2.out
427: test:
428: suffix: shiftinblocks
429: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type lu -pc_factor_shift_type INBLOCKS
430: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
432: test:
433: suffix: shiftinblocks2
434: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS
435: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
436: output_file: output/ex30_shiftinblocks.out
438: test:
439: suffix: shiftinblockssbaij
440: args: -f0 ${DATAFILESPATH}/matrices/small -mat_sigma -4.0 -ksp_type preonly -pc_type cholesky -pc_factor_shift_type INBLOCKS -mat_type sbaij
441: requires: datafilespath !complex double !defined(PETSC_USE_64BIT_INDICES)
442: output_file: output/ex30_shiftinblocks.out
444: TEST*/