Actual source code: bnls.c

  1: /*$Id: s.bnls.c 1.128 02/08/16 18:01:18-05:00 benson@rockies.mcs.anl.gov $*/

  3: #include "bnls.h"       /*I "tao_solver.h" I*/



  9: static int TaoSolve_BNLS(TAO_SOLVER tao, void*solver){

 11:   TAO_BNLS *bnls = (TAO_BNLS *)solver;
 12:   int info;
 13:   TaoInt lsflag,iter=0;
 14:   TaoTerminateReason reason=TAO_CONTINUE_ITERATING;
 15:   double f,f_full,gnorm,gdx,stepsize=1.0;
 16:   TaoTruth success;
 17:   TaoVec *XU, *XL;
 18:   TaoVec *X,  *G=bnls->G, *PG=bnls->PG;
 19:   TaoVec *R=bnls->R, *DXFree=bnls->DXFree;
 20:   TaoVec *DX=bnls->DX, *Work=bnls->Work;
 21:   TaoMat *H, *Hsub=bnls->Hsub;
 22:   TaoIndexSet *FreeVariables = bnls->FreeVariables;

 24:   TaoFunctionBegin;

 26:   /* Check if upper bound greater than lower bound. */
 27:   info = TaoGetSolution(tao,&X);CHKERRQ(info); bnls->X=X;
 28:   info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
 29:   info = TaoEvaluateVariableBounds(tao,XL,XU); CHKERRQ(info);
 30:   info = TaoGetHessian(tao,&H);CHKERRQ(info); bnls->H=H;

 32:   /*   Project the current point onto the feasible set */
 33:   info = X->Median(XL,X,XU); CHKERRQ(info);
 34:   
 35:   info = TaoComputeMeritFunctionGradient(tao,X,&f,G);CHKERRQ(info);
 36:   
 37:   while (reason==TAO_CONTINUE_ITERATING){
 38:     
 39:     /* Project the gradient and calculate the norm */
 40:     info = PG->BoundGradientProjection(G,XL,X,XU);CHKERRQ(info);
 41:     info = PG->Norm2(&gnorm); CHKERRQ(info);
 42:     
 43:     info = TaoMonitor(tao,iter++,f,gnorm,0.0,stepsize,&reason);
 44:     CHKERRQ(info);
 45:     if (reason!=TAO_CONTINUE_ITERATING) break;

 47:     info = FreeVariables->WhichEqual(PG,G); CHKERRQ(info);

 49:     info = TaoComputeHessian(tao,X,H);CHKERRQ(info);
 50:     
 51:     /* Create a reduced linear system */

 53:     info = R->SetReducedVec(G,FreeVariables);CHKERRQ(info);
 54:     info = R->Negate();CHKERRQ(info);

 56:     info = DXFree->SetReducedVec(DX,FreeVariables);CHKERRQ(info);
 57:     info = DXFree->SetToZero(); CHKERRQ(info);
 58:     
 59:     info = Hsub->SetReducedMatrix(H,FreeVariables,FreeVariables);CHKERRQ(info);

 61:     bnls->gamma_factor /= 2;
 62:     success = TAO_FALSE;

 64:     while (success==TAO_FALSE) {
 65:       
 66:       /* Approximately solve the reduced linear system */
 67:       info = TaoPreLinearSolve(tao,Hsub);CHKERRQ(info);
 68:       info = TaoLinearSolve(tao,Hsub,R,DXFree,&success);CHKERRQ(info);

 70:       info = DX->SetToZero(); CHKERRQ(info);
 71:       info = DX->ReducedXPY(DXFree,FreeVariables);CHKERRQ(info);
 72:       info = DX->Dot(G,&gdx); CHKERRQ(info);

 74:       if (gdx>=0 || success==TAO_FALSE) { /* Modify diagonal of Hessian if not a descent direction */
 75:         bnls->gamma_factor *= 2; 
 76:         bnls->gamma = bnls->gamma_factor*(gnorm); 
 77: #if !defined(PETSC_USE_COMPLEX)
 78:         info=PetscInfo2(tao,"TaoSolve_NLS:  modify diagonal (assume same nonzero structure), gamma_factor=%g, gamma=%g\n",bnls->gamma_factor,bnls->gamma);
 79:         CHKERRQ(info);
 80: #else
 81:         info=PetscInfo3(tao,"TaoSolve_NLS:  modify diagonal (asuume same nonzero structure), gamma_factor=%g, gamma=%g\n",
 82:              bnls->gamma_factor,PetscReal(bnls->gamma));CHKERRQ(info);
 83: #endif
 84:         info = Hsub->ShiftDiagonal(bnls->gamma);CHKERRQ(info);
 85:         success = TAO_FALSE;
 86:         
 87:       } else {
 88:         success = TAO_TRUE;
 89:       }

 91:     }
 92:     
 93:     stepsize=1.0;        
 94:     info = TaoLineSearchApply(tao,X,G,DX,Work,
 95:                               &f,&f_full,&stepsize,&lsflag);
 96:     CHKERRQ(info);

 98:     
 99:   }  /* END MAIN LOOP  */

101:   TaoFunctionReturn(0);
102: }


105: /*------------------------------------------------------------*/
108: static int TaoSetDown_BNLS(TAO_SOLVER tao, void*solver)
109: {
110:   TAO_BNLS *bnls = (TAO_BNLS *)solver;
111:   int      info;
112:   /* Free allocated memory in BNLS structure */
113:   TaoFunctionBegin;
114:   
115:   info = TaoVecDestroy(bnls->DX);CHKERRQ(info);bnls->DX=0;
116:   info = TaoVecDestroy(bnls->Work);CHKERRQ(info);
117:   info = TaoVecDestroy(bnls->DXFree);CHKERRQ(info);
118:   info = TaoVecDestroy(bnls->R);CHKERRQ(info);
119:   info = TaoVecDestroy(bnls->G);CHKERRQ(info);
120:   info = TaoVecDestroy(bnls->PG);CHKERRQ(info);
121:   info = TaoVecDestroy(bnls->XL);CHKERRQ(info);
122:   info = TaoVecDestroy(bnls->XU);CHKERRQ(info);
123:   
124:   info = TaoIndexSetDestroy(bnls->FreeVariables);CHKERRQ(info);
125:   info = TaoMatDestroy(bnls->Hsub);CHKERRQ(info);
126:   info = TaoDestroyLinearSolver(tao);CHKERRQ(info);

128:   TaoFunctionReturn(0);
129: }

131: /*------------------------------------------------------------*/
134: static int TaoSetOptions_BNLS(TAO_SOLVER tao, void*solver)
135: {
136:   int        info;
137:   TaoInt     ival;
138:   TaoTruth flg;

140:   TaoFunctionBegin;

142:   info = TaoOptionsHead("Newton Line Search Method for bound constrained optimization");CHKERRQ(info);

144:   info = TaoOptionInt("-redistribute","Redistribute Free variables (> 1 processors, only)","TaoPetscISType",1,&ival,&flg); CHKERRQ(info);

146:   info = TaoOptionName("-submatrixfree","Mask full matrix instead of extract submatrices","TaoPetscISType",&flg); CHKERRQ(info);

148:   info = TaoOptionsTail();CHKERRQ(info);
149:   info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);

151:   TaoFunctionReturn(0);
152: }

154: /*------------------------------------------------------------*/
157: static int TaoView_BNLS(TAO_SOLVER tao,void*solver)
158: {
159:   int        info;

161:   TaoFunctionBegin;
162:   info = TaoLineSearchView(tao);CHKERRQ(info);
163:   TaoFunctionReturn(0);
164: }


167: /* ---------------------------------------------------------- */
170: static int TaoSetUp_BNLS(TAO_SOLVER tao, void*solver){

172:   int info;
173:   TAO_BNLS *bnls = (TAO_BNLS *)solver;
174:   TaoVec* X;
175:   TaoMat *HH;

177:   TaoFunctionBegin;
178:   info = TaoGetSolution(tao,&bnls->X);CHKERRQ(info); X=bnls->X;
179:   info = TaoGetHessian(tao,&bnls->H);CHKERRQ(info);  HH=bnls->H;

181:   /* Allocate some arrays */
182:   info = X->Clone(&bnls->DX); CHKERRQ(info);
183:   info = X->Clone(&bnls->Work); CHKERRQ(info);
184:   info = X->Clone(&bnls->DXFree); CHKERRQ(info);
185:   info = X->Clone(&bnls->R); CHKERRQ(info);
186:   info = X->Clone(&bnls->G); CHKERRQ(info);
187:   info = X->Clone(&bnls->PG); CHKERRQ(info);
188:   info = X->Clone(&bnls->XL); CHKERRQ(info);
189:   info = X->Clone(&bnls->XU); CHKERRQ(info);

191:   info = TaoSetLagrangianGradientVector(tao,bnls->PG);CHKERRQ(info);
192:   info = TaoSetStepDirectionVector(tao,bnls->DX);CHKERRQ(info);
193:   info = TaoSetVariableBounds(tao,bnls->XL,bnls->XU);CHKERRQ(info);

195:   info = X->CreateIndexSet(&bnls->FreeVariables); CHKERRQ(info);
196:   info = bnls->H->CreateReducedMatrix(bnls->FreeVariables,bnls->FreeVariables,&bnls->Hsub); CHKERRQ(info);

198:   info = TaoCreateLinearSolver(tao,HH,100,0); CHKERRQ(info);

200:   info = TaoCheckFGH(tao);CHKERRQ(info);

202:   TaoFunctionReturn(0);
203: }

205: /*------------------------------------------------------------*/
208: int TaoGetDualVariables_BNLS(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void* solver)
209: {

211:   TAO_BNLS *bnls = (TAO_BNLS *) solver;
212:   TaoVec  *G=bnls->G,*GP=bnls->Work;
213:   int       info;

215:   TaoFunctionBegin;

217:   info = DXL->Waxpby(-1,G,1.0,GP); CHKERRQ(info);
218:   info = DXU->SetToZero(); CHKERRQ(info);
219:   info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);

221:   info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
222:   info = DXU->Axpy(1.0,DXL); CHKERRQ(info);

224:   TaoFunctionReturn(0);
225: }


228: /*------------------------------------------------------------*/
232: int TaoCreate_BNLS(TAO_SOLVER tao)
233: {
234:   TAO_BNLS *bnls;
235:   int      info;

237:   TaoFunctionBegin;

239:   info = TaoNew(TAO_BNLS,&bnls); CHKERRQ(info);
240:   info = PetscLogObjectMemory(tao,sizeof(TAO_BNLS)); CHKERRQ(info);

242:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_BNLS,(void*)bnls); CHKERRQ(info);
243:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_BNLS,TaoSetDown_BNLS); CHKERRQ(info);
244:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_BNLS); CHKERRQ(info);
245:   info=TaoSetTaoViewRoutine(tao,TaoView_BNLS); CHKERRQ(info);
246:   info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_BNLS); CHKERRQ(info);

248:   info = TaoSetMaximumIterates(tao,500); CHKERRQ(info);
249:   info = TaoSetTolerances(tao,1e-12,1e-12,0,0); CHKERRQ(info);

251:   /* Initialize pointers and variables */

253:   bnls->gamma = 0.0;
254:   bnls->gamma_factor = 0.01;
255:   bnls->DX=0;
256:   bnls->DXFree=0;
257:   bnls->R=0;
258:   bnls->Work=0;
259:   bnls->FreeVariables=0;
260:   bnls->Hsub=0;

262:   info = TaoCreateMoreThuenteBoundLineSearch(tao,0,0.9); CHKERRQ(info);

264:   TaoFunctionReturn(0);
265: }