Actual source code: snestest.c
petsc-3.8.4 2018-03-24
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: }