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:   MatCreateNormalHermitian(A,&N);
199:   MatMultHermitianTranspose(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,N);
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: }

285: /*TEST

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

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

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

323:    # Test least-square algorithms
324:    testset:
325:      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
326:      args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
327:      test:
328:         suffix: 4
329:         nsize: {{1 2 4}}
330:         args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
331:         args: -solve_normal -ksp_type cg
332:      test:
333:         suffix: 4a
334:         nsize: {{1 2 4}}
335:         args: -ksp_converged_reason -ksp_monitor_short -ksp_rtol 1e-5 -ksp_max_it 100
336:         args: -ksp_type {{cgls lsqr}separate output}
337:      test:
338:         # Test KSPLSQR-specific options
339:         suffix: 4b
340:         nsize: 2
341:         args: -ksp_converged_reason -ksp_rtol 1e-3 -ksp_max_it 200 -ksp_view
342:         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}

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

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

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

384:    # Test correct handling of local dimensions in PCApply
385:    testset:
386:      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
387:      requires: hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
388:      nsize: 3
389:      suffix: 7
390:      args: -f ${DATAFILESPATH}/matrices/matlab/small.mat -hdf5 -test_custom_layout 1 -ksp_type lsqr -pc_type jacobi

392:    # Test complex matrices
393:    testset:
394:      requires: double complex !defined(PETSC_USE_64BIT_INDICES)
395:      args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64
396:      output_file: output/ex27_8.out
397:      filter: grep -v "KSP type"
398:      test:
399:        suffix: 8
400:        args: -solve_normal 0 -ksp_type {{lsqr cgls}}
401:      test:
402:        suffix: 8_normal
403:        args: -solve_normal 1 -ksp_type {{cg bicg}}

405:    testset:
406:      requires: double suitesparse !defined(PETSC_USE_64BIT_INDICES)
407:      args: -solve_normal {{0 1}shared output} -pc_type qr
408:      output_file: output/ex27_9.out
409:      filter: grep -v "KSP type"
410:      test:
411:        suffix: 9_real
412:        requires: !complex
413:        args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64
414:      test:
415:        suffix: 9_complex
416:        requires: complex
417:        args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64

419: TEST*/