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;

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

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

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

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

102:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
103:                          Load system
104:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

185:   VecDuplicate(x,&Ab);

187:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
188:                     Setup solve for system
189:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

196:   MatCreateNormalHermitian(A,&N);
197:   MatMultHermitianTranspose(A,b,Ab);

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

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

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

225:   /*
226:                          Solve system
227:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

229:   /*
230:      Begin profiling next stage
231:   */
232:   PetscPreLoadStage("KSPSolve");

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

244:   /*
245:       Conclude profiling this stage
246:    */
247:   PetscPreLoadStage("Cleanup");

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

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

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

279:   PetscFinalize();
280:   return 0;
281: }

283: /*TEST

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

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

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

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

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

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

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

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

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

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

417: TEST*/