Actual source code: ex27.c
petsc-3.10.5 2019-03-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>
20: int main(int argc,char **args)
21: {
22: KSP ksp; /* linear solver context */
23: Mat A,N; /* matrix */
24: Vec x,b,r,Ab; /* approx solution, RHS, residual */
25: PetscViewer fd; /* viewer */
26: char file[PETSC_MAX_PATH_LEN]=""; /* input file name */
27: char file_x0[PETSC_MAX_PATH_LEN]=""; /* name of input file with initial guess */
28: KSPType ksptype;
29: PetscErrorCode ierr,ierrp;
30: PetscInt its,n,m;
31: PetscReal norm;
32: PetscBool nonzero_guess=PETSC_TRUE;
33: PetscBool solve_normal=PETSC_TRUE;
35: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
36: /*
37: Determine files from which we read the linear system
38: (matrix, right-hand-side and initial guess vector).
39: */
40: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);
41: PetscOptionsGetString(NULL,NULL,"-f_x0",file_x0,PETSC_MAX_PATH_LEN,NULL);
42: /*
43: Decide whether to solve the original system (-solve_normal 0)
44: or the normal equation (-solve_normal 1).
45: */
46: PetscOptionsGetBool(NULL,NULL,"-solve_normal",&solve_normal,NULL);
48: /* -----------------------------------------------------------
49: Beginning of linear solver loop
50: ----------------------------------------------------------- */
51: /*
52: Loop through the linear solve 2 times.
53: - The intention here is to preload and solve a small system;
54: then load another (larger) system and solve it as well.
55: This process preloads the instructions with the smaller
56: system so that more accurate performance monitoring (via
57: -log_view) can be done with the larger one (that actually
58: is the system of interest).
59: */
60: PetscPreLoadBegin(PETSC_FALSE,"Load system");
62: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
63: Load system
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
66: /*
67: Open binary file. Note that we use FILE_MODE_READ to indicate
68: reading from this file.
69: */
70: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
72: /*
73: Load the matrix and vector; then destroy the viewer.
74: */
75: MatCreate(PETSC_COMM_WORLD,&A);
76: MatSetType(A,MATMPIAIJ);
77: MatLoad(A,fd);
78: VecCreate(PETSC_COMM_WORLD,&b);
79: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
80: ierrp = VecLoad(b,fd);
81: PetscPopErrorHandler();
82: MatGetLocalSize(A,&m,&n);
83: if (ierrp) { /* if file contains no RHS, then use a vector of all ones */
84: PetscScalar one = 1.0;
85: PetscPrintf(PETSC_COMM_WORLD,"Failed to load RHS, so use a vector of all ones.\n");
86: VecSetSizes(b,m,PETSC_DECIDE);
87: VecSetFromOptions(b);
88: VecSet(b,one);
89: }
91: VecCreate(PETSC_COMM_WORLD,&x);
92: VecSetFromOptions(x);
94: /* load file_x0 if it is specified, otherwise try to reuse file */
95: if (file_x0[0]) {
96: PetscViewerDestroy(&fd);
97: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
98: }
99: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
100: ierrp = VecLoad(x,fd);
101: PetscPopErrorHandler();
102: PetscViewerDestroy(&fd);
103: if (ierrp) {
104: PetscPrintf(PETSC_COMM_WORLD,"Failed to load initial guess, so use a vector of all zeros.\n");
105: VecSetSizes(x,n,PETSC_DECIDE);
106: VecSet(x,0.0);
107: nonzero_guess=PETSC_FALSE;
108: }
110: VecDuplicate(x,&Ab);
112: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
113: Setup solve for system
114: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
116: /*
117: Conclude profiling last stage; begin profiling next stage.
118: */
119: PetscPreLoadStage("KSPSetUp");
121: MatCreateNormal(A,&N);
122: MatMultTranspose(A,b,Ab);
124: /*
125: Create linear solver; set operators; set runtime options.
126: */
127: KSPCreate(PETSC_COMM_WORLD,&ksp);
129: if (solve_normal) {
130: KSPSetOperators(ksp,N,N);
131: } else {
132: PC pc;
133: KSPSetType(ksp,KSPLSQR);
134: KSPGetPC(ksp,&pc);
135: PCSetType(pc,PCNONE);
136: KSPSetOperators(ksp,A,A);
137: }
138: KSPSetInitialGuessNonzero(ksp,nonzero_guess);
139: KSPSetFromOptions(ksp);
141: /*
142: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
143: enable more precise profiling of setting up the preconditioner.
144: These calls are optional, since both will be called within
145: KSPSolve() if they haven't been called already.
146: */
147: KSPSetUp(ksp);
148: KSPSetUpOnBlocks(ksp);
150: /*
151: Solve system
152: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: /*
155: Begin profiling next stage
156: */
157: PetscPreLoadStage("KSPSolve");
159: /*
160: Solve linear system
161: */
162: if (solve_normal) {
163: KSPSolve(ksp,Ab,x);
164: } else {
165: KSPSolve(ksp,b,x);
166: }
168: /*
169: Conclude profiling this stage
170: */
171: PetscPreLoadStage("Cleanup");
173: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
174: Check error, print output, free data structures.
175: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: /*
178: Check error
179: */
180: VecDuplicate(b,&r);
181: MatMult(A,x,r);
182: VecAXPY(r,-1.0,b);
183: VecNorm(r,NORM_2,&norm);
184: KSPGetIterationNumber(ksp,&its);
185: KSPGetType(ksp,&ksptype);
186: PetscPrintf(PETSC_COMM_WORLD,"KSP type: %s\n",ksptype);
187: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
188: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
190: /*
191: Free work space. All PETSc objects should be destroyed when they
192: are no longer needed.
193: */
194: MatDestroy(&A); VecDestroy(&b);
195: MatDestroy(&N); VecDestroy(&Ab);
196: VecDestroy(&r); VecDestroy(&x);
197: KSPDestroy(&ksp);
198: PetscPreLoadEnd();
199: /* -----------------------------------------------------------
200: End of linear solver loop
201: ----------------------------------------------------------- */
203: PetscFinalize();
204: return ierr;
205: }
209: /*TEST
211: test:
212: suffix: 1
213: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
214: args: -f ${DATAFILESPATH}/matrices/medium -ksp_view -ksp_monitor_short -ksp_max_it 100
216: test:
217: suffix: 2
218: nsize: 2
219: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
220: args: -f ${DATAFILESPATH}/matrices/shallow_water1 -ksp_view -ksp_monitor_short -ksp_max_it 100
222: # Test handling failing VecLoad without abort
223: test:
224: suffix: 3
225: nsize: {{1 2}separate output}
226: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
227: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
228: args: -f_x0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_x0
229: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
230: test:
231: suffix: 3a
232: nsize: {{1 2}separate output}
233: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
234: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
235: args: -f_x0 NONEXISTING_FILE
236: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
237: test:
238: suffix: 3b
239: nsize: {{1 2}separate output}
240: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
241: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0 # this file includes all A, b and x0
242: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
244: # Test least-square algorithms
245: test:
246: suffix: 4
247: nsize: {{1 2 4}}
248: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
249: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
250: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
251: args: -solve_normal 1 -ksp_type cg
252: test:
253: suffix: 4a
254: nsize: {{1 2 4}}
255: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
256: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
257: args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
258: args: -solve_normal 0 -ksp_type {{cgls lsqr}separate output}
259: test:
260: # Test KSPLSQR-specific options
261: suffix: 4b
262: nsize: 2
263: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
264: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
265: args: -ksp_converged_reason -ksp_rtol 1e-3 -ksp_max_it 200 -ksp_view
266: 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}
268: # Test for correct cgls convergence reason
269: test:
270: suffix: 5
271: nsize: 1
272: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
273: args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
274: args: -ksp_converged_reason -ksp_rtol 1e-2 -ksp_max_it 100
275: args: -solve_normal 0 -ksp_type cgls
277: TEST*/