Actual source code: ex27.c


  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_FALSE;
 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 && size > 1) {
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 == size-1) ? n+size-1 : n-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_PLIB,"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:   MatCreateVecs(A, &x, &b);
153:   PetscObjectSetName((PetscObject)b,b_name);
154:   VecSetFromOptions(b);
155:   VecLoadIfExists_Private(b,fd,&has);
156:   if (!has) {
157:     PetscScalar one = 1.0;
158:     PetscPrintf(PETSC_COMM_WORLD,"Failed to load RHS, so use a vector of all ones.\n");
159:     VecSetFromOptions(b);
160:     VecSet(b,one);
161:   }

163:   /*
164:      Load the initial guess vector if it is present in the file, otherwise use a vector of all zeros.
165:   */
166:   PetscObjectSetName((PetscObject)x,x0_name);
167:   VecSetFromOptions(x);
168:   /* load file_x0 if it is specified, otherwise try to reuse file */
169:   if (file_x0[0]) {
170:     PetscViewerDestroy(&fd);
171:     if (hdf5) {
172: #if defined(PETSC_HAVE_HDF5)
173:       PetscViewerHDF5Open(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
174: #endif
175:     } else {
176:       PetscViewerBinaryOpen(PETSC_COMM_WORLD,file_x0,FILE_MODE_READ,&fd);
177:     }
178:   }
179:   VecLoadIfExists_Private(x,fd,&has);
180:   if (!has) {
181:     PetscPrintf(PETSC_COMM_WORLD,"Failed to load initial guess, so use a vector of all zeros.\n");
182:     VecSet(x,0.0);
183:     nonzero_guess=PETSC_FALSE;
184:   }
185:   PetscViewerDestroy(&fd);

187:   VecDuplicate(x,&Ab);

189:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
190:                     Setup solve for system
191:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

193:   /*
194:      Conclude profiling last stage; begin profiling next stage.
195:   */
196:   PetscPreLoadStage("KSPSetUp");

198:   MatCreateNormal(A,&N);
199:   MatMultTranspose(A,b,Ab);

201:   /*
202:      Create linear solver; set operators; set runtime options.
203:   */
204:   KSPCreate(PETSC_COMM_WORLD,&ksp);

206:   if (solve_normal) {
207:     KSPSetOperators(ksp,N,N);
208:   } else {
209:     PC pc;
210:     KSPSetType(ksp,KSPLSQR);
211:     KSPGetPC(ksp,&pc);
212:     PCSetType(pc,PCNONE);
213:     KSPSetOperators(ksp,A,A);
214:   }
215:   KSPSetInitialGuessNonzero(ksp,nonzero_guess);
216:   KSPSetFromOptions(ksp);

218:   /*
219:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
220:      enable more precise profiling of setting up the preconditioner.
221:      These calls are optional, since both will be called within
222:      KSPSolve() if they haven't been called already.
223:   */
224:   KSPSetUp(ksp);
225:   KSPSetUpOnBlocks(ksp);

227:   /*
228:                          Solve system
229:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

231:   /*
232:      Begin profiling next stage
233:   */
234:   PetscPreLoadStage("KSPSolve");

236:   /*
237:      Solve linear system
238:   */
239:   if (solve_normal) {
240:     KSPSolve(ksp,Ab,x);
241:   } else {
242:     KSPSolve(ksp,b,x);
243:   }
244:   PetscObjectSetName((PetscObject)x,"x");

246:   /*
247:       Conclude profiling this stage
248:    */
249:   PetscPreLoadStage("Cleanup");

251:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
252:           Check error, print output, free data structures.
253:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

255:   /*
256:      Check error
257:   */
258:   VecDuplicate(b,&r);
259:   MatMult(A,x,r);
260:   VecAXPY(r,-1.0,b);
261:   VecNorm(r,NORM_2,&norm);
262:   KSPGetIterationNumber(ksp,&its);
263:   KSPGetType(ksp,&ksptype);
264:   PetscPrintf(PETSC_COMM_WORLD,"KSP type: %s\n",ksptype);
265:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
266:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);

268:   /*
269:      Free work space.  All PETSc objects should be destroyed when they
270:      are no longer needed.
271:   */
272:   MatDestroy(&A); VecDestroy(&b);
273:   MatDestroy(&N); VecDestroy(&Ab);
274:   VecDestroy(&r); VecDestroy(&x);
275:   KSPDestroy(&ksp);
276:   PetscPreLoadEnd();
277:   /* -----------------------------------------------------------
278:                       End of linear solver loop
279:      ----------------------------------------------------------- */

281:   PetscFinalize();
282:   return ierr;
283: }



287: /*TEST

289:    test:
290:       suffix: 1
291:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
292:       args: -f ${DATAFILESPATH}/matrices/medium -ksp_view -ksp_monitor_short -ksp_max_it 100 -solve_normal

294:    test:
295:       suffix: 2
296:       nsize: 2
297:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
298:       args: -f ${DATAFILESPATH}/matrices/shallow_water1 -ksp_view -ksp_monitor_short -ksp_max_it 100 -solve_normal

300:    # Test handling failing VecLoad without abort
301:    testset:
302:      requires: double !complex !define(PETSC_USE_64BIT_INDICES)
303:      args: -ksp_type cg -ksp_view -ksp_converged_reason -ksp_monitor_short -ksp_max_it 10
304:      test:
305:         suffix: 3
306:         nsize: {{1 2}separate output}
307:         args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
308:         args: -f_x0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_x0
309:      test:
310:         suffix: 3a
311:         nsize: {{1 2}separate output}
312:         args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system
313:         args: -f_x0 NONEXISTING_FILE
314:      test:
315:         suffix: 3b
316:         nsize: {{1 2}separate output}
317:         args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0  # this file includes all A, b and x0
318:      test:
319:         # Load square matrix, RHS and initial guess from HDF5 (Version 7.3 MAT-File)
320:         suffix: 3b_hdf5
321:         requires: hdf5 define(PETSC_HDF5_HAVE_ZLIB)
322:         nsize: {{1 2}separate output}
323:         args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/tiny_system_with_x0.mat -hdf5

325:    # Test least-square algorithms
326:    testset:
327:      requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
328:      args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
329:      test:
330:         suffix: 4
331:         nsize: {{1 2 4}}
332:         args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
333:         args: -solve_normal -ksp_type cg
334:      test:
335:         suffix: 4a
336:         nsize: {{1 2 4}}
337:         args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
338:         args: -ksp_type {{cgls lsqr}separate output}
339:      test:
340:         # Test KSPLSQR-specific options
341:         suffix: 4b
342:         nsize: 2
343:         args: -ksp_converged_reason -ksp_rtol 1e-3 -ksp_max_it 200 -ksp_view
344:         args: -ksp_type lsqr -ksp_convergence_test lsqr -ksp_lsqr_monitor -ksp_lsqr_compute_standard_error -ksp_lsqr_exact_mat_norm {{0 1}separate output}

346:    test:
347:       # Load rectangular matrix from HDF5 (Version 7.3 MAT-File)
348:       suffix: 4a_lsqr_hdf5
349:       nsize: {{1 2 4 8}}
350:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) hdf5 define(PETSC_HDF5_HAVE_ZLIB)
351:       args: -f ${DATAFILESPATH}/matrices/matlab/rectangular_ultrasound_4889x841.mat -hdf5
352:       args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
353:       args: -ksp_type lsqr
354:       args: -test_custom_layout {{0 1}}

356:    # Test for correct cgls convergence reason
357:    test:
358:       suffix: 5
359:       nsize: 1
360:       requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
361:       args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
362:       args: -ksp_converged_reason -ksp_rtol 1e-2 -ksp_max_it 100
363:       args: -ksp_type cgls


366:    # Load a matrix, RHS and solution from HDF5 (Version 7.3 MAT-File). Test immediate convergence.
367:    testset:
368:      nsize: {{1 2 4 8}}
369:      requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES) hdf5 define(PETSC_HDF5_HAVE_ZLIB)
370:      args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 10
371:      args: -ksp_type lsqr
372:      args: -test_custom_layout {{0 1}}
373:      args: -hdf5 -x0_name x
374:      test:
375:         suffix: 6_hdf5
376:         args: -f ${DATAFILESPATH}/matrices/matlab/small.mat
377:      test:
378:         suffix: 6_hdf5_rect
379:         args: -f ${DATAFILESPATH}/matrices/matlab/small_rect.mat
380:      test:
381:         suffix: 6_hdf5_dense
382:         args: -f ${DATAFILESPATH}/matrices/matlab/small_dense.mat -mat_type dense
383:      test:
384:         suffix: 6_hdf5_rect_dense
385:         args: -f ${DATAFILESPATH}/matrices/matlab/small_rect_dense.mat -mat_type dense

387:    # Test correct handling of local dimensions in PCApply
388:    testset:
389:      requires: datafilespath double !complex !define(PETSC_USE_64BIT_INDICES)
390:      requires: mumps hdf5 define(PETSC_HDF5_HAVE_ZLIB)
391:      nsize: 3
392:      args: -f ${DATAFILESPATH}/matrices/matlab/small.mat -hdf5
393:      args: -test_custom_layout 1
394:      test:
395:        suffix: 7a
396:        args: -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps
397:      test:
398:        suffix: 7b
399:        args: -ksp_type lsqr -pc_type none

401: TEST*/