Actual source code: ex27.c

petsc-master 2018-05-22
Report Typos and Errors

  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}

260:    # Test for correct cgls convergence reason
261:    test:
262:       suffix: 5
263:       nsize: 1
264:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
265:       args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
266:       args: -ksp_converged_reason -ksp_rtol 1e-2 -ksp_max_it 100
267:       args: -solve_normal 0 -ksp_type cgls

269: TEST*/