Actual source code: ex27.c
2: static char help[] = "Reads a PETSc matrix and vector from a file and solves the normal equations.\n\n";
4: /*
5: Include "petscksp.h" so that we can use KSP solvers. Note that this file
6: automatically includes:
7: petscsys.h - base PETSc routines petscvec.h - vectors
8: petscmat.h - matrices
9: petscis.h - index sets petscksp.h - Krylov subspace methods
10: petscviewer.h - viewers petscpc.h - preconditioners
11: */
12: #include <petscksp.h>
13: #include <petscviewerhdf5.h>
15: static PetscErrorCode VecLoadIfExists_Private(Vec b,PetscViewer fd,PetscBool *has)
16: {
17: PetscBool hdf5=PETSC_FALSE;
20: PetscObjectTypeCompare((PetscObject)fd,PETSCVIEWERHDF5,&hdf5);
21: if (hdf5) {
22: #if defined(PETSC_HAVE_HDF5)
23: PetscViewerHDF5HasObject(fd,(PetscObject)b,has);
24: if (*has) VecLoad(b,fd);
25: #else
26: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PETSc must be configured with HDF5 to use this feature");
27: #endif
28: } else {
29: PetscErrorCode ierrp;
30: PetscPushErrorHandler(PetscReturnErrorHandler,NULL);
31: ierrp = VecLoad(b,fd);
32: PetscPopErrorHandler();
33: *has = ierrp ? PETSC_FALSE : PETSC_TRUE;
34: }
35: return 0;
36: }
38: int main(int argc,char **args)
39: {
40: KSP ksp; /* linear solver context */
41: Mat A,N; /* matrix */
42: Vec x,b,r,Ab; /* approx solution, RHS, residual */
43: PetscViewer fd; /* viewer */
44: char file[PETSC_MAX_PATH_LEN]=""; /* input file name */
45: char file_x0[PETSC_MAX_PATH_LEN]=""; /* name of input file with initial guess */
46: char A_name[128]="A",b_name[128]="b",x0_name[128]="x0"; /* name of the matrix, RHS and initial guess */
47: KSPType ksptype;
48: PetscBool has;
49: PetscInt its,n,m;
50: PetscReal norm;
51: PetscBool nonzero_guess=PETSC_TRUE;
52: PetscBool solve_normal=PETSC_FALSE;
53: PetscBool hdf5=PETSC_FALSE;
54: PetscBool test_custom_layout=PETSC_FALSE;
55: PetscMPIInt rank,size;
57: PetscInitialize(&argc,&args,(char*)0,help);
58: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
59: MPI_Comm_size(PETSC_COMM_WORLD,&size);
60: /*
61: Determine files from which we read the linear system
62: (matrix, right-hand-side and initial guess vector).
63: */
64: PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),NULL);
65: PetscOptionsGetString(NULL,NULL,"-f_x0",file_x0,sizeof(file_x0),NULL);
66: PetscOptionsGetString(NULL,NULL,"-A_name",A_name,sizeof(A_name),NULL);
67: PetscOptionsGetString(NULL,NULL,"-b_name",b_name,sizeof(b_name),NULL);
68: PetscOptionsGetString(NULL,NULL,"-x0_name",x0_name,sizeof(x0_name),NULL);
69: /*
70: Decide whether to solve the original system (-solve_normal 0)
71: or the normal equation (-solve_normal 1).
72: */
73: PetscOptionsGetBool(NULL,NULL,"-solve_normal",&solve_normal,NULL);
74: /*
75: Decide whether to use the HDF5 reader.
76: */
77: PetscOptionsGetBool(NULL,NULL,"-hdf5",&hdf5,NULL);
78: /*
79: Decide whether custom matrix layout will be tested.
80: */
81: PetscOptionsGetBool(NULL,NULL,"-test_custom_layout",&test_custom_layout,NULL);
83: /* -----------------------------------------------------------
84: Beginning of linear solver loop
85: ----------------------------------------------------------- */
86: /*
87: Loop through the linear solve 2 times.
88: - The intention here is to preload and solve a small system;
89: then load another (larger) system and solve it as well.
90: This process preloads the instructions with the smaller
91: system so that more accurate performance monitoring (via
92: -log_view) can be done with the larger one (that actually
93: is the system of interest).
94: */
95: PetscPreLoadBegin(PETSC_FALSE,"Load system");
97: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
98: Load system
99: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
101: /*
102: Open binary file. Note that we use FILE_MODE_READ to indicate
103: reading from this file.
104: */
105: if (hdf5) {
106: #if defined(PETSC_HAVE_HDF5)
107: PetscViewerHDF5Open(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
108: PetscViewerPushFormat(fd,PETSC_VIEWER_HDF5_MAT);
109: #else
110: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PETSc must be configured with HDF5 to use this feature");
111: #endif
112: } else {
113: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
114: }
116: /*
117: Load the matrix.
118: Matrix type is set automatically but you can override it by MatSetType() prior to MatLoad().
119: Do that only if you really insist on the given type.
120: */
121: MatCreate(PETSC_COMM_WORLD,&A);
122: PetscObjectSetName((PetscObject)A,A_name);
123: MatSetFromOptions(A);
124: MatLoad(A,fd);
125: if (test_custom_layout && size > 1) {
126: /* Perturb the local sizes and create the matrix anew */
127: PetscInt m1,n1;
128: MatGetLocalSize(A,&m,&n);
129: m = rank ? m-1 : m+size-1;
130: n = (rank == size-1) ? n+size-1 : n-1;
131: MatDestroy(&A);
132: MatCreate(PETSC_COMM_WORLD,&A);
133: PetscObjectSetName((PetscObject)A,A_name);
134: MatSetSizes(A,m,n,PETSC_DECIDE,PETSC_DECIDE);
135: MatSetFromOptions(A);
136: MatLoad(A,fd);
137: MatGetLocalSize(A,&m1,&n1);
139: }
140: MatGetLocalSize(A,&m,&n);
142: /*
143: Load the RHS vector if it is present in the file, otherwise use a vector of all ones.
144: */
145: MatCreateVecs(A, &x, &b);
146: PetscObjectSetName((PetscObject)b,b_name);
147: VecSetFromOptions(b);
148: VecLoadIfExists_Private(b,fd,&has);
149: if (!has) {
150: PetscScalar one = 1.0;
151: PetscPrintf(PETSC_COMM_WORLD,"Failed to load RHS, so use a vector of all ones.\n");
152: VecSetFromOptions(b);
153: VecSet(b,one);
154: }
156: /*
157: Load the initial guess vector if it is present in the file, otherwise use a vector of all zeros.
158: */
159: PetscObjectSetName((PetscObject)x,x0_name);
160: VecSetFromOptions(x);
161: /* load file_x0 if it is specified, otherwise try to reuse file */
162: if (file_x0[0]) {
163: PetscViewerDestroy(&fd);
164: if (hdf5) {
165: #if defined(PETSC_HAVE_HDF5)
166: PetscViewerHDF5Open(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
167: #endif
168: } else {
169: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
170: }
171: }
172: VecLoadIfExists_Private(x,fd,&has);
173: if (!has) {
174: PetscPrintf(PETSC_COMM_WORLD,"Failed to load initial guess, so use a vector of all zeros.\n");
175: VecSet(x,0.0);
176: nonzero_guess=PETSC_FALSE;
177: }
178: PetscViewerDestroy(&fd);
180: VecDuplicate(x,&Ab);
182: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
183: Setup solve for system
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: /*
187: Conclude profiling last stage; begin profiling next stage.
188: */
189: PetscPreLoadStage("KSPSetUp");
191: MatCreateNormalHermitian(A,&N);
192: MatMultHermitianTranspose(A,b,Ab);
194: /*
195: Create linear solver; set operators; set runtime options.
196: */
197: KSPCreate(PETSC_COMM_WORLD,&ksp);
199: if (solve_normal) {
200: KSPSetOperators(ksp,N,N);
201: } else {
202: PC pc;
203: KSPSetType(ksp,KSPLSQR);
204: KSPGetPC(ksp,&pc);
205: PCSetType(pc,PCNONE);
206: KSPSetOperators(ksp,A,N);
207: }
208: KSPSetInitialGuessNonzero(ksp,nonzero_guess);
209: KSPSetFromOptions(ksp);
211: /*
212: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
213: enable more precise profiling of setting up the preconditioner.
214: These calls are optional, since both will be called within
215: KSPSolve() if they haven't been called already.
216: */
217: KSPSetUp(ksp);
218: KSPSetUpOnBlocks(ksp);
220: /*
221: Solve system
222: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
224: /*
225: Begin profiling next stage
226: */
227: PetscPreLoadStage("KSPSolve");
229: /*
230: Solve linear system
231: */
232: if (solve_normal) {
233: KSPSolve(ksp,Ab,x);
234: } else {
235: KSPSolve(ksp,b,x);
236: }
237: PetscObjectSetName((PetscObject)x,"x");
239: /*
240: Conclude profiling this stage
241: */
242: PetscPreLoadStage("Cleanup");
244: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
245: Check error, print output, free data structures.
246: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
248: /*
249: Check error
250: */
251: VecDuplicate(b,&r);
252: MatMult(A,x,r);
253: VecAXPY(r,-1.0,b);
254: VecNorm(r,NORM_2,&norm);
255: KSPGetIterationNumber(ksp,&its);
256: KSPGetType(ksp,&ksptype);
257: PetscPrintf(PETSC_COMM_WORLD,"KSP type: %s\n",ksptype);
258: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
259: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
261: /*
262: Free work space. All PETSc objects should be destroyed when they
263: are no longer needed.
264: */
265: MatDestroy(&A)); PetscCall(VecDestroy(&b);
266: MatDestroy(&N)); PetscCall(VecDestroy(&Ab);
267: VecDestroy(&r)); PetscCall(VecDestroy(&x);
268: KSPDestroy(&ksp);
269: PetscPreLoadEnd();
270: /* -----------------------------------------------------------
271: End of linear solver loop
272: ----------------------------------------------------------- */
274: PetscFinalize();
275: return 0;
276: }
278: /*TEST
280: test:
281: suffix: 1
282: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
283: args: -f ${DATAFILESPATH}/matrices/medium -ksp_view -ksp_monitor_short -ksp_max_it 100 -solve_normal
285: test:
286: suffix: 2
287: nsize: 2
288: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
289: args: -f ${DATAFILESPATH}/matrices/shallow_water1 -ksp_view -ksp_monitor_short -ksp_max_it 100 -solve_normal
291: # Test handling failing VecLoad without abort
292: testset:
293: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
294: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
295: test:
296: suffix: 3
297: nsize: {{1 2}separate output}
298: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
299: args: -f_x0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_x0
300: test:
301: suffix: 3a
302: nsize: {{1 2}separate output}
303: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
304: args: -f_x0 NONEXISTING_FILE
305: test:
306: suffix: 3b
307: nsize: {{1 2}separate output}
308: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0 # this file includes all A, b and x0
309: test:
310: # Load square matrix, RHS and initial guess from HDF5 (Version 7.3 MAT-File)
311: suffix: 3b_hdf5
312: requires: hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
313: nsize: {{1 2}separate output}
314: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0.mat -hdf5
316: # Test least-square algorithms
317: testset:
318: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
319: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
320: test:
321: suffix: 4
322: nsize: {{1 2 4}}
323: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
324: args: -solve_normal -ksp_type cg
325: test:
326: suffix: 4a
327: nsize: {{1 2 4}}
328: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
329: args: -ksp_type {{cgls lsqr}separate output}
330: test:
331: # Test KSPLSQR-specific options
332: suffix: 4b
333: nsize: 2
334: args: -ksp_converged_reason -ksp_rtol 1e-3 -ksp_max_it 200 -ksp_view
335: args: -ksp_type lsqr -ksp_convergence_test lsqr -ksp_lsqr_monitor -ksp_lsqr_compute_standard_error -ksp_lsqr_exact_mat_norm {{0 1}separate output}
337: test:
338: # Load rectangular matrix from HDF5 (Version 7.3 MAT-File)
339: suffix: 4a_lsqr_hdf5
340: nsize: {{1 2 4 8}}
341: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
342: args: -f ${DATAFILESPATH}/matrices/matlab/rectangular_ultrasound_4889x841.mat -hdf5
343: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
344: args: -ksp_type lsqr
345: args: -test_custom_layout {{0 1}}
347: # Test for correct cgls convergence reason
348: test:
349: suffix: 5
350: nsize: 1
351: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
352: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
353: args: -ksp_converged_reason -ksp_rtol 1e-2 -ksp_max_it 100
354: args: -ksp_type cgls
356: # Load a matrix, RHS and solution from HDF5 (Version 7.3 MAT-File). Test immediate convergence.
357: testset:
358: nsize: {{1 2 4 8}}
359: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES) hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
360: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 10
361: args: -ksp_type lsqr
362: args: -test_custom_layout {{0 1}}
363: args: -hdf5 -x0_name x
364: test:
365: suffix: 6_hdf5
366: args: -f ${DATAFILESPATH}/matrices/matlab/small.mat
367: test:
368: suffix: 6_hdf5_rect
369: args: -f ${DATAFILESPATH}/matrices/matlab/small_rect.mat
370: test:
371: suffix: 6_hdf5_dense
372: args: -f ${DATAFILESPATH}/matrices/matlab/small_dense.mat -mat_type dense
373: test:
374: suffix: 6_hdf5_rect_dense
375: args: -f ${DATAFILESPATH}/matrices/matlab/small_rect_dense.mat -mat_type dense
377: # Test correct handling of local dimensions in PCApply
378: testset:
379: requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
380: requires: hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
381: nsize: 3
382: suffix: 7
383: args: -f ${DATAFILESPATH}/matrices/matlab/small.mat -hdf5 -test_custom_layout 1 -ksp_type lsqr -pc_type jacobi
385: # Test complex matrices
386: testset:
387: requires: double complex !defined(PETSC_USE_64BIT_INDICES)
388: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64
389: output_file: output/ex27_8.out
390: filter: grep -v "KSP type"
391: test:
392: suffix: 8
393: args: -solve_normal 0 -ksp_type {{lsqr cgls}}
394: test:
395: suffix: 8_normal
396: args: -solve_normal 1 -ksp_type {{cg bicg}}
398: testset:
399: requires: double suitesparse !defined(PETSC_USE_64BIT_INDICES)
400: args: -solve_normal {{0 1}shared output} -pc_type qr
401: output_file: output/ex27_9.out
402: filter: grep -v "KSP type"
403: test:
404: suffix: 9_real
405: requires: !complex
406: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64
407: test:
408: suffix: 9_complex
409: requires: complex
410: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64
412: TEST*/