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: }