Actual source code: cg_prp.c
1: /*$Id: cg_prp.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_PRP(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);
40: if (cg->beta<0){
42: info = PetscLogInfo((tao,"TaoCG: Negative beta. Restart CG at iterate %d with gradient direction.\n",iter)); CHKERRQ(info);
43: cg->restarts++; cg->beta=0.0;
44: info = DX->Axpby(-1.0,G,cg->beta); CHKERRQ(info);
45: gdx=-gnorm2;
47: } else {
49: info = DX->Axpby(-1.0,G,cg->beta); CHKERRQ(info);
50:
51: info = DX->Dot(G,&gdx); CHKERRQ(info);
52: if (gdx>=0){ /* If not a descent direction, use gradient */
53: info = DX->ScaleCopyFrom(-1.0,G); CHKERRQ(info);
54: gdx=-gnorm2;
55: cg->restarts++; cg->beta=0.0;
56: info = PetscLogInfo((tao,"TaoCG: Not Descent Direction: Restart CG at iterate %d with gradient direction.\n",iter)); CHKERRQ(info);
57: }
58:
59: if ( fabs(gDotgprev)/(gnorm2) > cg->eta){
60: /* The gradients are somewhat orthogonal */
61: cg->beta=0;
62: gdx=-gnorm2;
63: info = DX->Axpby(-1.0,G,0.0); CHKERRQ(info);
64: info = PetscLogInfo((tao,"TaoCG: Restart CG at iterate %d with gradient direction.\n",iter)); CHKERRQ(info);
65: }
66: }
68: /* Line Search */
69: gnorm2Prev = gnorm2; step=1.0;
70: info = Gprev->CopyFrom(G); CHKERRQ(info);
71: info = TaoLineSearchApply(tao,X,G,DX,Work,&f,&step,&lsflag);CHKERRQ(info);
72: info = G->Norm2(&gnorm);CHKERRQ(info);
73: gnorm2 = gnorm*gnorm;
74: }
75:
76: TaoFunctionReturn(0);
77: }
79: /* ---------------------------------------------------------- */
80: /*
81: TAOCreate_CG_PRP - Creates the data structure for the nonlinear CG method
82: and sets the function pointers for all the routines it needs to call
83: (TAOSolve_CG_PRP() etc.)
86: */
90: int TaoCreate_CG_PRP(TAO_SOLVER tao)
91: {
92: TAO_CG *cg;
93: int info;
95: TaoFunctionBegin;
97: info = TaoNew(TAO_CG,&cg); CHKERRQ(info);
98: info = PetscLogObjectMemory(tao,sizeof(TAO_CG)); CHKERRQ(info);
100: info=TaoSetTaoSolveRoutine(tao,TaoSolve_CG_PRP,(void*)cg); CHKERRQ(info);
101: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_CG,TaoDestroy_CG); CHKERRQ(info);
102: info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_CG); CHKERRQ(info);
103: info=TaoSetTaoViewRoutine(tao,TaoView_CG); CHKERRQ(info);
105: info = TaoSetMaximumIterates(tao,2000); CHKERRQ(info);
106: info = TaoSetTolerances(tao,1e-4,1e-4,0,0); CHKERRQ(info);
107: info = TaoSetMaximumFunctionEvaluations(tao,4000); CHKERRQ(info);
109: info = TaoCreateMoreThuenteLineSearch(tao,0,0.4); CHKERRQ(info);
111: cg->eta = 0.5;
113: TaoFunctionReturn(0);
114: }