Actual source code: gpcg.c

      1: /*$Id: gpcg.c 1.93 05/05/09 11:09:17-05:00 sarich@zorak.(none) $*/

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

      5: char gpcgname[]="tao_gpcg";
      6: TaoMethod gpcgtypename = gpcgname;

      8: static int TaoGradProjections(TAO_SOLVER, TAO_GPCG *);
      9: static int GPCGCheckOptimalFace(TaoVec*, TaoVec*, TaoVec*, TaoVec*, TaoVec*,
     10:                                 TaoIndexSet*, TaoIndexSet*, TaoTruth *optimal);

     12: /*------------------------------------------------------------*/
     15: static int TaoSetDown_GPCG(TAO_SOLVER tao, void*solver)
     16: {
     17:   TAO_GPCG *gpcg = (TAO_GPCG *)solver;
     18:   int      info;
     19:   /* Free allocated memory in GPCG structure */
     20:   TaoFunctionBegin;
     21:
     22:   info = TaoVecDestroy(gpcg->X_New);CHKERRQ(info);
     23:   info = TaoVecDestroy(gpcg->G_New);CHKERRQ(info);
     24:   info = TaoVecDestroy(gpcg->DX);CHKERRQ(info);gpcg->DX=0;
     25:   info = TaoVecDestroy(gpcg->Work);CHKERRQ(info);
     26:   info = TaoVecDestroy(gpcg->DXFree);CHKERRQ(info);
     27:   info = TaoVecDestroy(gpcg->R);CHKERRQ(info);
     28:   info = TaoVecDestroy(gpcg->B);CHKERRQ(info);
     29:   info = TaoVecDestroy(gpcg->G);CHKERRQ(info);
     30:   info = TaoVecDestroy(gpcg->PG);CHKERRQ(info);
     31:   info = TaoVecDestroy(gpcg->XL);CHKERRQ(info);
     32:   info = TaoVecDestroy(gpcg->XU);CHKERRQ(info);
     33:
     34:   info = TaoIndexSetDestroy(gpcg->Free_Local);CHKERRQ(info);
     35:   info = TaoIndexSetDestroy(gpcg->TT);CHKERRQ(info);
     36:   info = TaoMatDestroy(gpcg->Hsub);CHKERRQ(info);
     37:
     38:   info = TaoDestroyLinearSolver(tao);CHKERRQ(info);

     40:   TaoFunctionReturn(0);
     41: }

     43: /*------------------------------------------------------------*/
     46: static int TaoSetFromOptions_GPCG(TAO_SOLVER tao, void*solver)
     47: {
     48:   TAO_GPCG *gpcg = (TAO_GPCG *)solver;
     49:   int      info,ival;
     50:   TaoTruth flg;

     52:   TaoFunctionBegin;
     53:   info = TaoOptionsHead("Gradient Projection, Conjugate Gradient method for bound constrained optimization");CHKERRQ(info);

     55:   info=TaoOptionInt("-gpcg_maxpgits","maximum number of gradient projections per GPCG iterate",0,gpcg->maxgpits,&gpcg->maxgpits,&flg);
     56:   CHKERRQ(info);

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

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

     63:   info = TaoOptionsTail();CHKERRQ(info);
     64:   info=TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
     65:   /*  info = TaoSetLinearSolverOptions(tao);CHKERRQ(info); */

     67:   TaoFunctionReturn(0);
     68: }

     70: /*------------------------------------------------------------*/
     73: static int TaoView_GPCG(TAO_SOLVER tao,void*solver)
     74: {
     75:   TAO_GPCG *gpcg = (TAO_GPCG *)solver;
     76:   int      info;

     78:   TaoFunctionBegin;

     80:   info = TaoPrintInt(tao," Total PG its: %d,",gpcg->total_gp_its);CHKERRQ(info);
     81:   info = TaoPrintDouble(tao," PG tolerance: %4.3f \n",gpcg->pg_ftol);CHKERRQ(info);
     82:   info = TaoLineSearchView(tao);CHKERRQ(info);

     84:   TaoFunctionReturn(0);
     85: }

     87: /* ---------------------------------------------------------- */
     90: int TaoGPCGComputeFunctionGradient(TAO_SOLVER tao, TaoVec *XX, double *ff, TaoVec *GG){
     91:   TAO_GPCG *gpcg;
     92:   TaoMat *HH;
     93:   int info;
     94:   double f1,f2;

     96:   TaoFunctionBegin;
     97:   info = TaoGetSolverContext(tao,gpcgtypename,(void**)&gpcg); CHKERRQ(info);
     98:   if (gpcg==0){TaoFunctionReturn(0);}
     99:   info = TaoGetHessian(tao,&HH);CHKERRQ(info);
    100:   info = HH->Multiply(XX,GG);CHKERRQ(info);
    101:   info = XX->Dot(GG,&f1);CHKERRQ(info);
    102:   info = XX->Dot(gpcg->B,&f2);CHKERRQ(info);
    103:   info = GG->Axpy(1.0,gpcg->B);CHKERRQ(info);
    104:   *ff=f1/2.0 + f2 + gpcg->c;
    105:   TaoFunctionReturn(0);
    106: }

    108: /* ---------------------------------------------------------- */
    111: static int TaoSetUp_GPCG(TAO_SOLVER tao,void*solver){

    113:   int      n,info;
    114:   TAO_GPCG *gpcg = (TAO_GPCG *) solver;
    115:   TaoVec   *X;
    116:   TaoMat   *HH;
    117:   TaoIndexSet *TIS;
    118:   TaoLinearSolver *ksp;

    120:   TaoFunctionBegin;

    122:   info = TaoGetSolution(tao,&X);CHKERRQ(info); gpcg->X=X;
    123:   info = TaoGetHessian(tao,&HH);CHKERRQ(info); gpcg->H=HH;

    125:   /* Allocate some arrays */
    126:   info=X->Clone(&gpcg->DX); CHKERRQ(info);
    127:   info=X->Clone(&gpcg->B); CHKERRQ(info);
    128:   info=X->Clone(&gpcg->Work); CHKERRQ(info);
    129:   info=X->Clone(&gpcg->X_New); CHKERRQ(info);
    130:   info=X->Clone(&gpcg->G_New); CHKERRQ(info);
    131:   info=X->Clone(&gpcg->DXFree); CHKERRQ(info);
    132:   info=X->Clone(&gpcg->R); CHKERRQ(info);
    133:   info=X->Clone(&gpcg->G); CHKERRQ(info);
    134:   info=X->Clone(&gpcg->PG); CHKERRQ(info);
    135:   info=X->Clone(&gpcg->XL); CHKERRQ(info);
    136:   info=X->Clone(&gpcg->XU); CHKERRQ(info);

    138:   info = TaoSetLagrangianGradientVector(tao,gpcg->PG);CHKERRQ(info);
    139:   info = TaoSetStepDirectionVector(tao,gpcg->DX);CHKERRQ(info);
    140:   info = TaoSetVariableBounds(tao,gpcg->XL,gpcg->XU);CHKERRQ(info);

    142:   info = X->GetDimension(&n); CHKERRQ(info);
    143:   gpcg->n=n;
    144:   info = TaoCreateLinearSolver(tao,HH,300,0); CHKERRQ(info);
    145:   info = TaoGetLinearSolver(tao,&ksp); CHKERRQ(info);
    146:   info = ksp->SetTolerances(1e-2,1e-30,1e30,n);CHKERRQ(info);

    148:   info = X->CreateIndexSet(&TIS); CHKERRQ(info);
    149:   gpcg->Free_Local = TIS;
    150:   info = gpcg->Free_Local->Duplicate(&gpcg->TT); CHKERRQ(info);

    152:   info = HH->CreateReducedMatrix(TIS,TIS,&gpcg->Hsub); CHKERRQ(info);

    154:   TaoFunctionReturn(0);
    155: }

    159: static int TaoSolve_GPCG(TAO_SOLVER tao, void *solver)
    160: {
    161:   TAO_GPCG *gpcg = (TAO_GPCG *)solver ;
    162:   int info,lsflag,iter=0;
    163:   TaoTruth optimal_face=TAO_FALSE,success;
    164:   double actred,f,f_new,gnorm,gdx,stepsize;
    165:   double c;
    166:   TaoVec *XU, *XL;
    167:   TaoVec *X,  *G=gpcg->G , *B=gpcg->B, *PG=gpcg->PG;
    168:   TaoVec *R=gpcg->R, *DXFree=gpcg->DXFree;
    169:   TaoVec *X_New=gpcg->X_New, *G_New=gpcg->G_New;
    170:   TaoVec *DX=gpcg->DX, *Work=gpcg->Work;
    171:   TaoMat *H, *Hsub=gpcg->Hsub;
    172:   TaoIndexSet *Free_Local = gpcg->Free_Local, *TIS=gpcg->TT;
    173:   TaoTerminateReason reason;

    175:   TaoFunctionBegin;

    177:   /* Check if upper bound greater than lower bound. */
    178:   info = TaoGetSolution(tao,&X);CHKERRQ(info);
    179:   info = TaoGetHessian(tao,&H);CHKERRQ(info);

    181:   info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
    182:   info = TaoEvaluateVariableBounds(tao,XL,XU); CHKERRQ(info);
    183:   info = X->Median(XL,X,XU); CHKERRQ(info);

    185:   info = TaoComputeHessian(tao,X,H); CHKERRQ(info);
    186:   info = TaoComputeFunctionGradient(tao,X,&f,B);
    187:   CHKERRQ(info);

    189:   /* Compute quadratic representation */
    190:   info = H->Multiply(X,Work); CHKERRQ(info);
    191:   info = X->Dot(Work,&c); CHKERRQ(info);
    192:   info = B->Axpy(-1.0,Work); CHKERRQ(info);
    193:   info = X->Dot(B,&stepsize); CHKERRQ(info);
    194:   gpcg->c=f-c/2.0-stepsize;

    196:   info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
    197:
    198:   info = TaoGPCGComputeFunctionGradient(tao, X, &gpcg->f , G);
    199:
    200:   /* Project the gradient and calculate the norm */
    201:   info = G_New->CopyFrom(G);CHKERRQ(info);
    202:   info = PG->BoundGradientProjection(G,XL,X,XU);CHKERRQ(info);
    203:   info = PG->Norm2(&gpcg->gnorm); CHKERRQ(info);
    204:   gpcg->step=1.0;

    206:     /* Check Stopping Condition      */
    207:   info=TaoMonitor(tao,iter++,gpcg->f,gpcg->gnorm,0,gpcg->step,&reason); CHKERRQ(info);

    209:   while (reason == TAO_CONTINUE_ITERATING){

    211:     info = TaoGradProjections(tao, gpcg); CHKERRQ(info);

    213:     info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
    214:     info = Free_Local->GetSize(&gpcg->n_free); CHKERRQ(info);
    215:     f=gpcg->f; gnorm=gpcg->gnorm;

    217:     if (gpcg->n_free > 0){
    218:
    219:       /* Create a reduced linear system */
    220:       info = R->SetReducedVec(G,Free_Local);CHKERRQ(info);
    221:       info = R->Negate(); CHKERRQ(info);
    222:       info = DXFree->SetReducedVec(DX,Free_Local);CHKERRQ(info);
    223:       info = DXFree->SetToZero(); CHKERRQ(info);

    225:       info = Hsub->SetReducedMatrix(H,Free_Local,Free_Local);CHKERRQ(info);

    227:       info = TaoPreLinearSolve(tao,Hsub);CHKERRQ(info);

    229:       while (1){
    230:
    231:         /* Approximately solve the reduced linear system */
    232:         info = TaoLinearSolve(tao,Hsub,R,DXFree,&success);CHKERRQ(info);
    233:
    234:         info=DX->SetToZero(); CHKERRQ(info);
    235:         info=DX->ReducedXPY(DXFree,Free_Local);CHKERRQ(info);

    237:         info = G->Dot(DX,&gdx); CHKERRQ(info);

    239:         stepsize=1.0; f_new=f;
    240:         info = X_New->CopyFrom(X); CHKERRQ(info);
    241:         info = G_New->CopyFrom(G); CHKERRQ(info);
    242:
    243:         info = TaoLineSearchApply(tao,X_New,G_New,DX,Work,
    244:                                   &f_new,&stepsize,&lsflag);
    245:         CHKERRQ(info);

    247:         actred = f_new - f;
    248:
    249:         /* Evaluate the function and gradient at the new point */
    250:         info =  PG->BoundGradientProjection(G_New,XL,X_New,XU);
    251:         CHKERRQ(info);
    252:         info = PG->Norm2(&gnorm);  CHKERRQ(info);

    254:         info = GPCGCheckOptimalFace(X_New,XL,XU,PG,Work, Free_Local, TIS,
    255:                                     &optimal_face); CHKERRQ(info);

    257:         info = TaoMonitor(tao,iter,f_new,gnorm,0.0,stepsize,&reason); CHKERRQ(info);

    259:         if (stepsize < 1 || optimal_face==TAO_FALSE || reason!=TAO_CONTINUE_ITERATING){
    260:           f=f_new;
    261:           info = X->CopyFrom(X_New); CHKERRQ(info);
    262:           info = G->CopyFrom(G_New); CHKERRQ(info);
    263:           break;
    264:         }

    266:       } /* end linear solve loop */
    267:
    268:
    269:     } else {
    270:
    271:       actred = 0;
    272:       /* if there were no free variables, no cg method */
    273:       info = TaoMonitor(tao,iter,f,gnorm,0.0,1.0,&reason); CHKERRQ(info);

    275:     }

    277:     gpcg->f=f;gpcg->gnorm=gnorm; gpcg->actred=actred;
    278:     if (reason!=TAO_CONTINUE_ITERATING) break;
    279:     iter++;

    282:   }  /* END MAIN LOOP  */

    284:   TaoFunctionReturn(0);
    285: }

    289: static int TaoGradProjections(TAO_SOLVER tao, TAO_GPCG *gpcg)
    290: {
    291:   int i,lsflag=0,info;
    292:   TaoTruth optimal_face=TAO_FALSE;
    293:   double actred=-1.0,actred_max=0.0, gAg,gtg=gpcg->gnorm,alpha;
    294:   double f_new,gdx;
    295:   TaoMat *H;
    296:   TaoVec *DX=gpcg->DX,*XL=gpcg->XL,*XU=gpcg->XU,*Work=gpcg->Work;
    297:   TaoVec *X=gpcg->X,*G=gpcg->G;
    298:   /*
    299:      The gradient and function value passed into and out of this
    300:      routine should be current and correct.
    301:
    302:      The free, active, and binding variables should be already identified
    303:   */
    304:
    305:   TaoFunctionBegin;
    306:
    307:   info = TaoGetSolution(tao,&X);CHKERRQ(info);
    308:   info = TaoGetHessian(tao,&H);CHKERRQ(info);
    309:   info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);

    311:   for (i=0;i<gpcg->maxgpits;i++){

    313:     if ( -actred <= (gpcg->pg_ftol)*actred_max) break;
    314:
    315:     info = DX->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);
    316:     info = DX->Negate(); CHKERRQ(info);
    317:     info = DX->Dot(G,&gdx); CHKERRQ(info);

    319:     info= H->Multiply(DX,Work); CHKERRQ(info);
    320:     info= DX->Dot(Work,&gAg); CHKERRQ(info);
    321:
    322:     gpcg->gp_iterates++; gpcg->total_gp_its++;
    323:
    324:     gtg=-gdx;
    325:     alpha = TaoAbsDouble(gtg/gAg);
    326:     gpcg->stepsize = alpha; f_new=gpcg->f;

    328:     info = TaoLineSearchApply(tao,X,G,DX,Work,
    329:                               &f_new,&gpcg->stepsize,&lsflag);
    330:     CHKERRQ(info);

    332:     /* Update the iterate */
    333:     actred = f_new - gpcg->f;
    334:     actred_max = TaoMax(actred_max,-(f_new - gpcg->f));
    335:     gpcg->f = f_new;
    336:     info = GPCGCheckOptimalFace(X,XL,XU,G,Work,gpcg->Free_Local,gpcg->TT,
    337:                                 &optimal_face); CHKERRQ(info);

    339:     if ( optimal_face == TAO_TRUE ) break;

    341:   }
    342:
    343:   gpcg->gnorm=gtg;
    344:   TaoFunctionReturn(0);

    346: } /* End gradient projections */

    351: static int GPCGCheckOptimalFace(TaoVec *X,TaoVec *XL,TaoVec*XU,TaoVec *PG,TaoVec*W,
    352:                                 TaoIndexSet*Free_Local, TaoIndexSet*TT,
    353:                                 TaoTruth *optimal)
    354: {
    355:   int info,n_free;
    356:   double rr;
    357:   TaoTruth same;

    359:   TaoFunctionBegin;
    360:   *optimal = TAO_FALSE;

    362:   /* Also check to see if the active set is the same */

    364:   info = TT->WhichBetween(XL,X,XU); CHKERRQ(info);
    365:   info = Free_Local->IsSame(TT,&same); CHKERRQ(info);
    366:   info = Free_Local->GetSize(&n_free); CHKERRQ(info);
    367:   if (same == TAO_FALSE){
    368:     info = Free_Local->WhichBetween(XL,X,XU); CHKERRQ(info);
    369:     *optimal = TAO_FALSE;
    370:     TaoFunctionReturn(0);
    371:   } else {
    372:     *optimal = TAO_TRUE;
    373:   }

    376:   info = W->CopyFrom(PG); CHKERRQ(info);
    377:   info = W->Negate(); CHKERRQ(info);

    379:   info = W->BoundGradientProjection(W,XL,X,XU); CHKERRQ(info);
    380:   info = W->Axpy(1.0,PG); CHKERRQ(info);

    382:   info = W->Norm2(&rr); CHKERRQ(info);
    383:   if (rr>0) *optimal = TAO_FALSE;

    385:   *optimal = TAO_FALSE;
    386:   /*
    387:   info = gpcg->TT->whichNonNegative(W); CHKERRQ(info);
    388:   info = gpcg->TT->GetSize(&n); CHKERRQ(info);
    389:   if (n==0) *optimal = TAO_TRUE;
    390:   */
    391:   TaoFunctionReturn(0);
    392: }

    394: int TaoDefaultMonitor_GPCG(TAO_SOLVER tao,void *dummy)
    395: {
    396:   TAO_GPCG *gpcg;
    397:   double   fct,gnorm;
    398:   int      its,nfree,info;

    400:   TaoFunctionBegin;
    401:   info = TaoGetSolutionStatus(tao,&its,&fct,&gnorm,0,0,0);CHKERRQ(info);
    402:   info = TaoGetSolverContext(tao,gpcgtypename,(void**)&gpcg); CHKERRQ(info);
    403:   if (gpcg==0){TaoFunctionReturn(0);}
    404:   nfree=gpcg->n_free;
    405:   info = TaoPrintInt(tao,"iter: %d,",its);CHKERRQ(info);
    406:   info = TaoPrintDouble(tao," Fcn value: %g,",fct);CHKERRQ(info);
    407:   info = TaoPrintDouble(tao," PGrad. norm: %g, ",gnorm);CHKERRQ(info);
    408:   info = TaoPrintInt(tao,"free vars:%d \n",nfree);CHKERRQ(info);
    409:   TaoFunctionReturn(0);
    410: }

    414: int TaoGetDualVariables_GPCG(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void* solver)
    415: {

    417:   TAO_GPCG *gpcg = (TAO_GPCG *) solver;
    418:   TaoVec  *G=gpcg->G,*GP=gpcg->Work;
    419:   TaoVec  *X,*XL,*XU;
    420:   int       info;

    422:   TaoFunctionBegin;
    423:   info = TaoGetVariableBounds(tao,&XL,&XU); CHKERRQ(info);
    424:   info = TaoGetSolution(tao,&X); CHKERRQ(info);
    425:   info = GP->BoundGradientProjection(G,XL,X,XU); CHKERRQ(info);

    427:   info = DXL->Waxpby(-1.0,G,1.0,GP); CHKERRQ(info);
    428:   info = DXU->SetToZero(); CHKERRQ(info);
    429:   info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);

    431:   info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
    432:   info = GP->SetToZero(); CHKERRQ(info);
    433:   info = DXU->PointwiseMinimum(GP,DXU); CHKERRQ(info);

    435:   TaoFunctionReturn(0);
    436: }

    438: /*------------------------------------------------------------*/
    442: int TaoCreate_GPCG(TAO_SOLVER tao)
    443: {
    444:   TAO_GPCG *gpcg;
    445:   int      info;

    447:   TaoFunctionBegin;

    449:   info = TaoNew(TAO_GPCG,&gpcg); CHKERRQ(info);
    450:   info = PetscLogObjectMemory(tao,sizeof(TAO_GPCG)); CHKERRQ(info);

    452:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_GPCG,(void*)gpcg); CHKERRQ(info);
    453:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_GPCG,TaoSetDown_GPCG); CHKERRQ(info);
    454:   info=TaoSetTaoOptionsRoutine(tao,TaoSetFromOptions_GPCG); CHKERRQ(info);
    455:   info=TaoSetTaoViewRoutine(tao,TaoView_GPCG); CHKERRQ(info);
    456:   info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_GPCG); CHKERRQ(info);

    458:   info = TaoSetMaximumIterates(tao,500); CHKERRQ(info);
    459:   info = TaoSetMaximumFunctionEvaluations(tao,100000); CHKERRQ(info);
    460:   info = TaoSetTolerances(tao,1e-12,1e-12,0,0); CHKERRQ(info);

    462:   /* Initialize pointers and variables */
    463:   gpcg->n=0;
    464:   gpcg->maxgpits = 8;
    465:   gpcg->pg_ftol = 0.1;

    467:   gpcg->gp_iterates=0; /* Cumulative number */
    468:   gpcg->total_gp_its = 0;
    469:
    470:   /* Initialize pointers and variables */
    471:   gpcg->n_bind=0;
    472:   gpcg->n_free = 0;
    473:   gpcg->n_upper=0;
    474:   gpcg->n_lower=0;

    476:   //  info = TaoCreateProjectedLineSearch(tao); CHKERRQ(info);
    477:   info = TaoGPCGCreateLineSearch(tao); CHKERRQ(info);

    479:   TaoFunctionReturn(0);
    480: }