Actual source code: ex27.c
petsc-main 2021-04-20
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_TRUE;
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) {
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 ? n-1 : n+size-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_SUP,"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: VecCreate(PETSC_COMM_WORLD,&b);
153: PetscObjectSetName((PetscObject)b,b_name);
154: VecSetSizes(b,m,PETSC_DECIDE);
155: VecSetFromOptions(b);
156: VecLoadIfExists_Private(b,fd,&has);
157: if (!has) {
158: PetscScalar one = 1.0;
159: PetscPrintf(PETSC_COMM_WORLD,"Failed to load RHS, so use a vector of all ones.\n");
160: VecSetFromOptions(b);
161: VecSet(b,one);
162: }
164: /*
165: Load the initial guess vector if it is present in the file, otherwise use a vector of all zeros.
166: */
167: VecCreate(PETSC_COMM_WORLD,&x);
168: PetscObjectSetName((PetscObject)x,x0_name);
169: VecSetSizes(x,n,PETSC_DECIDE);
170: VecSetFromOptions(x);
171: /* load file_x0 if it is specified, otherwise try to reuse file */
172: if (file_x0[0]) {
173: PetscViewerDestroy(&fd);
174: if (hdf5) {
175: #if defined(PETSC_HAVE_HDF5)
176: PetscViewerHDF5Open(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
177: #endif
178: } else {
179: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
180: }
181: }
182: VecLoadIfExists_Private(x,fd,&has);
183: if (!has) {
184: PetscPrintf(PETSC_COMM_WORLD,"Failed to load initial guess, so use a vector of all zeros.\n");
185: VecSet(x,0.0);
186: nonzero_guess=PETSC_FALSE;
187: }
188: PetscViewerDestroy(&fd);
190: VecDuplicate(x,&Ab);
192: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
193: Setup solve for system
194: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
196: /*
197: Conclude profiling last stage; begin profiling next stage.
198: */
199: PetscPreLoadStage("KSPSetUp");
201: MatCreateNormal(A,&N);
202: MatMultTranspose(A,b,Ab);
204: /*
205: Create linear solver; set operators; set runtime options.
206: */
207: KSPCreate(PETSC_COMM_WORLD,&ksp);
209: if (solve_normal) {
210: KSPSetOperators(ksp,N,N);
211: } else {
212: PC pc;
213: KSPSetType(ksp,KSPLSQR);
214: KSPGetPC(ksp,&pc);
215: PCSetType(pc,PCNONE);
216: KSPSetOperators(ksp,A,A);
217: }
218: KSPSetInitialGuessNonzero(ksp,nonzero_guess);
219: KSPSetFromOptions(ksp);
221: /*
222: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
223: enable more precise profiling of setting up the preconditioner.
224: These calls are optional, since both will be called within
225: KSPSolve() if they haven't been called already.
226: */
227: KSPSetUp(ksp);
228: KSPSetUpOnBlocks(ksp);
230: /*
231: Solve system
232: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
234: /*
235: Begin profiling next stage
236: */
237: PetscPreLoadStage("KSPSolve");
239: /*
240: Solve linear system
241: */
242: if (solve_normal) {
243: KSPSolve(ksp,Ab,x);
244: } else {
245: KSPSolve(ksp,b,x);
246: }
247: PetscObjectSetName((PetscObject)x,"x");
249: /*
250: Conclude profiling this stage
251: */
252: PetscPreLoadStage("Cleanup");
254: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
255: Check error, print output, free data structures.
256: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
258: /*
259: Check error
260: */
261: VecDuplicate(b,&r);
262: MatMult(A,x,r);
263: VecAXPY(r,-1.0,b);
264: VecNorm(r,NORM_2,&norm);
265: KSPGetIterationNumber(ksp,&its);
266: KSPGetType(ksp,&ksptype);
267: PetscPrintf(PETSC_COMM_WORLD,"KSP type: %s\n",ksptype);
268: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
269: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
271: /*
272: Free work space. All PETSc objects should be destroyed when they
273: are no longer needed.
274: */
275: MatDestroy(&A); VecDestroy(&b);
276: MatDestroy(&N); VecDestroy(&Ab);
277: VecDestroy(&r); VecDestroy(&x);
278: KSPDestroy(&ksp);
279: PetscPreLoadEnd();
280: /* -----------------------------------------------------------
281: End of linear solver loop
282: ----------------------------------------------------------- */
284: PetscFinalize();
285: return ierr;
286: }
290: /*TEST
292: test:
293: suffix: 1
294: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
295: args: -f ${DATAFILESPATH}/matrices/medium -ksp_view -ksp_monitor_short -ksp_max_it 100
297: test:
298: suffix: 2
299: nsize: 2
300: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
301: args: -f ${DATAFILESPATH}/matrices/shallow_water1 -ksp_view -ksp_monitor_short -ksp_max_it 100
303: # Test handling failing VecLoad without abort
304: testset:
305: requires: double !complex !define(PETSC_USE_64BIT_INDICES)
306: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
307: test:
308: suffix: 3
309: nsize: {{1 2}separate output}
310: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
311: args: -f_x0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_x0
312: test:
313: suffix: 3a
314: nsize: {{1 2}separate output}
315: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
316: args: -f_x0 NONEXISTING_FILE
317: test:
318: suffix: 3b
319: nsize: {{1 2}separate output}
320: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0 # this file includes all A, b and x0
321: test:
322: # Load square matrix, RHS and initial guess from HDF5 (Version 7.3 MAT-File)
323: suffix: 3b_hdf5
324: requires: hdf5 define(PETSC_HDF5_HAVE_ZLIB)
325: nsize: {{1 2}separate output}
326: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0.mat -hdf5
328: # Test least-square algorithms
329: testset:
330: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
331: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
332: test:
333: suffix: 4
334: nsize: {{1 2 4}}
335: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
336: args: -solve_normal 1 -ksp_type cg
337: test:
338: suffix: 4a
339: nsize: {{1 2 4}}
340: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
341: args: -solve_normal 0 -ksp_type {{cgls lsqr}separate output}
342: test:
343: # Test KSPLSQR-specific options
344: suffix: 4b
345: nsize: 2
346: args: -ksp_converged_reason -ksp_rtol 1e-3 -ksp_max_it 200 -ksp_view
347: 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}
349: test:
350: # Load rectangular matrix from HDF5 (Version 7.3 MAT-File)
351: suffix: 4a_lsqr_hdf5
352: nsize: {{1 2 4 8}}
353: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) hdf5 define(PETSC_HDF5_HAVE_ZLIB)
354: args: -f ${DATAFILESPATH}/matrices/matlab/rectangular_ultrasound_4889x841.mat -hdf5
355: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
356: args: -solve_normal 0 -ksp_type lsqr
357: args: -test_custom_layout {{0 1}}
359: # Test for correct cgls convergence reason
360: test:
361: suffix: 5
362: nsize: 1
363: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
364: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
365: args: -ksp_converged_reason -ksp_rtol 1e-2 -ksp_max_it 100
366: args: -solve_normal 0 -ksp_type cgls
369: # Load a matrix, RHS and solution from HDF5 (Version 7.3 MAT-File). Test immediate convergence.
370: testset:
371: nsize: {{1 2 4 8}}
372: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) hdf5 define(PETSC_HDF5_HAVE_ZLIB)
373: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 10
374: args: -solve_normal 0 -ksp_type lsqr
375: args: -test_custom_layout {{0 1}}
376: args: -hdf5 -x0_name x
377: test:
378: suffix: 6_hdf5
379: args: -f ${DATAFILESPATH}/matrices/matlab/small.mat
380: test:
381: suffix: 6_hdf5_rect
382: args: -f ${DATAFILESPATH}/matrices/matlab/small_rect.mat
383: test:
384: suffix: 6_hdf5_dense
385: args: -f ${DATAFILESPATH}/matrices/matlab/small_dense.mat -mat_type dense
386: test:
387: suffix: 6_hdf5_rect_dense
388: args: -f ${DATAFILESPATH}/matrices/matlab/small_rect_dense.mat -mat_type dense
390: TEST*/