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