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