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); 

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


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

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

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

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

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

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

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

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

 92:     if (ls_status<0) {
 93:       /* Linesearch failed
 94:          Reset factors and use scaled (projected) gradient step */
 95:       ++blmP->reset;

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

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

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

119:       if ((int) ls_status < 0) {

121:       }
122:     }

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


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



138:   return(0);
139: }

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

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

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

178:    return(0);
179: }

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

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

199:   return(0);
200: }

202: /*------------------------------------------------------------*/
205: static PetscErrorCode TaoSetFromOptions_BLMVM(TaoSolver tao)
206: {


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


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

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


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

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


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

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

267:   return(0);
268: }

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

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

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

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

300: }