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,lsflag,iter=0;
     13:   TaoTerminateReason reason=TAO_CONTINUE_ITERATING;
     14:   double f,gnorm,gdx,stepsize=1.0;
     15:   TaoTruth success;
     16:   TaoVec *XU, *XL;
     17:   TaoVec *X,  *G=bnls->G, *PG=bnls->PG;
     18:   TaoVec *R=bnls->R, *DXFree=bnls->DXFree;
     19:   TaoVec *DX=bnls->DX, *Work=bnls->Work;
     20:   TaoMat *H, *Hsub=bnls->Hsub;
     21:   TaoIndexSet *FreeVariables = bnls->FreeVariables;

     23:   TaoFunctionBegin;

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

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

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

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

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

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

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

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

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

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

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

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

    100:   TaoFunctionReturn(0);
    101: }

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

    127:   TaoFunctionReturn(0);
    128: }

    130: /*------------------------------------------------------------*/
    133: static int TaoSetOptions_BNLS(TAO_SOLVER tao, void*solver)
    134: {
    135:   int        info,ival;
    136:   TaoTruth flg;

    138:   TaoFunctionBegin;

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

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

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

    146:   info = TaoOptionsTail();CHKERRQ(info);
    147:   info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);

    149:   TaoFunctionReturn(0);
    150: }

    152: /*------------------------------------------------------------*/
    155: static int TaoView_BNLS(TAO_SOLVER tao,void*solver)
    156: {
    157:   int        info;

    159:   TaoFunctionBegin;
    160:   info = TaoLineSearchView(tao);CHKERRQ(info);
    161:   TaoFunctionReturn(0);
    162: }

    165: /* ---------------------------------------------------------- */
    168: static int TaoSetUp_BNLS(TAO_SOLVER tao, void*solver){

    170:   int n,info;
    171:   TAO_BNLS *bnls = (TAO_BNLS *)solver;
    172:   TaoVec* X;
    173:   TaoMat *HH;
    174:   TaoLinearSolver *ksp;

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

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

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

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

    197:   info = TaoCreateLinearSolver(tao,HH,100,0); CHKERRQ(info);
    198:   info = TaoGetLinearSolver(tao,&ksp); CHKERRQ(info);
    199:   info = X->GetDimension(&n); CHKERRQ(info);
    200:   info = ksp->SetTolerances(1e-2,1e-30,1e30,n);CHKERRQ(info);

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

    204:   TaoFunctionReturn(0);
    205: }

    207: /*------------------------------------------------------------*/
    210: int TaoGetDualVariables_BNLS(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void* solver)
    211: {

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

    217:   TaoFunctionBegin;

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

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

    226:   TaoFunctionReturn(0);
    227: }

    230: /*------------------------------------------------------------*/
    234: int TaoCreate_BNLS(TAO_SOLVER tao)
    235: {
    236:   TAO_BNLS *bnls;
    237:   int      info;

    239:   TaoFunctionBegin;

    241:   info = TaoNew(TAO_BNLS,&bnls); CHKERRQ(info);
    242:   info = PetscLogObjectMemory(tao,sizeof(TAO_BNLS)); CHKERRQ(info);

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

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

    253:   /* Initialize pointers and variables */

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

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

    266:   TaoFunctionReturn(0);
    267: }