Actual source code: fdtest.c
petsc-3.7.3 2016-08-01
1: #include <petsc/private/taoimpl.h> /*I "petsctao.h" I*/
3: typedef struct {
4: PetscBool check_gradient;
5: PetscBool check_hessian;
6: PetscBool complete_print;
7: } Tao_Test;
9: /*
10: TaoSolve_Test - Tests whether a hand computed Hessian
11: matches one compute via finite differences.
12: */
15: PetscErrorCode TaoSolve_Test(Tao tao)
16: {
17: Mat A = tao->hessian,B;
18: Vec x = tao->solution,g1,g2;
20: PetscInt i;
21: PetscReal nrm,gnorm,hcnorm,fdnorm;
22: MPI_Comm comm;
23: Tao_Test *fd = (Tao_Test*)tao->data;
26: comm = ((PetscObject)tao)->comm;
27: if (fd->check_gradient) {
28: VecDuplicate(x,&g1);
29: VecDuplicate(x,&g2);
31: PetscPrintf(comm,"Testing hand-coded gradient (hc) against finite difference gradient (fd), if the ratio ||fd - hc|| / ||hc|| is\n");
32: PetscPrintf(comm,"0 (1.e-8), the hand-coded gradient is probably correct.\n");
34: if (!fd->complete_print) {
35: PetscPrintf(comm,"Run with -tao_test_display to show difference\n");
36: PetscPrintf(comm,"between hand-coded and finite difference gradient.\n");
37: }
38: for (i=0; i<3; i++) {
39: if (i == 1) {VecSet(x,-1.0);}
40: else if (i == 2) {VecSet(x,1.0);}
42: /* Compute both version of gradient */
43: TaoComputeGradient(tao,x,g1);
44: TaoDefaultComputeGradient(tao,x,g2,NULL);
45: if (fd->complete_print) {
46: MPI_Comm gcomm;
47: PetscViewer viewer;
48: PetscPrintf(comm,"Finite difference gradient\n");
49: PetscObjectGetComm((PetscObject)g2,&gcomm);
50: PetscViewerASCIIGetStdout(gcomm,&viewer);
51: VecView(g2,viewer);
52: PetscPrintf(comm,"Hand-coded gradient\n");
53: PetscObjectGetComm((PetscObject)g1,&gcomm);
54: PetscViewerASCIIGetStdout(gcomm,&viewer);
55: VecView(g1,viewer);
56: PetscPrintf(comm,"\n");
57: }
59: VecAXPY(g2,-1.0,g1);
60: VecNorm(g1,NORM_2,&hcnorm);
61: VecNorm(g2,NORM_2,&fdnorm);
63: if (!hcnorm) hcnorm=1.0e-20;
64: PetscPrintf(comm,"ratio ||fd-hc||/||hc|| = %g, difference ||fd-hc|| = %g\n", (double)(fdnorm/hcnorm), (double)fdnorm);
66: }
67: VecDestroy(&g1);
68: VecDestroy(&g2);
69: }
71: if (fd->check_hessian) {
72: if (A != tao->hessian_pre) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Cannot test with alternative preconditioner");
74: PetscPrintf(comm,"Testing hand-coded Hessian (hc) against finite difference Hessian (fd). If the ratio is\n");
75: PetscPrintf(comm,"O (1.e-8), the hand-coded Hessian is probably correct.\n");
77: if (!fd->complete_print) {
78: PetscPrintf(comm,"Run with -tao_test_display to show difference\n");
79: PetscPrintf(comm,"of hand-coded and finite difference Hessian.\n");
80: }
81: for (i=0;i<3;i++) {
82: /* compute both versions of Hessian */
83: TaoComputeHessian(tao,x,A,A);
84: if (!i) {MatConvert(A,MATSAME,MAT_INITIAL_MATRIX,&B);}
85: TaoDefaultComputeHessian(tao,x,B,B,tao->user_hessP);
86: if (fd->complete_print) {
87: MPI_Comm bcomm;
88: PetscViewer viewer;
89: PetscPrintf(comm,"Finite difference Hessian\n");
90: PetscObjectGetComm((PetscObject)B,&bcomm);
91: PetscViewerASCIIGetStdout(bcomm,&viewer);
92: MatView(B,viewer);
93: }
94: /* compare */
95: MatAYPX(B,-1.0,A,DIFFERENT_NONZERO_PATTERN);
96: MatNorm(B,NORM_FROBENIUS,&nrm);
97: MatNorm(A,NORM_FROBENIUS,&gnorm);
98: if (fd->complete_print) {
99: MPI_Comm hcomm;
100: PetscViewer viewer;
101: PetscPrintf(comm,"Hand-coded Hessian\n");
102: PetscObjectGetComm((PetscObject)B,&hcomm);
103: PetscViewerASCIIGetStdout(hcomm,&viewer);
104: MatView(A,viewer);
105: PetscPrintf(comm,"Hand-coded minus finite difference Hessian\n");
106: MatView(B,viewer);
107: }
108: if (!gnorm) gnorm = 1.0e-20;
109: PetscPrintf(comm,"ratio ||fd-hc||/||hc|| = %g, difference ||fd-hc|| = %g\n",(double)(nrm/gnorm),(double)nrm);
110: }
112: MatDestroy(&B);
113: }
114: tao->reason = TAO_CONVERGED_USER;
115: return(0);
116: }
117: /* ------------------------------------------------------------ */
120: PetscErrorCode TaoDestroy_Test(Tao tao)
121: {
125: PetscFree(tao->data);
126: return(0);
127: }
131: static PetscErrorCode TaoSetFromOptions_Test(PetscOptionItems *PetscOptionsObject,Tao tao)
132: {
133: Tao_Test *fd = (Tao_Test *)tao->data;
137: PetscOptionsHead(PetscOptionsObject,"Hand-coded Hessian tester options");
138: PetscOptionsBool("-tao_test_display","Display difference between hand-coded and finite difference Hessians","None",fd->complete_print,&fd->complete_print,NULL);
139: PetscOptionsBool("-tao_test_gradient","Test Hand-coded gradient against finite-difference gradient","None",fd->check_gradient,&fd->check_gradient,NULL);
140: PetscOptionsBool("-tao_test_hessian","Test Hand-coded hessian against finite-difference hessian","None",fd->check_hessian,&fd->check_hessian,NULL);
141: if (fd->check_gradient == PETSC_FALSE && fd->check_hessian == PETSC_FALSE) {
142: fd->check_gradient = PETSC_TRUE;
143: }
144: PetscOptionsTail();
145: return(0);
146: }
148: /* ------------------------------------------------------------ */
149: /*C
150: FD_TEST - Test hand-coded Hessian against finite difference Hessian
152: Options Database:
153: . -tao_test_display Display difference between approximate and hand-coded Hessian
155: Level: intermediate
157: .seealso: TaoCreate(), TaoSetType()
159: */
162: PETSC_EXTERN PetscErrorCode TaoCreate_Test(Tao tao)
163: {
164: Tao_Test *fd;
168: tao->ops->setup = 0;
169: tao->ops->solve = TaoSolve_Test;
170: tao->ops->destroy = TaoDestroy_Test;
171: tao->ops->setfromoptions = TaoSetFromOptions_Test;
172: tao->ops->view = 0;
173: PetscNewLog(tao,&fd);
174: tao->data = (void*)fd;
175: fd->complete_print = PETSC_FALSE;
176: fd->check_gradient = PETSC_TRUE;
177: fd->check_hessian = PETSC_FALSE;
178: return(0);
179: }