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";

  4: /*
  5:   Include "petscksp.h" so that we can use KSP solvers.  Note that this file
  6:   automatically includes:
  7:      petscsys.h       - base PETSc routines   petscvec.h - vectors
  8:      petscmat.h - matrices
  9:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 10:      petscviewer.h - viewers               petscpc.h  - preconditioners
 11: */
 12: #include <petscksp.h>
 13: #include <petscviewerhdf5.h>

 15: static PetscErrorCode VecLoadIfExists_Private(Vec b,PetscViewer fd,PetscBool *has)
 16: {
 17:   PetscBool      hdf5=PETSC_FALSE;

 20:   PetscObjectTypeCompare((PetscObject)fd,PETSCVIEWERHDF5,&hdf5);
 21:   if (hdf5) {
 22: #if defined(PETSC_HAVE_HDF5)
 23:     PetscViewerHDF5HasObject(fd,(PetscObject)b,has);
 24:     if (*has) VecLoad(b,fd);
 25: #else
 26:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PETSc must be configured with HDF5 to use this feature");
 27: #endif
 28:   } else {
 29:     PetscErrorCode ierrp;
 30:     PetscPushErrorHandler(PetscReturnErrorHandler,NULL);
 31:     ierrp = VecLoad(b,fd);
 32:     PetscPopErrorHandler();
 33:     *has  = ierrp ? PETSC_FALSE : PETSC_TRUE;
 34:   }
 35:   return 0;
 36: }

 38: int main(int argc,char **args)
 39: {
 40:   KSP            ksp;             /* linear solver context */
 41:   Mat            A,N;                /* matrix */
 42:   Vec            x,b,r,Ab;          /* approx solution, RHS, residual */
 43:   PetscViewer    fd;               /* viewer */
 44:   char           file[PETSC_MAX_PATH_LEN]="";     /* input file name */
 45:   char           file_x0[PETSC_MAX_PATH_LEN]="";  /* name of input file with initial guess */
 46:   char           A_name[128]="A",b_name[128]="b",x0_name[128]="x0";  /* name of the matrix, RHS and initial guess */
 47:   KSPType        ksptype;
 48:   PetscBool      has;
 49:   PetscInt       its,n,m;
 50:   PetscReal      norm;
 51:   PetscBool      nonzero_guess=PETSC_TRUE;
 52:   PetscBool      solve_normal=PETSC_FALSE;
 53:   PetscBool      hdf5=PETSC_FALSE;
 54:   PetscBool      test_custom_layout=PETSC_FALSE;
 55:   PetscMPIInt    rank,size;

 57:   PetscInitialize(&argc,&args,(char*)0,help);
 58:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
 59:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 60:   /*
 61:      Determine files from which we read the linear system
 62:      (matrix, right-hand-side and initial guess vector).
 63:   */
 64:   PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),NULL);
 65:   PetscOptionsGetString(NULL,NULL,"-f_x0",file_x0,sizeof(file_x0),NULL);
 66:   PetscOptionsGetString(NULL,NULL,"-A_name",A_name,sizeof(A_name),NULL);
 67:   PetscOptionsGetString(NULL,NULL,"-b_name",b_name,sizeof(b_name),NULL);
 68:   PetscOptionsGetString(NULL,NULL,"-x0_name",x0_name,sizeof(x0_name),NULL);
 69:   /*
 70:      Decide whether to solve the original system (-solve_normal 0)
 71:      or the normal equation (-solve_normal 1).
 72:   */
 73:   PetscOptionsGetBool(NULL,NULL,"-solve_normal",&solve_normal,NULL);
 74:   /*
 75:      Decide whether to use the HDF5 reader.
 76:   */
 77:   PetscOptionsGetBool(NULL,NULL,"-hdf5",&hdf5,NULL);
 78:   /*
 79:      Decide whether custom matrix layout will be tested.
 80:   */
 81:   PetscOptionsGetBool(NULL,NULL,"-test_custom_layout",&test_custom_layout,NULL);

 83:   /* -----------------------------------------------------------
 84:                   Beginning of linear solver loop
 85:      ----------------------------------------------------------- */
 86:   /*
 87:      Loop through the linear solve 2 times.
 88:       - The intention here is to preload and solve a small system;
 89:         then load another (larger) system and solve it as well.
 90:         This process preloads the instructions with the smaller
 91:         system so that more accurate performance monitoring (via
 92:         -log_view) can be done with the larger one (that actually
 93:         is the system of interest).
 94:   */
 95:   PetscPreLoadBegin(PETSC_FALSE,"Load system");

 97:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
 98:                          Load system
 99:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

101:   /*
102:      Open binary file.  Note that we use FILE_MODE_READ to indicate
103:      reading from this file.
104:   */
105:   if (hdf5) {
106: #if defined(PETSC_HAVE_HDF5)
107:     PetscViewerHDF5Open(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
108:     PetscViewerPushFormat(fd,PETSC_VIEWER_HDF5_MAT);
109: #else
110:     SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"PETSc must be configured with HDF5 to use this feature");
111: #endif
112:   } else {
113:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
114:   }

116:   /*
117:      Load the matrix.
118:      Matrix type is set automatically but you can override it by MatSetType() prior to MatLoad().
119:      Do that only if you really insist on the given type.
120:   */
121:   MatCreate(PETSC_COMM_WORLD,&A);
122:   PetscObjectSetName((PetscObject)A,A_name);
123:   MatSetFromOptions(A);
124:   MatLoad(A,fd);
125:   if (test_custom_layout && size > 1) {
126:     /* Perturb the local sizes and create the matrix anew */
127:     PetscInt m1,n1;
128:     MatGetLocalSize(A,&m,&n);
129:     m = rank ? m-1 : m+size-1;
130:     n = (rank == size-1) ? n+size-1 : n-1;
131:     MatDestroy(&A);
132:     MatCreate(PETSC_COMM_WORLD,&A);
133:     PetscObjectSetName((PetscObject)A,A_name);
134:     MatSetSizes(A,m,n,PETSC_DECIDE,PETSC_DECIDE);
135:     MatSetFromOptions(A);
136:     MatLoad(A,fd);
137:     MatGetLocalSize(A,&m1,&n1);
139:   }
140:   MatGetLocalSize(A,&m,&n);

142:   /*
143:      Load the RHS vector if it is present in the file, otherwise use a vector of all ones.
144:   */
145:   MatCreateVecs(A, &x, &b);
146:   PetscObjectSetName((PetscObject)b,b_name);
147:   VecSetFromOptions(b);
148:   VecLoadIfExists_Private(b,fd,&has);
149:   if (!has) {
150:     PetscScalar one = 1.0;
151:     PetscPrintf(PETSC_COMM_WORLD,"Failed to load RHS, so use a vector of all ones.\n");
152:     VecSetFromOptions(b);
153:     VecSet(b,one);
154:   }

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

180:   VecDuplicate(x,&Ab);

182:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
183:                     Setup solve for system
184:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

186:   /*
187:      Conclude profiling last stage; begin profiling next stage.
188:   */
189:   PetscPreLoadStage("KSPSetUp");

191:   MatCreateNormalHermitian(A,&N);
192:   MatMultHermitianTranspose(A,b,Ab);

194:   /*
195:      Create linear solver; set operators; set runtime options.
196:   */
197:   KSPCreate(PETSC_COMM_WORLD,&ksp);

199:   if (solve_normal) {
200:     KSPSetOperators(ksp,N,N);
201:   } else {
202:     PC pc;
203:     KSPSetType(ksp,KSPLSQR);
204:     KSPGetPC(ksp,&pc);
205:     PCSetType(pc,PCNONE);
206:     KSPSetOperators(ksp,A,N);
207:   }
208:   KSPSetInitialGuessNonzero(ksp,nonzero_guess);
209:   KSPSetFromOptions(ksp);

211:   /*
212:      Here we explicitly call KSPSetUp() and KSPSetUpOnBlocks() to
213:      enable more precise profiling of setting up the preconditioner.
214:      These calls are optional, since both will be called within
215:      KSPSolve() if they haven't been called already.
216:   */
217:   KSPSetUp(ksp);
218:   KSPSetUpOnBlocks(ksp);

220:   /*
221:                          Solve system
222:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

224:   /*
225:      Begin profiling next stage
226:   */
227:   PetscPreLoadStage("KSPSolve");

229:   /*
230:      Solve linear system
231:   */
232:   if (solve_normal) {
233:     KSPSolve(ksp,Ab,x);
234:   } else {
235:     KSPSolve(ksp,b,x);
236:   }
237:   PetscObjectSetName((PetscObject)x,"x");

239:   /*
240:       Conclude profiling this stage
241:    */
242:   PetscPreLoadStage("Cleanup");

244:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
245:           Check error, print output, free data structures.
246:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

248:   /*
249:      Check error
250:   */
251:   VecDuplicate(b,&r);
252:   MatMult(A,x,r);
253:   VecAXPY(r,-1.0,b);
254:   VecNorm(r,NORM_2,&norm);
255:   KSPGetIterationNumber(ksp,&its);
256:   KSPGetType(ksp,&ksptype);
257:   PetscPrintf(PETSC_COMM_WORLD,"KSP type: %s\n",ksptype);
258:   PetscPrintf(PETSC_COMM_WORLD,"Number of iterations = %3D\n",its);
259:   PetscPrintf(PETSC_COMM_WORLD,"Residual norm %g\n",(double)norm);

261:   /*
262:      Free work space.  All PETSc objects should be destroyed when they
263:      are no longer needed.
264:   */
265:   MatDestroy(&A)); PetscCall(VecDestroy(&b);
266:   MatDestroy(&N)); PetscCall(VecDestroy(&Ab);
267:   VecDestroy(&r)); PetscCall(VecDestroy(&x);
268:   KSPDestroy(&ksp);
269:   PetscPreLoadEnd();
270:   /* -----------------------------------------------------------
271:                       End of linear solver loop
272:      ----------------------------------------------------------- */

274:   PetscFinalize();
275:   return 0;
276: }

278: /*TEST

280:    test:
281:       suffix: 1
282:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
283:       args: -f ${DATAFILESPATH}/matrices/medium -ksp_view -ksp_monitor_short -ksp_max_it 100 -solve_normal

285:    test:
286:       suffix: 2
287:       nsize: 2
288:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
289:       args: -f ${DATAFILESPATH}/matrices/shallow_water1 -ksp_view -ksp_monitor_short -ksp_max_it 100 -solve_normal

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

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

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

347:    # Test for correct cgls convergence reason
348:    test:
349:       suffix: 5
350:       nsize: 1
351:       requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
352:       args: -f ${DATAFILESPATH}/matrices/rectangular_ultrasound_4889x841
353:       args: -ksp_converged_reason -ksp_rtol 1e-2 -ksp_max_it 100
354:       args: -ksp_type cgls

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

377:    # Test correct handling of local dimensions in PCApply
378:    testset:
379:      requires: datafilespath double !complex !defined(PETSC_USE_64BIT_INDICES)
380:      requires: hdf5 defined(PETSC_HDF5_HAVE_ZLIB)
381:      nsize: 3
382:      suffix: 7
383:      args: -f ${DATAFILESPATH}/matrices/matlab/small.mat -hdf5 -test_custom_layout 1 -ksp_type lsqr -pc_type jacobi

385:    # Test complex matrices
386:    testset:
387:      requires: double complex !defined(PETSC_USE_64BIT_INDICES)
388:      args: -f ${wPETSC_DIR}/share/petsc/datafiles/matrices/nh-complex-int32-float64
389:      output_file: output/ex27_8.out
390:      filter: grep -v "KSP type"
391:      test:
392:        suffix: 8
393:        args: -solve_normal 0 -ksp_type {{lsqr cgls}}
394:      test:
395:        suffix: 8_normal
396:        args: -solve_normal 1 -ksp_type {{cg bicg}}

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

412: TEST*/