Actual source code: cg_pr.c

  1: /*$Id: cg_pr.c 1.21 05/05/10 15:21:38-05:00 sarich@zorak.(none) $*/

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

  7: static int TaoSolve_CG_PR(TAO_SOLVER tao, void *solver)
  8: {
  9:   TAO_CG *cg = (TAO_CG *) solver;
 10:   TaoVec*    X,*G; /* solution vector,  gradient vector */
 11:   TaoVec*    Gprev=cg->Gprev;
 12:   TaoVec*    DX=cg->DX, *Work=cg->Work;
 13:   int    iter=0,lsflag=0,info;
 14:   double gnorm2Prev,gdx,gDotgprev,f,gnorm,gnorm2,step=0;
 15:   TaoTerminateReason reason;

 17:   TaoFunctionBegin;
 18:   info=TaoGetSolution(tao,&X);CHKERRQ(info);
 19:   info=TaoGetGradient(tao,&G);CHKERRQ(info);

 21:   info = TaoComputeMeritFunctionGradient(tao,X,&f,G);CHKERRQ(info);
 22:   info = G->Norm2(&gnorm);  CHKERRQ(info);

 24:   info = DX->SetToZero(); CHKERRQ(info);
 25:   info = Gprev->CopyFrom(G); CHKERRQ(info);

 27:   cg->beta=0; cg->restarts=0;
 28:   gnorm2=gnorm*gnorm; gnorm2Prev = gnorm2;

 30:   /* Enter loop */
 31:   while (1){

 33:     /* Test for convergence */
 34:     info = TaoMonitor(tao,iter++,f,gnorm,0.0,step,&reason);CHKERRQ(info);
 35:     if (reason!=TAO_CONTINUE_ITERATING) break;

 37:     info = G->Dot(Gprev,&gDotgprev); CHKERRQ(info);
 38:     cg->beta=( gnorm2-gDotgprev )/(gnorm2Prev);
 39:     info = DX->Axpby(-1.0,G,cg->beta); CHKERRQ(info);
 40: 
 41:     info = DX->Dot(G,&gdx); CHKERRQ(info);
 42:     if (gdx>=0){ /* If not a descent direction, use gradient */
 43:       info = DX->ScaleCopyFrom(-1.0,G); CHKERRQ(info);
 44:       gdx=-gnorm2;
 45:       cg->restarts++; cg->beta=0.0;
 46:       info = PetscLogInfo((tao,"TaoCG: Not Descent Direction: Restart CG at iterate %d with gradient direction.\n",iter)); CHKERRQ(info);
 47:     }

 49:     if ( fabs(gDotgprev)/(gnorm2) > cg->eta){
 50:       /* The gradients are somewhat orthogonal */
 51:       gdx=-gnorm2;
 52:       cg->beta=0;
 53:       info = DX->Axpby(-1.0,G,0.0); CHKERRQ(info);
 54:       info = PetscLogInfo((tao,"TaoCG: Restart CG at iterate %d with gradient direction.\n",iter)); CHKERRQ(info);
 55:     }

 57:     /* Line Search */
 58:     gnorm2Prev = gnorm2;  step=1.0;
 59:     info = Gprev->CopyFrom(G); CHKERRQ(info);
 60:     info = TaoLineSearchApply(tao,X,G,DX,Work,&f,&step,&lsflag);CHKERRQ(info);
 61:     info = G->Norm2(&gnorm);CHKERRQ(info);
 62:     gnorm2=gnorm*gnorm;
 63:   }
 64: 
 65:   TaoFunctionReturn(0);
 66: }

 68: /* ---------------------------------------------------------- */
 69: /*
 70:     TAOCreate_CG_PR - Creates the data structure for the nonlinear CG method
 71:     and sets the function pointers for all the routines it needs to call
 72:     (TAOSolve_CG() etc.)

 75: */
 79: int TaoCreate_CG_PR(TAO_SOLVER tao)
 80: {
 81:   TAO_CG *cg;
 82:   int    info;

 84:   TaoFunctionBegin;

 86:   info = TaoNew(TAO_CG,&cg); CHKERRQ(info);
 87:   info = PetscLogObjectMemory(tao,sizeof(TAO_CG)); CHKERRQ(info);

 89:   info=TaoSetTaoSolveRoutine(tao,TaoSolve_CG_PR,(void*)cg); CHKERRQ(info);
 90:   info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_CG,TaoDestroy_CG); CHKERRQ(info);
 91:   info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_CG); CHKERRQ(info);
 92:   info=TaoSetTaoViewRoutine(tao,TaoView_CG); CHKERRQ(info);

 94:   info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
 95:   info = TaoSetTolerances(tao,1e-4,1e-4,0,0); CHKERRQ(info);
 96:   info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);

 98:   info = TaoCreateMoreThuenteLineSearch(tao,0,0.4); CHKERRQ(info);

100:   cg->eta = 0.5;

102:   TaoFunctionReturn(0);
103: }