Actual source code: ex27.c
petsc-3.9.4 2018-09-11
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,u,Ab; /* approx solution, RHS, exact solution */
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: PetscErrorCode ierr,ierrp;
29: PetscInt its,n,m;
30: PetscReal norm;
31: PetscBool nonzero_guess=PETSC_TRUE;
33: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
34: /*
35: Determine files from which we read the linear system
36: (matrix, right-hand-side and initial guess vector).
37: */
38: PetscOptionsGetString(NULL,NULL,"-f",file,PETSC_MAX_PATH_LEN,NULL);
39: PetscOptionsGetString(NULL,NULL,"-f_x0",file_x0,PETSC_MAX_PATH_LEN,NULL);
41: /* -----------------------------------------------------------
42: Beginning of linear solver loop
43: ----------------------------------------------------------- */
44: /*
45: Loop through the linear solve 2 times.
46: - The intention here is to preload and solve a small system;
47: then load another (larger) system and solve it as well.
48: This process preloads the instructions with the smaller
49: system so that more accurate performance monitoring (via
50: -log_view) can be done with the larger one (that actually
51: is the system of interest).
52: */
53: PetscPreLoadBegin(PETSC_FALSE,"Load system");
55: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
56: Load system
57: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
59: /*
60: Open binary file. Note that we use FILE_MODE_READ to indicate
61: reading from this file.
62: */
63: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
65: /*
66: Load the matrix and vector; then destroy the viewer.
67: */
68: MatCreate(PETSC_COMM_WORLD,&A);
69: MatSetType(A,MATMPIAIJ);
70: MatLoad(A,fd);
71: VecCreate(PETSC_COMM_WORLD,&b);
72: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
73: ierrp = VecLoad(b,fd);
74: PetscPopErrorHandler();
75: MatGetLocalSize(A,&m,&n);
76: if (ierrp) { /* if file contains no RHS, then use a vector of all ones */
77: PetscScalar one = 1.0;
78: PetscPrintf(PETSC_COMM_WORLD,"Failed to load RHS, so use a vector of all ones.\n");
79: VecSetSizes(b,m,PETSC_DECIDE);
80: VecSetFromOptions(b);
81: VecSet(b,one);
82: }
84: VecCreate(PETSC_COMM_WORLD,&x);
85: VecSetFromOptions(x);
87: /* load file_x0 if it is specified, otherwise try to reuse file */
88: if (file_x0[0]) {
89: PetscViewerDestroy(&fd);
90: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
91: }
92: PetscPushErrorHandler(PetscIgnoreErrorHandler,NULL);
93: ierrp = VecLoad(x,fd);
94: PetscPopErrorHandler();
95: PetscViewerDestroy(&fd);
96: if (ierrp) {
97: PetscPrintf(PETSC_COMM_WORLD,"Failed to load initial guess, so use a vector of all zeros.\n");
98: VecSetSizes(x,n,PETSC_DECIDE);
99: VecSet(x,0.0);
100: nonzero_guess=PETSC_FALSE;
101: }
103: VecDuplicate(x,&u);
104: VecDuplicate(x,&Ab);
106: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
107: Setup solve for system
108: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: /*
111: Conclude profiling last stage; begin profiling next stage.
112: */
113: PetscPreLoadStage("KSPSetUp");
115: MatCreateNormal(A,&N);
116: MatMultTranspose(A,b,Ab);
118: /*
119: Create linear solver; set operators; set runtime options.
120: */
121: KSPCreate(PETSC_COMM_WORLD,&ksp);
123: KSPSetOperators(ksp,N,N);
124: KSPSetInitialGuessNonzero(ksp,nonzero_guess);
125: KSPSetFromOptions(ksp);
127: /*
128: Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
129: enable more precise profiling of setting up the preconditioner.
130: These calls are optional, since both will be called within
131: KSPSolve() if they haven't been called already.
132: */
133: KSPSetUp(ksp);
134: KSPSetUpOnBlocks(ksp);
136: /*
137: Solve system
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: /*
141: Begin profiling next stage
142: */
143: PetscPreLoadStage("KSPSolve");
145: /*
146: Solve linear system
147: */
148: KSPSolve(ksp,Ab,x);
150: /*
151: Conclude profiling this stage
152: */
153: PetscPreLoadStage("Cleanup");
155: /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
156: Check error, print output, free data structures.
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159: /*
160: Check error
161: */
162: MatMult(A,x,u);
163: VecAXPY(u,-1.0,b);
164: VecNorm(u,NORM_2,&norm);
165: KSPGetIterationNumber(ksp,&its);
166: PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
167: PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);
169: /*
170: Free work space. All PETSc objects should be destroyed when they
171: are no longer needed.
172: */
173: MatDestroy(&A); VecDestroy(&b);
174: MatDestroy(&N); VecDestroy(&Ab);
175: VecDestroy(&u); VecDestroy(&x);
176: KSPDestroy(&ksp);
177: PetscPreLoadEnd();
178: /* -----------------------------------------------------------
179: End of linear solver loop
180: ----------------------------------------------------------- */
182: PetscFinalize();
183: return ierr;
184: }
188: /*TEST
190: test:
191: suffix: 1
192: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
193: args: -f ${DATAFILESPATH}/matrices/medium -ksp_view -ksp_monitor_short -ksp_max_it 100
195: test:
196: suffix: 2
197: nsize: 2
198: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
199: args: -f ${DATAFILESPATH}/matrices/shallow_water1 -ksp_view -ksp_monitor_short -ksp_max_it 100
201: # Test handling failing VecLoad without abort
202: test:
203: suffix: 3
204: nsize: {{1 2}separate output}
205: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
206: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
207: args: -f_x0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_x0
208: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
209: test:
210: suffix: 3a
211: nsize: {{1 2}separate output}
212: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
213: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
214: args: -f_x0 NONEXISTING_FILE
215: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
216: test:
217: suffix: 3b
218: nsize: {{1 2}separate output}
219: requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
220: args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0 # this file includes all A, b and x0
221: args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
224: TEST*/