Actual source code: ex27.c
petsc-3.11.4 2019-09-28
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(PetscIgnoreErrorHandler,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: KSPType ksptype;
54: PetscBool has;
55: PetscInt its,n,m;
56: PetscReal norm;
57: PetscBool nonzero_guess=PETSC_TRUE;
58: PetscBool solve_normal=PETSC_TRUE;
59: PetscBool hdf5=PETSC_FALSE;
60: PetscBool test_custom_layout=PETSC_FALSE;
61: PetscMPIInt rank,size;
63: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
64: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
65: MPI_Comm_size(PETSC_COMM_WORLD,&size);
66: /*
67: Determine files from which we read the linear system
68: (matrix, right-hand-side and initial guess vector).
69: */
70: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);
71: PetscOptionsGetString(NULL,NULL,"-f_x0",file_x0,PETSC_MAX_PATH_LEN,NULL);
72: /*
73: Decide whether to solve the original system (-solve_normal 0)
74: or the normal equation (-solve_normal 1).
75: */
76: PetscOptionsGetBool(NULL,NULL,"-solve_normal",&solve_normal,NULL);
77: /*
78: Decide whether to use the HDF5 reader.
79: */
80: PetscOptionsGetBool(NULL,NULL,"-hdf5",&hdf5,NULL);
81: /*
82: Decide whether custom matrix layout will be tested.
83: */
84: PetscOptionsGetBool(NULL,NULL,"-test_custom_layout",&test_custom_layout,NULL);
86: /* -----------------------------------------------------------
87: Beginning of linear solver loop
88: ----------------------------------------------------------- */
89: /*
90: Loop through the linear solve 2 times.
91: - The intention here is to preload and solve a small system;
92: then load another (larger) system and solve it as well.
93: This process preloads the instructions with the smaller
94: system so that more accurate performance monitoring (via
95: -log_view) can be done with the larger one (that actually
96: is the system of interest).
97: */
98: PetscPreLoadBegin(PETSC_FALSE,"Load system");
100: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
101: Load system
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104: /*
105: Open binary file. Note that we use FILE_MODE_READ to indicate
106: reading from this file.
107: */
108: if (hdf5) {
109: #if defined(PETSC_HAVE_HDF5)
110: PetscViewerHDF5Open(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
111: PetscViewerPushFormat(fd,PETSC_VIEWER_HDF5_MAT);
112: #else
113: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PETSc must be configured with HDF5 to use this feature");
114: #endif
115: } else {
116: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
117: }
119: /*
120: Load the matrix.
121: Matrix type is set automatically but you can override it by MatSetType() prior to MatLoad().
122: Do that only if you really insist on the given type.
123: */
124: MatCreate(PETSC_COMM_WORLD,&A);
125: PetscObjectSetName((PetscObject)A,"A");
126: MatSetFromOptions(A);
127: MatLoad(A,fd);
128: if (test_custom_layout) {
129: /* Perturb the local sizes and create the matrix anew */
130: PetscInt m1,n1;
131: MatGetLocalSize(A,&m,&n);
132: m = rank ? m-1 : m+size-1;
133: n = rank ? n-1 : n+size-1;
134: MatDestroy(&A);
135: MatCreate(PETSC_COMM_WORLD,&A);
136: PetscObjectSetName((PetscObject)A,"A");
137: MatSetSizes(A,m,n,PETSC_DECIDE,PETSC_DECIDE);
138: MatSetFromOptions(A);
139: MatLoad(A,fd);
140: MatGetLocalSize(A,&m1,&n1);
141: if (m1 != m || n1 != n) SETERRQ4(PETSC_COMM_WORLD,PETSC_ERR_SUP,"resulting sizes differ from demanded ones: %D %D != %D %D",m1,n1,m,n);
142: }
143: MatGetLocalSize(A,&m,&n);
145: /*
146: Load the RHS vector if it is present in the file, otherwise use a vector of all ones.
147: */
148: VecCreate(PETSC_COMM_WORLD,&b);
149: PetscObjectSetName((PetscObject)b,"b");
150: VecSetSizes(b,m,PETSC_DECIDE);
151: VecSetFromOptions(b);
152: VecLoadIfExists_Private(b,fd,&has);
153: if (!has) {
154: PetscScalar one = 1.0;
155: PetscPrintf(PETSC_COMM_WORLD,"Failed to load RHS, so use a vector of all ones.\n");
156: VecSetFromOptions(b);
157: VecSet(b,one);
158: }
160: /*
161: Load the initial guess vector if it is present in the file, otherwise use a vector of all zeros.
162: */
163: VecCreate(PETSC_COMM_WORLD,&x);
164: PetscObjectSetName((PetscObject)x,"x0");
165: VecSetSizes(x,n,PETSC_DECIDE);
166: VecSetFromOptions(x);
167: /* load file_x0 if it is specified, otherwise try to reuse file */
168: if (file_x0[0]) {
169: PetscViewerDestroy(&fd);
170: if (hdf5) {
171: #if defined(PETSC_HAVE_HDF5)
172: PetscViewerHDF5Open(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
173: #endif
174: } else {
175: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
176: }
177: }
178: VecLoadIfExists_Private(x,fd,&has);
179: if (!has) {
180: PetscPrintf(PETSC_COMM_WORLD,"Failed to load initial guess, so use a vector of all zeros.\n");
181: VecSet(x,0.0);
182: nonzero_guess=PETSC_FALSE;
183: }
184: PetscViewerDestroy(&fd);
186: VecDuplicate(x,&Ab);
188: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
189: Setup solve for system
190: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192: /*
193: Conclude profiling last stage; begin profiling next stage.
194: */
195: PetscPreLoadStage("KSPSetUp");
197: MatCreateNormal(A,&N);
198: MatMultTranspose(A,b,Ab);
200: /*
201: Create linear solver; set operators; set runtime options.
202: */
203: KSPCreate(PETSC_COMM_WORLD,&ksp);
205: if (solve_normal) {
206: KSPSetOperators(ksp,N,N);
207: } else {
208: PC pc;
209: KSPSetType(ksp,KSPLSQR);
210: KSPGetPC(ksp,&pc);
211: PCSetType(pc,PCNONE);
212: KSPSetOperators(ksp,A,A);
213: }
214: KSPSetInitialGuessNonzero(ksp,nonzero_guess);
215: KSPSetFromOptions(ksp);
217: /*
218: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
219: enable more precise profiling of setting up the preconditioner.
220: These calls are optional, since both will be called within
221: KSPSolve() if they haven't been called already.
222: */
223: KSPSetUp(ksp);
224: KSPSetUpOnBlocks(ksp);
226: /*
227: Solve system
228: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
230: /*
231: Begin profiling next stage
232: */
233: PetscPreLoadStage("KSPSolve");
235: /*
236: Solve linear system
237: */
238: if (solve_normal) {
239: KSPSolve(ksp,Ab,x);
240: } else {
241: KSPSolve(ksp,b,x);
242: }
243: PetscObjectSetName((PetscObject)x,"x");
245: /*
246: Conclude profiling this stage
247: */
248: PetscPreLoadStage("Cleanup");
250: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
251: Check error, print output, free data structures.
252: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
254: /*
255: Check error
256: */
257: VecDuplicate(b,&r);
258: MatMult(A,x,r);
259: VecAXPY(r,-1.0,b);
260: VecNorm(r,NORM_2,&norm);
261: KSPGetIterationNumber(ksp,&its);
262: KSPGetType(ksp,&ksptype);
263: PetscPrintf(PETSC_COMM_WORLD,"KSP type: %s\n",ksptype);
264: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
265: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
267: /*
268: Free work space. All PETSc objects should be destroyed when they
269: are no longer needed.
270: */
271: MatDestroy(&A); VecDestroy(&b);
272: MatDestroy(&N); VecDestroy(&Ab);
273: VecDestroy(&r); VecDestroy(&x);
274: KSPDestroy(&ksp);
275: PetscPreLoadEnd();
276: /* -----------------------------------------------------------
277: End of linear solver loop
278: ----------------------------------------------------------- */
280: PetscFinalize();
281: return ierr;
282: }
286: /*TEST
288: test:
289: suffix: 1
290: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
291: args: -f ${DATAFILESPATH}/matrices/medium -ksp_view -ksp_monitor_short -ksp_max_it 100
293: test:
294: suffix: 2
295: nsize: 2
296: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
297: args: -f ${DATAFILESPATH}/matrices/shallow_water1 -ksp_view -ksp_monitor_short -ksp_max_it 100
299: # Test handling failing VecLoad without abort
300: testset:
301: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
302: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
303: test:
304: suffix: 3
305: nsize: {{1 2}separate output}
306: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
307: args: -f_x0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_x0
308: test:
309: suffix: 3a
310: nsize: {{1 2}separate output}
311: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
312: args: -f_x0 NONEXISTING_FILE
313: test:
314: suffix: 3b
315: nsize: {{1 2}separate output}
316: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0 # this file includes all A, b and x0
317: test:
318: # Load square matrix, RHS and initial guess from HDF5 (Version 7.3 MAT-File)
319: suffix: 3b_hdf5
320: requires: hdf5 zlib
321: nsize: {{1 2}separate output}
322: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0.mat -hdf5
324: # Test least-square algorithms
325: testset:
326: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
327: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
328: test:
329: suffix: 4
330: nsize: {{1 2 4}}
331: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
332: args: -solve_normal 1 -ksp_type cg
333: test:
334: suffix: 4a
335: nsize: {{1 2 4}}
336: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
337: args: -solve_normal 0 -ksp_type {{cgls lsqr}separate output}
338: test:
339: # Test KSPLSQR-specific options
340: suffix: 4b
341: nsize: 2
342: args: -ksp_converged_reason -ksp_rtol 1e-3 -ksp_max_it 200 -ksp_view
343: args: -solve_normal 0 -ksp_type lsqr -ksp_convergence_test lsqr -ksp_lsqr_monitor -ksp_lsqr_compute_standard_error -ksp_lsqr_exact_mat_norm {{0 1}separate output}
345: test:
346: # Load rectangular matrix from HDF5 (Version 7.3 MAT-File)
347: suffix: 4a_lsqr_hdf5
348: nsize: {{1 2 4 8}}
349: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) hdf5 zlib
350: args: -f ${DATAFILESPATH}/matrices/matlab/rectangular_ultrasound_4889x841.mat -hdf5
351: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
352: args: -solve_normal 0 -ksp_type lsqr
353: args: -test_custom_layout {{0 1}}
355: # Test for correct cgls convergence reason
356: test:
357: suffix: 5
358: nsize: 1
359: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
360: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
361: args: -ksp_converged_reason -ksp_rtol 1e-2 -ksp_max_it 100
362: args: -solve_normal 0 -ksp_type cgls
364: TEST*/