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