Actual source code: fdtest.c
1: #include "private/taosolver_impl.h" /*I "taosolver.h" I*/
3: typedef struct {
4: PetscBool check_gradient;
5: PetscBool check_hessian;
6: PetscBool complete_print;
7: } FD_Test;
9: /*
10: TaoSolve_FD - Tests whether a hand computed Hessian
11: matches one compute via finite differences.
12: */
15: PetscErrorCode TaoSolve_FD(TaoSolver tao)
16: {
17: Mat A = tao->hessian,B;
18: Vec x = tao->solution,g1,g2;
20: PetscInt i;
21: MatStructure flg;
22: PetscReal nrm,gnorm,hcnorm,fdnorm;
23: MPI_Comm comm;
24: FD_Test *fd = (FD_Test*)tao->data;
27: comm = ((PetscObject)tao)->comm;
28: if (fd->check_gradient) {
29: VecDuplicate(x,&g1);
30: VecDuplicate(x,&g2);
32: PetscPrintf(comm,"Testing hand-coded gradient (hc) against finite difference gradient (fd), if the ratio ||fd - hc|| / ||hc|| is\n");
33: PetscPrintf(comm,"0 (1.e-8), the hand-coded gradient is probably correct.\n");
34:
35: if (!fd->complete_print) {
36: PetscPrintf(comm,"Run with -tao_fd_test_display to show difference\n");
37: PetscPrintf(comm,"between hand-coded and finite difference gradient.\n");
38: }
39: for (i=0; i<3; i++) {
40: if (i == 1) {VecSet(x,-1.0);}
41: else if (i == 2) {VecSet(x,1.0);}
42:
43: /* Compute both version of gradient */
44: TaoComputeGradient(tao,x,g1);
45: TaoDefaultComputeGradient(tao,x,g2,PETSC_NULL);
46: if (fd->complete_print) {
47: MPI_Comm gcomm;
48: PetscViewer viewer;
49: PetscPrintf(comm,"Finite difference gradient\n");
50: PetscObjectGetComm((PetscObject)g2,&gcomm);
51: PetscViewerASCIIGetStdout(gcomm,&viewer);
52: VecView(g2,viewer);
53: PetscPrintf(comm,"Hand-coded gradient\n");
54: PetscObjectGetComm((PetscObject)g1,&gcomm);
55: PetscViewerASCIIGetStdout(gcomm,&viewer);
56: VecView(g1,viewer);
57: PetscPrintf(comm,"\n");
58: }
59:
60: VecAXPY(g2,-1.0,g1);
61: VecNorm(g1,NORM_2,&hcnorm);
62: VecNorm(g2,NORM_2,&fdnorm);
63:
64: if (!hcnorm) hcnorm=1.0e-20;
65: PetscPrintf(comm,"ratio ||fd-hc||/||hc|| = %G, difference ||fd-hc|| = %G\n", fdnorm/hcnorm, fdnorm);
67: }
68: VecDestroy(&g1);
69: VecDestroy(&g2);
70: }
75: if (fd->check_hessian) {
76: if (A != tao->hessian_pre) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot test with alternative preconditioner");
78: PetscPrintf(comm,"Testing hand-coded Hessian (hc) against finite difference Hessian (fd). If the ratio is\n");
79: PetscPrintf(comm,"O (1.e-8), the hand-coded Hessian is probably correct.\n");
80:
81: if (!fd->complete_print) {
82: PetscPrintf(comm,"Run with -tao_fd_test_display to show difference\n");
83: PetscPrintf(comm,"of hand-coded and finite difference Hessian.\n");
84: }
85: for (i=0;i<3;i++) {
86: /* compute both versions of Hessian */
87: TaoComputeHessian(tao,x,&A,&A,&flg);
88: if (!i) {MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&B);}
89: TaoDefaultComputeHessian(tao,x,&B,&B,&flg,tao->user_hessP);
90: if (fd->complete_print) {
91: MPI_Comm bcomm;
92: PetscViewer viewer;
93: PetscPrintf(comm,"Finite difference Hessian\n");
94: PetscObjectGetComm((PetscObject)B,&bcomm);
95: PetscViewerASCIIGetStdout(bcomm,&viewer);
96: MatView(B,viewer);
97: }
98: /* compare */
99: MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
100: MatNorm(B,NORM_FROBENIUS,&nrm);
101: MatNorm(A,NORM_FROBENIUS,&gnorm);
102: if (fd->complete_print) {
103: MPI_Comm hcomm;
104: PetscViewer viewer;
105: PetscPrintf(comm,"Hand-coded Hessian\n");
106: PetscObjectGetComm((PetscObject)B,&hcomm);
107: PetscViewerASCIIGetStdout(hcomm,&viewer);
108: MatView(A,viewer);
109: PetscPrintf(comm,"Hand-coded minus finite difference Hessian\n");
110: MatView(B,viewer);
111: }
112: if (!gnorm) gnorm = 1.0e-20;
113: PetscPrintf(comm,"ratio ||fd-hc||/||hc|| = %G, difference ||fd-hc|| = %G\n",nrm/gnorm,nrm);
114: }
116: MatDestroy(&B);
117: }
118: tao->reason = TAO_CONVERGED_USER;
119: return(0);
120: }
121: /* ------------------------------------------------------------ */
124: PetscErrorCode TaoDestroy_FD(TaoSolver tao)
125: {
128: PetscFree(tao->data);
129: return(0);
130: }
134: static PetscErrorCode TaoSetFromOptions_FD(TaoSolver tao)
135: {
136: FD_Test *fd = (FD_Test *)tao->data;
141: PetscOptionsHead("Hand-coded Hessian tester options");
142: PetscOptionsBool("-tao_fd_test_display","Display difference between hand-coded and finite difference Hessians","None",fd->complete_print,&fd->complete_print,PETSC_NULL);
143: PetscOptionsBool("-tao_fd_test_gradient","Test Hand-coded gradient against finite-difference gradient","None",fd->check_gradient,&fd->check_gradient,PETSC_NULL);
144: PetscOptionsBool("-tao_fd_test_hessian","Test Hand-coded hessian against finite-difference hessian","None",fd->check_hessian,&fd->check_hessian,PETSC_NULL);
145: if (fd->check_gradient == PETSC_FALSE && fd->check_hessian == PETSC_FALSE) {
146: fd->check_gradient = PETSC_TRUE;
147: }
148: PetscOptionsTail();
149:
150: return(0);
151: }
153: /* ------------------------------------------------------------ */
154: /*C
155: FD_TEST - Test hand-coded Hessian against finite difference Hessian
157: Options Database:
158: . -tao_fd_test_display Display difference between approximate and hand-coded Hessian
160: Level: intermediate
162: .seealso: TaoCreate(), TaoSetType()
164: */
168: PetscErrorCode TaoCreate_FD(TaoSolver tao)
169: {
170: FD_Test *fd;
174: tao->ops->setup = 0;
175: tao->ops->solve = TaoSolve_FD;
176: tao->ops->destroy = TaoDestroy_FD;
177: tao->ops->setfromoptions = TaoSetFromOptions_FD;
178: tao->ops->view = 0;
180:
181: ierr = PetscNewLog(tao,FD_Test,&fd);
182: tao->data = (void*)fd;
183: fd->complete_print = PETSC_FALSE;
184: fd->check_gradient = PETSC_TRUE;
185: fd->check_hessian = PETSC_FALSE;
186: return(0);
187: }