Actual source code: rscs.c

  1: #include "src/complementarity/impls/ssls/ssls.h"
  2: // #include "src/tao_impl.h"

  4: int Tao_RSCS_FunctionGradient(TAO_SOLVER, TaoVec*, double*, TaoVec*, void *);
  5: int Tao_RSCS_Function(TAO_SOLVER, TaoVec *, double *, void *);
  6: static int TaoSolve_RSCS(TAO_SOLVER tao, void*solver);

  8: typedef struct {
  9:   TaoVec *f;
 10:   TaoMat *J;
 11: } FGMeritCtx;


 16: static int TaoSolve_RSCS(TAO_SOLVER tao, void*solver){

 18:   TAO_SSLS *asls = (TAO_SSLS *)solver;
 19:   FGMeritCtx  meritctx;
 20:   TaoTerminateReason reason=TAO_CONTINUE_ITERATING;
 21:   TaoVec *x, *l, *u, *ff, *d, *w, *g;
 22:   TaoMat *J,*Jsub;
 23:   TaoVec *dxfree,*r1;
 24:   TaoIndexSet *FreeVariableSet;
 25:   int iter=0, info;
 26:   int lsflag;
 27:   double ndpsi,stepsize=1.0;
 28:   //  double psi,psi1,psi2;
 29:   TaoTruth success;
 30:   double fff,gdx;
 31:   double gamma = 0.0, gamma_factor =1.0e-12;
 32:   //  double fref;

 34:   TaoFunctionBegin;

 36:   ff = asls->dpsi;
 37:   d = asls->d;
 38:   g = asls->t1;
 39:   w = asls->w;

 41:   /*
 42:   f = asls->f;
 43:   ff = asls->ff;
 44:   t2 = asls->t2;
 45:   da = asls->da;
 46:   db = asls->db;
 47:   */
 48:   /* Check if upper bound greater than lower bound. */
 49:   info = TaoGetSolution(tao, &x); CHKERRQ(info);
 50:   info = TaoGetVariableBounds(tao, &l, &u); CHKERRQ(info);
 51:   info = TaoEvaluateVariableBounds(tao,l,u); CHKERRQ(info);
 52:   info = TaoGetJacobian(tao, &J); CHKERRQ(info);
 53:   meritctx.J=J;
 54:   meritctx.f=ff;
 55:   info = TaoSetMeritFunction(tao, Tao_RSCS_Function, Tao_RSCS_FunctionGradient,
 56:                              TAO_NULL, TAO_NULL, (void*)&meritctx ); CHKERRQ(info);

 58:   /*   Project the current point onto the feasible set */
 59:   info = x->Median(l,x,u); CHKERRQ(info);
 60: 
 61:   info = TaoComputeMeritFunctionGradient(tao, x, &fff, g); CHKERRQ(info);

 63:   info = x->CreateIndexSet(&FreeVariableSet); CHKERRQ(info);
 64:   info = J->CreateReducedMatrix(FreeVariableSet,FreeVariableSet,&Jsub); CHKERRQ(info);
 65:   info = x->Clone(&dxfree); CHKERRQ(info);
 66:   info = x->Clone(&r1); CHKERRQ(info);

 68:   while (reason==TAO_CONTINUE_ITERATING){
 69: 
 70:     /* Project the gradient and calculate the norm */
 71:     info = w->BoundGradientProjection(ff,l,x,u);CHKERRQ(info);
 72:     info = w->Norm2(&ndpsi); CHKERRQ(info);

 74:     info = TaoMonitor(tao,iter++,fff,ndpsi,0.0,stepsize,&reason);
 75:     CHKERRQ(info);

 77:     if (reason!=TAO_CONTINUE_ITERATING) break;

 79:     info = FreeVariableSet->WhichEqual(w,ff); CHKERRQ(info);

 81:     /* Create a reduced linear system */
 82:     info = r1->SetReducedVec(ff,FreeVariableSet);CHKERRQ(info);
 83:     info = r1->Negate();CHKERRQ(info);

 85:     info = dxfree->SetReducedVec(d,FreeVariableSet);CHKERRQ(info);
 86:     info = dxfree->SetToZero(); CHKERRQ(info);
 87: 
 88:     info = TaoComputeJacobian(tao,x,J);CHKERRQ(info);
 89:     info = Jsub->SetReducedMatrix(J,FreeVariableSet,FreeVariableSet);CHKERRQ(info);
 90: 
 91:     success = TAO_FALSE;
 92:     gamma = gamma_factor*(ndpsi);
 93:     while (success==TAO_FALSE) {
 94: 
 95:       /* Approximately solve the reduced linear system */
 96:       info = TaoPreLinearSolve(tao,Jsub);CHKERRQ(info);
 97:       info = TaoLinearSolve(tao,Jsub,r1,dxfree,&success);CHKERRQ(info);
 98: 
 99:       info = d->SetToZero(); CHKERRQ(info);
100:       info = d->ReducedXPY(dxfree,FreeVariableSet);CHKERRQ(info);
101: 
102:       info = d->Dot(ff,&gdx); CHKERRQ(info);
103: 

105:       if (success==TAO_FALSE) { /* Modify diagonal of Hessian if not a descent direction */

107:         info = Jsub->SetReducedMatrix(J,FreeVariableSet,FreeVariableSet);CHKERRQ(info);
108:         gamma *=10;
109:         //          printf("Shift diagonal: %4.2e\n",gamma);
110:         info = PetscLogInfo((tao,"TaoSolve_NLS:  modify diagonal (asuume same nonzero structure), gamma_factor=%g, gamma=%g\n",gamma_factor,gamma));CHKERRQ(info);
111:         info = Jsub->ShiftDiagonal(gamma);CHKERRQ(info);
112:         success = TAO_FALSE;
113: 
114:       } else {
115:         success = TAO_TRUE;
116:       }
117: 
118:     }

120:     //    fref=fff;
121:     info = g->ScaleCopyFrom(-1.0,d);CHKERRQ(info);
122:     stepsize=1.0;
123:     info = TaoLineSearchApply(tao,x,g,d,w,
124:                               &fff,&stepsize,&lsflag);
125:     CHKERRQ(info);


128:     if (lsflag!=0){
129:       int attempts = 0;

131:       printf("Try Again \n");
132:       info = FreeVariableSet->WhichBetween(l,x,u); CHKERRQ(info);
133: 
134:       /* Create a reduced linear system */
135:       info = r1->SetReducedVec(ff,FreeVariableSet);CHKERRQ(info);
136:       info = r1->Negate();CHKERRQ(info);
137: 
138:       info = dxfree->SetReducedVec(d,FreeVariableSet);CHKERRQ(info);
139:       info = dxfree->SetToZero(); CHKERRQ(info);
140: 
141:       info = TaoComputeJacobian(tao,x,J);CHKERRQ(info);
142:       info = Jsub->SetReducedMatrix(J,FreeVariableSet,FreeVariableSet);CHKERRQ(info);
143: 
144:       success = TAO_FALSE;
145:       gamma = gamma_factor*(ndpsi);
146:       while (success==TAO_FALSE && attempts < 10) {
147: 
148:         /* Approximately solve the reduced linear system */
149:         info = TaoPreLinearSolve(tao,Jsub);CHKERRQ(info);
150:         info = TaoLinearSolve(tao,Jsub,r1,dxfree,&success);CHKERRQ(info);
151: 
152:         info = d->SetToZero(); CHKERRQ(info);
153:         info = d->ReducedXPY(dxfree,FreeVariableSet);CHKERRQ(info);
154: 
155:         info = d->Dot(ff,&gdx); CHKERRQ(info);
156: 
157:         if (success==TAO_FALSE) { /* Modify diagonal of Hessian if not a descent direction */
158: 
159:           info = Jsub->SetReducedMatrix(J,FreeVariableSet,FreeVariableSet);CHKERRQ(info);
160:           gamma *= 10;
161:           //          printf("Shift diagonal: %4.2e\n",gamma);
162:           info = PetscLogInfo((tao,"TaoSolve_NLS:  modify diagonal (asuume same nonzero structure), gamma_factor=%g, gamma=%g\n",
163:                                gamma_factor,gamma)); CHKERRQ(info);
164:           info = Jsub->ShiftDiagonal(gamma);CHKERRQ(info);
165:           success = TAO_FALSE;
166: 
167:         } else {
168:           success = TAO_TRUE;
169:         }
170:         ++attempts;
171:       }
172: 
173:       //      fref=fff;
174:       info = g->ScaleCopyFrom(-1.0,d);CHKERRQ(info);
175:       stepsize=1.0;
176:       info = TaoLineSearchApply(tao,x,g,d,w,
177:                                 &fff,&stepsize,&lsflag);
178:       CHKERRQ(info);
179: 
180:     }

182:     if ( iter>=40 && (iter%10==0) ){
183:       printf("Steepest Descent \n");
184:       info = d->ScaleCopyFrom(-1.0,ff);CHKERRQ(info);
185:       info = g->CopyFrom(ff);CHKERRQ(info);
186:       stepsize=1.0;
187:       info = TaoLineSearchApply(tao,x,g,d,w,
188:                                 &fff,&stepsize,&lsflag);
189:       CHKERRQ(info);
190:       stepsize = 1.0e-6;
191:     }

193:     if (lsflag!=0 ){
194:       printf("Steepest Descent \n");
195:       info = d->ScaleCopyFrom(-1.0,ff);CHKERRQ(info);
196:       info = g->CopyFrom(ff);CHKERRQ(info);
197:       stepsize=1.0;
198:       info = TaoLineSearchApply(tao,x,g,d,w,
199:                                 &fff,&stepsize,&lsflag);
200:       CHKERRQ(info);
201:       stepsize = 1.0e-6;
202: 
203:     }

205:     printf("Stepsize: %4.2e\n",stepsize);
206:     /*
207:     if (fff>= fref){
208:       if (rref<2){ tao->lsflag=0;} 
209:       rref++;
210:     } else {
211:       rref=0;
212:     }
213:     */
214: 
215:   }  /* END MAIN LOOP  */


218:   info = TaoMatDestroy(Jsub);CHKERRQ(info);
219:   info = TaoVecDestroy(dxfree);CHKERRQ(info);
220:   info = TaoVecDestroy(r1);CHKERRQ(info);
221:   info = TaoIndexSetDestroy(FreeVariableSet);CHKERRQ(info);

223:   TaoFunctionReturn(0);
224: }


227: /* ---------------------------------------------------------- */
231: int TaoCreate_RSCS(TAO_SOLVER tao)
232: {
233:   TAO_SSLS *asls;
234:   int        info;

236:   TaoFunctionBegin;

238:   info = TaoNew(TAO_SSLS, &asls); CHKERRQ(info);
239:   info = PetscLogObjectMemory(tao, sizeof(TAO_SSLS)); CHKERRQ(info);

241:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_RSCS,(void*)asls); CHKERRQ(info);
242:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_SSLS,TaoSetDown_SSLS); CHKERRQ(info);
243:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_SSLS); CHKERRQ(info);
244:   info=TaoSetTaoViewRoutine(tao,TaoView_SSLS); CHKERRQ(info);

246:   // info = TaoCreateMoreThuenteBoundLineSearch(tao,0,0.9); CHKERRQ(info);
247:   info = TaoCreateNDProjectedArmijoLineSearch(tao); CHKERRQ(info);
248:   //   info = TaoCreateProjectedArmijoLineSearch(tao); CHKERRQ(info);
249:   info = TaoSetMeritFunction(tao, Tao_RSCS_Function, Tao_RSCS_FunctionGradient,
250:                              TAO_NULL, TAO_NULL, asls); CHKERRQ(info);

252:   info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
253:   info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);

255:   info = TaoSetTolerances(tao,0,0,0,0); CHKERRQ(info);
256:   info = TaoSetGradientTolerances(tao,1.0e-16,0.0,0.0); CHKERRQ(info);
257:   info = TaoSetFunctionLowerBound(tao,1.0e-8); CHKERRQ(info);

259:   TaoFunctionReturn(0);
260: }

263: /*------------------------------------------------------------*/
266: int Tao_RSCS_Function(TAO_SOLVER tao, TaoVec *X, double *fcn, void *solver)
267: {
268:   int info;
269:   double ndpsi;
270:   TaoVec *TT,*XL,*XU;
271:   FGMeritCtx*ctx = (FGMeritCtx*)solver;

273:   TaoFunctionBegin;
274:   info = TaoGetVariableBounds(tao, &XL, &XU); CHKERRQ(info);
275:   info = X->Clone(&TT); CHKERRQ(info);
276:   info = TaoComputeConstraints(tao, X, ctx->f); CHKERRQ(info);
277:   info = TT->Fischer(X, ctx->f, XL, XU); CHKERRQ(info);
278:   //  info = TT->BoundGradientProjection(ctx->f,XL,X,XU);CHKERRQ(info);
279:   info = TT->Norm2(&ndpsi); CHKERRQ(info);
280:   *fcn=ndpsi;
281:   info = PetscLogInfo((tao,"RSCS Function value: %4.2e\n",*fcn)); CHKERRQ(info);
282:   info = TaoVecDestroy(TT); CHKERRQ(info);
283:   TaoFunctionReturn(0);
284: }

286: /*------------------------------------------------------------*/
289: int Tao_RSCS_FunctionGradient(TAO_SOLVER tao, TaoVec *X, double *fcn, 
290:                               TaoVec *G, void *solver)
291: {
292:   int info;
293:   double ndpsi;
294:   TaoVec *XL,*XU,*TT,*f;
295:   TaoMat *J;
296:   FGMeritCtx*ctx = (FGMeritCtx*)solver;

298:   TaoFunctionBegin;
299:   J=ctx->J;
300:   f=ctx->f;
301:   info = TaoGetVariableBounds(tao, &XL, &XU); CHKERRQ(info);
302:   info = TaoComputeConstraints(tao, X, f); CHKERRQ(info);
303:   //  info = TaoComputeJacobian(tao,X,J); CHKERRQ(info);

305:   if (0==1){
306:     info = G->BoundGradientProjection(f,XL,X,XU);CHKERRQ(info);
307:     info = G->Norm2(&ndpsi); CHKERRQ(info);
308:     info = G->CopyFrom(f);CHKERRQ(info);
309:     *fcn=ndpsi;
310:   } else if (2==2){
311:     info = G->Fischer(X, f, XL, XU); CHKERRQ(info);
312:     info = G->Norm2(fcn); CHKERRQ(info);
313:     info = G->CopyFrom(f);CHKERRQ(info);
314:   } else {
315:     info = G->Clone(&TT); CHKERRQ(info);
316:     info = TT->BoundGradientProjection(f,XL,X,XU);CHKERRQ(info);
317:     info = TT->Norm2(&ndpsi); CHKERRQ(info);
318:     *fcn=ndpsi*ndpsi/2;
319:     *fcn=ndpsi;
320:     info = J->MultiplyTranspose(TT,G);CHKERRQ(info);
321:     info = TaoVecDestroy(TT); CHKERRQ(info);
322:   }

324:   info = PetscLogInfo((tao,"RSCS Function value: %4.2e\n",*fcn)); CHKERRQ(info);

326:   TaoFunctionReturn(0);
327: }