Actual source code: blmvm.c

  1: #include "taolinesearch.h"
  2: #include "src/matrix/lmvmmat.h"
  3: #include "blmvm.h"

  5: /*------------------------------------------------------------*/
  8: static PetscErrorCode TaoSolve_BLMVM(TaoSolver tao)
  9: {
 11:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;

 13:   TaoSolverTerminationReason reason = TAO_CONTINUE_ITERATING;
 14:   TaoLineSearchTerminationReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;

 16:   PetscReal f, fold, gdx, gnorm;
 17:   PetscReal stepsize = 1.0,delta;

 19:   PetscInt iter = 0;
 20:   

 23:   
 24:   /*  Project initial point onto bounds */
 25:   TaoComputeVariableBounds(tao); 
 26:   VecMedian(tao->XL,tao->solution,tao->XU,tao->solution); 
 27:   TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU); 

 29:   /* Check convergence criteria */
 30:   TaoComputeObjectiveAndGradient(tao, tao->solution,&f,blmP->unprojected_gradient); 
 31:   VecBoundGradientProjection(blmP->unprojected_gradient,tao->solution, tao->XL,tao->XU,tao->gradient); 


 34:   VecNorm(tao->gradient,NORM_2,&gnorm); 
 35:   if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
 36:     SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN");
 37:   }

 39:   TaoMonitor(tao, iter, f, gnorm, 0.0, stepsize, &reason); 
 40:   if (reason != TAO_CONTINUE_ITERATING) {
 41:     return(0);
 42:   }

 44:   /* Set initial scaling for the function */
 45:   if (f != 0.0) {
 46:     delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm);
 47:   }
 48:   else {
 49:     delta = 2.0 / (gnorm*gnorm);
 50:   }
 51:   MatLMVMSetDelta(blmP->M,delta); 

 53:   /* Set counter for gradient/reset steps */
 54:   blmP->grad = 0;
 55:   blmP->reset = 0;

 57:   /* Have not converged; continue with Newton method */
 58:   while (reason == TAO_CONTINUE_ITERATING) {
 59:     
 60:     /* Compute direction */
 61:     MatLMVMUpdate(blmP->M, tao->solution, tao->gradient); 
 62:     MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection); 
 63:     VecBoundGradientProjection(tao->stepdirection,tao->solution,tao->XL,tao->XU,tao->gradient); 

 65:     /* Check for success (descent direction) */
 66:     VecDot(blmP->unprojected_gradient, tao->gradient, &gdx); 
 67:     if (gdx <= 0) {
 68:       /* Step is not descent or solve was not successful
 69:          Use steepest descent direction (scaled) */
 70:       ++blmP->grad;

 72:       if (f != 0.0) {
 73:         delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm);
 74:       }
 75:       else {
 76:         delta = 2.0 / (gnorm*gnorm);
 77:       }
 78:       MatLMVMSetDelta(blmP->M,delta); 
 79:       MatLMVMReset(blmP->M); 
 80:       MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient); 
 81:       MatLMVMSolve(blmP->M,blmP->unprojected_gradient, tao->stepdirection); 
 82:     } 
 83:     VecScale(tao->stepdirection,-1.0); 

 85:     /* Perform the linesearch */
 86:     fold = f;
 87:     VecCopy(tao->solution, blmP->Xold); 
 88:     VecCopy(blmP->unprojected_gradient, blmP->Gold); 
 89:     TaoLineSearchSetInitialStepLength(tao->linesearch,1.0); 
 90:     TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status); 
 91:     TaoAddLineSearchCounts(tao); 

 93:     if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
 94:       /* Linesearch failed
 95:          Reset factors and use scaled (projected) gradient step */
 96:       ++blmP->reset;

 98:       f = fold;
 99:       VecCopy(blmP->Xold, tao->solution); 
100:       VecCopy(blmP->Gold, blmP->unprojected_gradient); 

102:       if (f != 0.0) {
103:         delta = 2.0* PetscAbsScalar(f) / (gnorm*gnorm);
104:       }
105:       else {
106:         delta = 2.0/ (gnorm*gnorm);
107:       }
108:       MatLMVMSetDelta(blmP->M,delta); 
109:       MatLMVMReset(blmP->M); 
110:       MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient); 
111:       MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection); 
112:       VecScale(tao->stepdirection, -1.0); 

114:       /* This may be incorrect; linesearch has values fo stepmax and stepmin
115:          that should be reset. */
116:       TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
117:       TaoLineSearchApply(tao->linesearch,tao->solution,&f, blmP->unprojected_gradient, tao->stepdirection,  &stepsize, &ls_status); 
118:       TaoAddLineSearchCounts(tao); 

120:       if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
121:         tao->reason = TAO_DIVERGED_LS_FAILURE;
122:         break;
123:       }
124:     }

126:     /* Check for termination */
127:     VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient); 
128:     VecNorm(tao->gradient, NORM_2, &gnorm); 


131:     if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) {
132:       SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Not-a-Number");
133:     }
134:     iter++;
135:     TaoMonitor(tao, iter, f, gnorm, 0.0, stepsize, &reason); 
136:   }



140:   return(0);
141: }

145: static PetscErrorCode TaoSetup_BLMVM(TaoSolver tao)
146: {
147:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
148:   PetscInt n,N;

152:   /* Existence of tao->solution checked in TaoSetup() */
153:   
154:   VecDuplicate(tao->solution,&blmP->Xold); 
155:   VecDuplicate(tao->solution,&blmP->Gold); 
156:   VecDuplicate(tao->solution, &blmP->unprojected_gradient);

158:   if (!tao->stepdirection) {
159:       VecDuplicate(tao->solution, &tao->stepdirection);
160:       
161:   }
162:   if (!tao->gradient) {
163:       VecDuplicate(tao->solution,&tao->gradient); 
164:   }
165:   if (!tao->XL) {
166:       VecDuplicate(tao->solution,&tao->XL); 
167:       VecSet(tao->XL,TAO_NINFINITY); 
168:   }
169:   if (!tao->XU) {
170:       VecDuplicate(tao->solution,&tao->XU); 
171:       VecSet(tao->XU,TAO_INFINITY); 
172:   }
173:   /* Create matrix for the limited memory approximation */
174:   VecGetLocalSize(tao->solution,&n); 
175:   VecGetSize(tao->solution,&N); 
176:   MatCreateLMVM(((PetscObject)tao)->comm,n,N,&blmP->M); 
177:   MatLMVMAllocateVectors(blmP->M,tao->solution); 

179:    return(0);
180: }

182: /* ---------------------------------------------------------- */
185: static PetscErrorCode TaoDestroy_BLMVM(TaoSolver tao)
186: {
187:   TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;

191:   if (tao->setupcalled) {
192:     MatDestroy(&blmP->M); 
193:     VecDestroy(&blmP->unprojected_gradient); 
194:     VecDestroy(&blmP->Xold); 
195:     VecDestroy(&blmP->Gold); 
196:   }
197:   PetscFree(tao->data); 
198:   tao->data = PETSC_NULL;

200:   return(0);
201: }

203: /*------------------------------------------------------------*/
206: static PetscErrorCode TaoSetFromOptions_BLMVM(TaoSolver tao)
207: {


212:   PetscOptionsHead("Limited-memory variable-metric method for bound constrained optimization"); 
213:   TaoLineSearchSetFromOptions(tao->linesearch);
214:   PetscOptionsTail();
215:   return(0);
216: }


219: /*------------------------------------------------------------*/
222: static int TaoView_BLMVM(TaoSolver tao, PetscViewer viewer)
223: {

225:     
226:     TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
227:     PetscBool isascii;


232:     PetscTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii); 
233:     if (isascii) {
234:       PetscViewerASCIIPushTab(viewer); 
235:       PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lmP->grad); 
236:       PetscViewerASCIIPopTab(viewer); 
237:     } else {
238:       SETERRQ1(((PetscObject)tao)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for TAO BLMVM",((PetscObject)viewer)->type_name);
239:     }
240:     return(0);
241: }

245: static PetscErrorCode TaoComputeDual_BLMVM(TaoSolver tao, Vec DXL, Vec DXU)
246: {
247:   TAO_BLMVM *blm = (TAO_BLMVM *) tao->data;


255:   if (!tao->gradient || !blm->unprojected_gradient) {
256:       SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");
257:   }
258:   
259:   VecCopy(tao->gradient,DXL); 
260:   VecAXPY(DXL,-1.0,blm->unprojected_gradient); 
261:   VecSet(DXU,0.0); 
262:   VecPointwiseMax(DXL,DXL,DXU); 

264:   VecCopy(blm->unprojected_gradient,DXU); 
265:   VecAXPY(DXU,-1.0,tao->gradient); 
266:   VecAXPY(DXU,1.0,DXL); 

268:   return(0);
269: }

271: /* ---------------------------------------------------------- */
275: PetscErrorCode TaoCreate_BLMVM(TaoSolver tao)
276: {
277:   TAO_BLMVM *blmP;
278:   const char *morethuente_type = TAOLINESEARCH_MT;

282:   tao->ops->setup = TaoSetup_BLMVM;
283:   tao->ops->solve = TaoSolve_BLMVM;
284:   tao->ops->view = TaoView_BLMVM;
285:   tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
286:   tao->ops->destroy = TaoDestroy_BLMVM;
287:   tao->ops->computedual = TaoComputeDual_BLMVM;

289:   PetscNewLog(tao, TAO_BLMVM, &blmP); 
290:   tao->data = (void*)blmP;
291:   tao->max_it = 2000;
292:   tao->max_funcs = 4000;
293:   tao->fatol = 1e-4;
294:   tao->frtol = 1e-4;

296:   TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch); 
297:   TaoLineSearchSetType(tao->linesearch, morethuente_type); 
298:   TaoLineSearchUseTaoSolverRoutines(tao->linesearch,tao); 
299:   return(0);

301: }