Actual source code: ssils.c

  1: #include "ssls.h"

  5: PetscErrorCode TaoSetUp_SSILS(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->ff); 
 15:   VecDuplicate(tao->solution,&ssls->dpsi); 
 16:   VecDuplicate(tao->solution,&ssls->da); 
 17:   VecDuplicate(tao->solution,&ssls->db); 
 18:   VecDuplicate(tao->solution,&ssls->t1); 
 19:   VecDuplicate(tao->solution,&ssls->t2); 
 20:   return(0);
 21: }


 26: PetscErrorCode TaoDestroy_SSILS(TaoSolver tao)
 27: {
 28:   TAO_SSLS *ssls = (TAO_SSLS *)tao->data;


 33:   VecDestroy(&ssls->ff); 
 34:   VecDestroy(&ssls->dpsi); 
 35:   VecDestroy(&ssls->da); 
 36:   VecDestroy(&ssls->db); 
 37:   VecDestroy(&ssls->t1); 
 38:   VecDestroy(&ssls->t2); 
 39:   PetscFree(tao->data); 
 40:   tao->data = PETSC_NULL;
 41:   return(0);
 42: }

 46: static PetscErrorCode TaoSolve_SSILS(TaoSolver tao)
 47: {
 48:   TAO_SSLS *ssls = (TAO_SSLS *)tao->data;
 49:   PetscReal psi, ndpsi, normd, innerd, t=0;
 50:   PetscReal delta, rho;
 51:   PetscInt iter=0,kspits;
 52:   TaoSolverTerminationReason reason;
 53:   TaoLineSearchTerminationReason ls_reason;


 58:   /* Assume that Setup has been called!
 59:      Set the structure for the Jacobian and create a linear solver. */
 60:  
 61:   delta = ssls->delta;
 62:   rho = ssls->rho;

 64:   TaoComputeVariableBounds(tao); 
 65:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution); 
 66:   TaoLineSearchSetObjectiveAndGradientRoutine(tao->linesearch,Tao_SSLS_FunctionGradient,tao);
 67:   
 68:   TaoLineSearchSetObjectiveRoutine(tao->linesearch,Tao_SSLS_Function,tao); 

 70:   /* Calculate the function value and fischer function value at the 
 71:      current iterate */
 72:   TaoLineSearchComputeObjectiveAndGradient(tao->linesearch,tao->solution,&psi,ssls->dpsi); 
 73:   VecNorm(ssls->dpsi,NORM_2,&ndpsi); 


 76:   while (1) {
 77:     ierr=PetscInfo3(tao, "iter: %D, merit: %G, ndpsi: %G\n",
 78:                        iter, ssls->merit, ndpsi); 
 79:     /* Check the termination criteria */
 80:     TaoMonitor(tao,iter++,ssls->merit,ndpsi,0.0,t,&reason); 
 81:            
 82:     if (reason!=TAO_CONTINUE_ITERATING) break;

 84:     /* Calculate direction.  (Really negative of newton direction.  Therefore,
 85:        rest of the code uses -d.) */
 86:     KSPSetOperators(tao->ksp,tao->jacobian,tao->jacobian_pre,ssls->matflag); 
 87:     KSPSolve(tao->ksp,ssls->ff,tao->stepdirection); 
 88:     KSPGetIterationNumber(tao->ksp,&kspits);  
 89:     tao->ksp_its+=kspits;
 90:     VecNorm(tao->stepdirection,NORM_2,&normd); 
 91:     VecDot(tao->stepdirection,ssls->dpsi,&innerd); 

 93:     /* Make sure that we have a descent direction */
 94:     if (innerd <= delta*pow(normd, rho)) {
 95:       PetscInfo(tao, "newton direction not descent\n"); 
 96:       VecCopy(ssls->dpsi,tao->stepdirection); 
 97:       VecDot(tao->stepdirection,ssls->dpsi,&innerd); 
 98:     }

100:     VecScale(tao->stepdirection, -1.0); 
101:     innerd = -innerd;

103:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
104:     TaoLineSearchApply(tao->linesearch,tao->solution,&psi,ssls->dpsi,tao->stepdirection,&t,&ls_reason); 
105:     VecNorm(ssls->dpsi,NORM_2,&ndpsi); 
106:   }

108:   return(0);
109: }

111: /* ---------------------------------------------------------- */
115: PetscErrorCode TaoCreate_SSILS(TaoSolver tao)
116: {
117:   TAO_SSLS *ssls;
119:   const char *armijo_type = TAOLINESEARCH_ARMIJO;

122:   PetscNewLog(tao,TAO_SSLS,&ssls); 
123:   tao->data = (void*)ssls;
124:   tao->ops->solve=TaoSolve_SSILS;
125:   tao->ops->setup=TaoSetUp_SSILS;
126:   tao->ops->view=TaoView_SSLS;
127:   tao->ops->setfromoptions=TaoSetFromOptions_SSLS;
128:   tao->ops->destroy = TaoDestroy_SSILS;

130:   ssls->delta = 1e-10;
131:   ssls->rho = 2.1;

133:   TaoLineSearchCreate(((PetscObject)tao)->comm,&tao->linesearch); 
134:   TaoLineSearchSetType(tao->linesearch,armijo_type); 
135:   TaoLineSearchSetFromOptions(tao->linesearch);
136:   /* Note: linesearch objective and objectivegradient routines are set in solve routine */
137:   KSPCreate(((PetscObject)tao)->comm,&tao->ksp); 
138:   
139:   tao->max_it = 2000;
140:   tao->max_funcs = 4000;
141:   tao->fatol = 0; tao->frtol = 0; tao->gttol=0; tao->grtol=0;
142:   tao->gatol = 1.0e-16;
143:   tao->fmin = 1.0e-8;

145:   return(0);
146: }