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