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