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: }