Actual source code: ex27.c

petsc-3.11.3 2019-06-26
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>
 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(PetscIgnoreErrorHandler,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:   KSPType        ksptype;
 54:   PetscBool      has;
 55:   PetscInt       its,n,m;
 56:   PetscReal      norm;
 57:   PetscBool      nonzero_guess=PETSC_TRUE;
 58:   PetscBool      solve_normal=PETSC_TRUE;
 59:   PetscBool      hdf5=PETSC_FALSE;
 60:   PetscBool      test_custom_layout=PETSC_FALSE;
 61:   PetscMPIInt    rank,size;

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

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

100:   /* - - - - - - - - - - - New Stage - - - - - - - - - - - - -
101:                          Load system
102:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

186:   VecDuplicate(x,&Ab);

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

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

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

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

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

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

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

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

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

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

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

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

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

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



286: /*TEST

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

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

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

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

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

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

364: TEST*/