Actual source code: ssfls.c

  1: #include "ssls.h"

  5: PetscErrorCode TaoSetUp_SSFLS(TaoSolver tao)
  6: {
  7:   TAO_SSLS *ssls = (TAO_SSLS *)tao->data;


 12:   VecDuplicate(tao->solution,&tao->gradient); 
 13:   VecDuplicate(tao->solution,&tao->stepdirection); 
 14:   VecDuplicate(tao->solution,&ssls->w); 
 15:   VecDuplicate(tao->solution,&ssls->ff); 
 16:   VecDuplicate(tao->solution,&ssls->dpsi); 
 17:   VecDuplicate(tao->solution,&ssls->da); 
 18:   VecDuplicate(tao->solution,&ssls->db); 
 19:   VecDuplicate(tao->solution,&ssls->t1); 
 20:   VecDuplicate(tao->solution,&ssls->t2); 
 21:   if (!tao->XL) {
 22:     VecDuplicate(tao->solution,&tao->XL); 
 23:   }
 24:   if (!tao->XU) {
 25:     VecDuplicate(tao->solution,&tao->XU); 
 26:   }
 27:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU); 
 28:   return(0);
 29: }


 34: static PetscErrorCode TaoSolve_SSFLS(TaoSolver tao)
 35: {
 36:   TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
 37:   PetscReal psi, ndpsi, normd, innerd, t=0;
 38:   PetscReal delta, rho;
 39:   PetscInt iter=0,kspits;
 40:   TaoSolverTerminationReason reason;
 41:   TaoLineSearchTerminationReason ls_reason;


 46:   /* Assume that Setup has been called!
 47:      Set the structure for the Jacobian and create a linear solver. */
 48:  
 49:   delta = ssls->delta;
 50:   rho = ssls->rho;

 52:   TaoComputeVariableBounds(tao); 
 53:   /* Project solution inside bounds */
 54:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution); 
 55:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);
 56:   
 57:   TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao); 

 59:   /* Calculate the function value and fischer function value at the 
 60:      current iterate */
 61:   TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi); 
 62:   VecNorm(ssls->dpsi,NORM_2,&ndpsi); 


 65:   while (1) {
 66:     ierr=PetscInfo3(tao, "iter: %D, merit: %G, ndpsi: %G\n",
 67:                        iter, ssls->merit, ndpsi); 
 68:     /* Check the termination criteria */
 69:     TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason); 
 70:            
 71:     if (reason!=TAO_CONTINUE_ITERATING) break;

 73:     /* Calculate direction.  (Really negative of newton direction.  Therefore,
 74:        rest of the code uses -d.) */
 75:     KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre,ssls->matflag); 
 76:     KSPSolve(tao->ksp,ssls->ff,tao->stepdirection); 
 77:     KSPGetIterationNumber(tao->ksp,&kspits);  
 78:     tao->ksp_its+=kspits;
 79:     
 80:     VecCopy(tao->stepdirection,ssls->w); 
 81:     VecScale(ssls->w,-1.0); 
 82:     VecBoundGradientProjection(ssls->w,tao->solution,tao->XL,tao->XU,ssls->w); 
 83:     
 84:     VecNorm(ssls->w,NORM_2,&normd); 
 85:     VecDot(ssls->w,ssls->dpsi,&innerd); 

 87:     /* Make sure that we have a descent direction */
 88:     if (innerd >= -delta*pow(normd, rho)) {
 89:       PetscInfo(tao, "newton direction not descent\n"); 
 90:       VecCopy(ssls->dpsi,tao->stepdirection); 
 91:       VecDot(ssls->w,ssls->dpsi,&innerd); 
 92:     }

 94:     VecScale(tao->stepdirection, -1.0); 
 95:     innerd = -innerd;

 97:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
 98:     TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason); 
 99:     VecNorm(ssls->dpsi,NORM_2,&ndpsi); 
100:   }

102:   return(0);
103: }

107: PetscErrorCode TaoDestroy_SSFLS(TaoSolver tao)
108: {
109:   TAO_SSLS *ssls = (TAO_SSLS *)tao->data;


114:   VecDestroy(&ssls->ff); 
115:   VecDestroy(&ssls->w); 
116:   VecDestroy(&ssls->dpsi); 
117:   VecDestroy(&ssls->da); 
118:   VecDestroy(&ssls->db); 
119:   VecDestroy(&ssls->t1); 
120:   VecDestroy(&ssls->t2); 
121:   PetscFree(tao->data); 
122:   tao->data = PETSC_NULL;
123:   return(0);
124: }

126: /* ---------------------------------------------------------- */
130: PetscErrorCode TaoCreate_SSFLS(TaoSolver tao)
131: {
132:   TAO_SSLS *ssls;
134:   const char *armijo_type = TAOLINESEARCH_ARMIJO;

137:   PetscNewLog(tao,TAO_SSLS,&ssls); 
138:   tao->data = (void*)ssls;
139:   tao->ops->solve=TaoSolve_SSFLS;
140:   tao->ops->setup=TaoSetUp_SSFLS;
141:   tao->ops->view=TaoView_SSLS;
142:   tao->ops->setfromoptions=TaoSetFromOptions_SSLS;
143:   tao->ops->destroy = TaoDestroy_SSFLS;

145:   ssls->delta = 1e-10;
146:   ssls->rho = 2.1;

148:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch); 
149:   TaoLineSearchSetType(tao->linesearch,armijo_type); 
150:   TaoLineSearchSetFromOptions(tao->linesearch);
151:   /* Linesearch objective and objectivegradient routines are
152:      set in solve routine */
153:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp); 
154:   
155:   tao->max_it = 2000;
156:   tao->max_funcs = 4000;
157:   tao->fatol = 0; tao->frtol = 0; tao->grtol=0; tao->grtol=0;
158:   tao->gatol = 1.0e-16;
159:   tao->fmin = 1.0e-8;

161:   return(0);
162: }