Actual source code: snestest.c

petsc-3.8.3 2017-12-09
Report Typos and Errors

  2:  #include <petsc/private/snesimpl.h>

  4: typedef struct {
  5:   PetscBool     complete_print;
  6:   PetscBool     threshold_print;
  7:   PetscScalar   threshold;
  8: } SNES_Test;


 11: PetscErrorCode SNESSolve_Test(SNES snes)
 12: {
 13:   Mat            A = snes->jacobian,B,C;
 14:   Vec            x = snes->vec_sol,f = snes->vec_func,f1 = snes->vec_sol_update;
 16:   PetscInt       i;
 17:   PetscReal      nrm,gnorm;
 18:   SNES_Test      *neP = (SNES_Test*)snes->data;
 19:   PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
 20:   void           *ctx;
 21:   PetscReal      fnorm,f1norm,dnorm;

 24:   if (A != snes->jacobian_pre) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot test with alternative preconditioner");

 26:   PetscPrintf(PetscObjectComm((PetscObject)snes),"Testing hand-coded Jacobian, if the ratio is\n");
 27:   PetscPrintf(PetscObjectComm((PetscObject)snes),"O(1.e-8), the hand-coded Jacobian is probably correct.\n");
 28:   if (!neP->complete_print) {
 29:     PetscPrintf(PetscObjectComm((PetscObject)snes),"Run with -snes_test_display to show difference\n");
 30:     PetscPrintf(PetscObjectComm((PetscObject)snes),"of hand-coded and finite difference Jacobian.\n");
 31:   }


 34:   for (i=0; i<3; i++) {
 35:     void                     *functx;
 36:     static const char *const loc[] = {"user-defined state","constant state -1.0","constant state 1.0"};
 37:     PetscInt                 m,n,M,N;

 39:     if (i == 1) {
 40:       VecSet(x,-1.0);
 41:     } else if (i == 2) {
 42:       VecSet(x,1.0);
 43:     }

 45:     /* evaluate the function at this point because SNESComputeJacobianDefaultColor() assumes that the function has been evaluated and put into snes->vec_func */
 46:     SNESComputeFunction(snes,x,f);
 47:     if (snes->domainerror) {
 48:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Domain error at %s\n",loc[i]);
 49:       snes->domainerror = PETSC_FALSE;
 50:       continue;
 51:     }

 53:     /* compute both versions of Jacobian */
 54:     SNESComputeJacobian(snes,x,A,A);

 56:     MatCreate(PetscObjectComm((PetscObject)A),&B);
 57:     MatGetSize(A,&M,&N);
 58:     MatGetLocalSize(A,&m,&n);
 59:     MatSetSizes(B,m,n,M,N);
 60:     MatSetType(B,((PetscObject)A)->type_name);
 61:     MatSetUp(B);
 62:     MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

 64:     SNESGetFunction(snes,NULL,NULL,&functx);
 65:     SNESComputeJacobianDefault(snes,x,B,B,functx);
 66:     if (neP->complete_print) {
 67:       MPI_Comm    comm;
 68:       PetscViewer viewer;
 69:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Finite difference Jacobian (%s)\n",loc[i]);
 70:       PetscObjectGetComm((PetscObject)B,&comm);
 71:       PetscViewerASCIIGetStdout(comm,&viewer);
 72:       MatView(B,viewer);
 73:     }
 74:     /* compare */
 75:     MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
 76:     MatNorm(B,NORM_FROBENIUS,&nrm);
 77:     MatNorm(A,NORM_FROBENIUS,&gnorm);
 78:     if (neP->complete_print) {
 79:       MPI_Comm    comm;
 80:       PetscViewer viewer;
 81:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Hand-coded Jacobian (%s)\n",loc[i]);
 82:       PetscObjectGetComm((PetscObject)B,&comm);
 83:       PetscViewerASCIIGetStdout(comm,&viewer);
 84:       MatView(A,viewer);
 85:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Hand-coded minus finite-difference Jacobian (%s)\n",loc[i]);
 86:       MatView(B,viewer);
 87:     }


 90:     if (neP->threshold_print) {
 91:       MPI_Comm          comm;
 92:       PetscViewer       viewer;
 93:       PetscInt          Istart, Iend, *ccols, bncols, cncols, j, row;
 94:       PetscScalar       *cvals;
 95:       const PetscInt    *bcols;
 96:       const PetscScalar *bvals;
 97: 
 98:       MatCreate(PetscObjectComm((PetscObject)A),&C);
 99:       MatSetSizes(C,m,n,M,N);
100:       MatSetType(C,((PetscObject)A)->type_name);
101:       MatSetUp(C);
102:       MatSetOption(C,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
103:       MatGetOwnershipRange(B,&Istart,&Iend);

105:       for (row = Istart; row < Iend; row++) {
106:         MatGetRow(B,row,&bncols,&bcols,&bvals);
107:         PetscMalloc2(bncols,&ccols,bncols,&cvals);
108:         for (j = 0, cncols = 0; j < bncols; j++) {
109:           if (PetscAbsScalar(bvals[j]) > PetscAbsScalar(neP->threshold)) {
110:             ccols[cncols] = bcols[j];
111:             cvals[cncols] = bvals[j];
112:             cncols += 1;
113:           }
114:         }
115:         if(cncols) {
116:           MatSetValues(C,1,&row,cncols,ccols,cvals,INSERT_VALUES);
117:         }
118:         MatRestoreRow(B,row,&bncols,&bcols,&bvals);
119:         PetscFree2(ccols,cvals);
120:       }
121: 
122:       MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
123:       MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
124: 
125:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Entries where difference is over threshold (%s)\n",loc[i]);
126:       PetscObjectGetComm((PetscObject)C,&comm);
127:       PetscViewerASCIIGetStdout(comm,&viewer);

129:       MatView(C,viewer);
130:       MatDestroy(&C);
131:     }

133:     if (!gnorm) gnorm = 1; /* just in case */
134:     PetscPrintf(PetscObjectComm((PetscObject)snes),"Norm of matrix ratio %g, difference %g (%s)\n",(double)(nrm/gnorm),(double)nrm,loc[i]);

136:     SNESGetObjective(snes,&objective,&ctx);
137:     if (objective) {
138:       SNESComputeFunction(snes,x,f);
139:       VecNorm(f,NORM_2,&fnorm);
140:       if (neP->complete_print) {
141:         PetscViewer viewer;
142:         PetscPrintf(PetscObjectComm((PetscObject)snes),"Hand-coded Function (%s)\n",loc[i]);
143:         PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
144:         VecView(f,viewer);
145:       }
146:       SNESObjectiveComputeFunctionDefaultFD(snes,x,f1,NULL);
147:       VecNorm(f1,NORM_2,&f1norm);
148:       if (neP->complete_print) {
149:         PetscViewer viewer;
150:         PetscPrintf(PetscObjectComm((PetscObject)snes),"Finite-difference Function (%s)\n",loc[i]);
151:         PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
152:         VecView(f1,viewer);
153:       }
154:       /* compare the two */
155:       VecAXPY(f,-1.0,f1);
156:       VecNorm(f,NORM_2,&dnorm);
157:       if (!fnorm) fnorm = 1.;
158:       PetscPrintf(PetscObjectComm((PetscObject)snes),"Norm of function ratio %g, difference %g (%s)\n",dnorm/fnorm,dnorm,loc[i]);
159:       if (neP->complete_print) {
160:         PetscViewer viewer;
161:         PetscPrintf(PetscObjectComm((PetscObject)snes),"Difference (%s)\n",loc[i]);
162:         PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)snes),&viewer);
163:         VecView(f,viewer);
164:       }
165:     }
166:     MatDestroy(&B);
167:   }

169:   /*
170:    Abort after the first iteration due to the jacobian not being valid.
171:   */

173:   SETERRQ(PetscObjectComm((PetscObject)snes),PETSC_ERR_ARG_WRONGSTATE,"SNESTest aborts after Jacobian test: it is NORMAL behavior.");
174:   return(0);
175: }


178: /* ------------------------------------------------------------ */
179: PetscErrorCode SNESDestroy_Test(SNES snes)
180: {

184:   PetscFree(snes->data);
185:   return(0);
186: }

188: static PetscErrorCode SNESSetFromOptions_Test(PetscOptionItems *PetscOptionsObject,SNES snes)
189: {
190:   SNES_Test      *ls = (SNES_Test*)snes->data;

194:   PetscOptionsHead(PetscOptionsObject,"Hand-coded Jacobian tester options");
195:   PetscOptionsBool("-snes_test_display","Display difference between hand-coded and finite difference Jacobians","None",ls->complete_print,&ls->complete_print,NULL);
196:   PetscOptionsScalar("-snes_test_display_threshold", "Display difference between hand-coded and finite difference Jacobians which exceed input threshold", "None", ls->threshold, &ls->threshold, &ls->threshold_print);
197:   PetscOptionsTail();
198:   return(0);
199: }

201: PetscErrorCode SNESSetUp_Test(SNES snes)
202: {

206:   SNESSetUpMatrices(snes);
207:   return(0);
208: }

210: /* ------------------------------------------------------------ */
211: /*MC
212:       SNESTEST - Test hand-coded Jacobian against finite difference Jacobian

214:    Options Database:
215: +  -snes_type test    - use a SNES solver that evaluates the difference between hand-code and finite-difference Jacobians
216: -  -snes_test_display - display the elements of the matrix, the difference between the Jacobian approximated by finite-differencing and hand-coded Jacobian

218:    Level: intermediate

220:    Notes: This solver is not a solver and does not converge to a solution.  SNESTEST checks the Jacobian at three
221:    points: the 0, 1, and -1 solution vectors.  At each point the following is reported.

223:    Output:
224: +  difference - ||J - Jd||, the norm of the difference of the hand-coded Jacobian J and the approximate Jacobian Jd obtained by finite-differencing
225:    the residual,
226: -  ratio      - ||J - Jd||/||J||, the ratio of the norms of the above difference and the hand-coded Jacobian.

228:    Frobenius norm is used in the above throughout. After doing these three tests, it always aborts with the error message
229:    "SNESTest aborts after Jacobian test".  No other behavior is to be expected.  It may be similarly used to check if a
230:    SNES function is the gradient of an objective function set with SNESSetObjective().

232: .seealso:  SNESCreate(), SNES, SNESSetType(), SNESUpdateCheckJacobian(), SNESNEWTONLS, SNESNEWTONTR

234: M*/
235: PETSC_EXTERN PetscErrorCode SNESCreate_Test(SNES snes)
236: {
237:   SNES_Test      *neP;

241:   snes->ops->solve          = SNESSolve_Test;
242:   snes->ops->destroy        = SNESDestroy_Test;
243:   snes->ops->setfromoptions = SNESSetFromOptions_Test;
244:   snes->ops->view           = 0;
245:   snes->ops->setup          = SNESSetUp_Test;
246:   snes->ops->reset          = 0;

248:   snes->usesksp = PETSC_FALSE;

250:   snes->alwayscomputesfinalresidual = PETSC_FALSE;

252:   PetscNewLog(snes,&neP);
253:   snes->data          = (void*)neP;
254:   neP->complete_print = PETSC_FALSE;
255:   return(0);
256: }

258: /*@C
259:     SNESUpdateCheckJacobian - Checks each Jacobian computed by the nonlinear solver comparing the users function with a finite difference computation.

261:    Options Database:
262: +    -snes_check_jacobian - use this every time SNESSolve() is called
263: -    -snes_check_jacobian_view -  Display difference between Jacobian approximated by finite-differencing and the hand-coded Jacobian

265:    Output:
266: +  difference - ||J - Jd||, the norm of the difference of the hand-coded Jacobian J and the approximate Jacobian Jd obtained by finite-differencing
267:    the residual,
268: -  ratio      - ||J - Jd||/||J||, the ratio of the norms of the above difference and the hand-coded Jacobian.

270:    Notes:
271:    Frobenius norm is used in the above throughout.  This check is carried out every SNES iteration.

273:    Level: intermediate

275: .seealso:  SNESTEST, SNESCreate(), SNES, SNESSetType(), SNESNEWTONLS, SNESNEWTONTR, SNESSolve()

277: @*/
278: PetscErrorCode SNESUpdateCheckJacobian(SNES snes,PetscInt it)
279: {
280:   Mat            A = snes->jacobian,B;
281:   Vec            x = snes->vec_sol,f = snes->vec_func,f1 = snes->vec_sol_update;
283:   PetscReal      nrm,gnorm;
284:   PetscErrorCode (*objective)(SNES,Vec,PetscReal*,void*);
285:   void           *ctx;
286:   PetscReal      fnorm,f1norm,dnorm;
287:   PetscInt       m,n,M,N;
288:   PetscBool      complete_print = PETSC_FALSE;
289:   void           *functx;
290:   PetscViewer    viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm((PetscObject)snes));

293:   PetscOptionsHasName(((PetscObject)snes)->options,((PetscObject)snes)->prefix,"-snes_check_jacobian_view",&complete_print);
294:   if (A != snes->jacobian_pre) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot check Jacobian with alternative preconditioner");

296:   PetscViewerASCIIAddTab(viewer,((PetscObject)snes)->tablevel);
297:   PetscViewerASCIIPrintf(viewer,"      Testing hand-coded Jacobian, if the ratio is O(1.e-8), the hand-coded Jacobian is probably correct.\n");
298:   if (!complete_print) {
299:     PetscViewerASCIIPrintf(viewer,"      Run with -snes_check_jacobian_view [viewer][:filename][:format] to show difference of hand-coded and finite difference Jacobian.\n");
300:   }

302:   /* compute both versions of Jacobian */
303:   SNESComputeJacobian(snes,x,A,A);

305:   MatCreate(PetscObjectComm((PetscObject)A),&B);
306:   MatGetSize(A,&M,&N);
307:   MatGetLocalSize(A,&m,&n);
308:   MatSetSizes(B,m,n,M,N);
309:   MatSetType(B,((PetscObject)A)->type_name);
310:   MatSetUp(B);
311:   MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
312:   SNESGetFunction(snes,NULL,NULL,&functx);
313:   SNESComputeJacobianDefault(snes,x,B,B,functx);

315:   if (complete_print) {
316:     PetscViewerASCIIPrintf(viewer,"    Finite difference Jacobian\n");
317:     MatViewFromOptions(B,(PetscObject)snes,"-snes_check_jacobian_view");
318:   }
319:   /* compare */
320:   MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
321:   MatNorm(B,NORM_FROBENIUS,&nrm);
322:   MatNorm(A,NORM_FROBENIUS,&gnorm);
323:   if (complete_print) {
324:     PetscViewerASCIIPrintf(viewer,"    Hand-coded Jacobian\n");
325:     MatViewFromOptions(A,(PetscObject)snes,"-snes_check_jacobian_view");
326:     PetscViewerASCIIPrintf(viewer,"    Hand-coded minus finite difference Jacobian\n");
327:     MatViewFromOptions(B,(PetscObject)snes,"-snes_check_jacobian_view");
328:   }
329:   if (!gnorm) gnorm = 1; /* just in case */
330:   PetscViewerASCIIPrintf(viewer,"    %g = ||J - Jfd||/||J|| %g  = ||J - Jfd||\n",(double)(nrm/gnorm),(double)nrm);

332:   SNESGetObjective(snes,&objective,&ctx);
333:   if (objective) {
334:     SNESComputeFunction(snes,x,f);
335:     VecNorm(f,NORM_2,&fnorm);
336:     if (complete_print) {
337:       PetscViewerASCIIPrintf(viewer,"    Hand-coded Objective Function \n");
338:       VecView(f,viewer);
339:     }
340:     SNESObjectiveComputeFunctionDefaultFD(snes,x,f1,NULL);
341:     VecNorm(f1,NORM_2,&f1norm);
342:     if (complete_print) {
343:       PetscViewerASCIIPrintf(viewer,"    Finite-Difference Objective Function\n");
344:       VecView(f1,viewer);
345:     }
346:     /* compare the two */
347:     VecAXPY(f,-1.0,f1);
348:     VecNorm(f,NORM_2,&dnorm);
349:     if (!fnorm) fnorm = 1.;
350:     PetscViewerASCIIPrintf(viewer,"    %g = Norm of objective function ratio %g = difference\n",dnorm/fnorm,dnorm);
351:     if (complete_print) {
352:       PetscViewerASCIIPrintf(viewer,"    Difference\n");
353:       VecView(f,viewer);
354:     }
355:   }
356:   PetscViewerASCIISubtractTab(viewer,((PetscObject)snes)->tablevel);

358:   MatDestroy(&B);
359:   return(0);
360: }