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;
13: TaoInt lsflag,iter=0;
14: TaoTerminateReason reason=TAO_CONTINUE_ITERATING;
15: double f,f_full,gnorm,gdx,stepsize=1.0;
16: TaoTruth success;
17: TaoVec *XU, *XL;
18: TaoVec *X, *G=bnls->G, *PG=bnls->PG;
19: TaoVec *R=bnls->R, *DXFree=bnls->DXFree;
20: TaoVec *DX=bnls->DX, *Work=bnls->Work;
21: TaoMat *H, *Hsub=bnls->Hsub;
22: TaoIndexSet *FreeVariables = bnls->FreeVariables;
24: TaoFunctionBegin;
26: /* Check if upper bound greater than lower bound. */
27: info = TaoGetSolution(tao,&X);CHKERRQ(info); bnls->X=X;
28: info = TaoGetVariableBounds(tao,&XL,&XU);CHKERRQ(info);
29: info = TaoEvaluateVariableBounds(tao,XL,XU); CHKERRQ(info);
30: info = TaoGetHessian(tao,&H);CHKERRQ(info); bnls->H=H;
32: /* Project the current point onto the feasible set */
33: info = X->Median(XL,X,XU); CHKERRQ(info);
34:
35: info = TaoComputeMeritFunctionGradient(tao,X,&f,G);CHKERRQ(info);
36:
37: while (reason==TAO_CONTINUE_ITERATING){
38:
39: /* Project the gradient and calculate the norm */
40: info = PG->BoundGradientProjection(G,XL,X,XU);CHKERRQ(info);
41: info = PG->Norm2(&gnorm); CHKERRQ(info);
42:
43: info = TaoMonitor(tao,iter++,f,gnorm,0.0,stepsize,&reason);
44: CHKERRQ(info);
45: if (reason!=TAO_CONTINUE_ITERATING) break;
47: info = FreeVariables->WhichEqual(PG,G); CHKERRQ(info);
49: info = TaoComputeHessian(tao,X,H);CHKERRQ(info);
50:
51: /* Create a reduced linear system */
53: info = R->SetReducedVec(G,FreeVariables);CHKERRQ(info);
54: info = R->Negate();CHKERRQ(info);
56: info = DXFree->SetReducedVec(DX,FreeVariables);CHKERRQ(info);
57: info = DXFree->SetToZero(); CHKERRQ(info);
58:
59: info = Hsub->SetReducedMatrix(H,FreeVariables,FreeVariables);CHKERRQ(info);
61: bnls->gamma_factor /= 2;
62: success = TAO_FALSE;
64: while (success==TAO_FALSE) {
65:
66: /* Approximately solve the reduced linear system */
67: info = TaoPreLinearSolve(tao,Hsub);CHKERRQ(info);
68: info = TaoLinearSolve(tao,Hsub,R,DXFree,&success);CHKERRQ(info);
70: info = DX->SetToZero(); CHKERRQ(info);
71: info = DX->ReducedXPY(DXFree,FreeVariables);CHKERRQ(info);
72: info = DX->Dot(G,&gdx); CHKERRQ(info);
74: if (gdx>=0 || success==TAO_FALSE) { /* Modify diagonal of Hessian if not a descent direction */
75: bnls->gamma_factor *= 2;
76: bnls->gamma = bnls->gamma_factor*(gnorm);
77: #if !defined(PETSC_USE_COMPLEX)
78: info=PetscInfo2(tao,"TaoSolve_NLS: modify diagonal (assume same nonzero structure), gamma_factor=%g, gamma=%g\n",bnls->gamma_factor,bnls->gamma);
79: CHKERRQ(info);
80: #else
81: info=PetscInfo3(tao,"TaoSolve_NLS: modify diagonal (asuume same nonzero structure), gamma_factor=%g, gamma=%g\n",
82: bnls->gamma_factor,PetscReal(bnls->gamma));CHKERRQ(info);
83: #endif
84: info = Hsub->ShiftDiagonal(bnls->gamma);CHKERRQ(info);
85: success = TAO_FALSE;
86:
87: } else {
88: success = TAO_TRUE;
89: }
91: }
92:
93: stepsize=1.0;
94: info = TaoLineSearchApply(tao,X,G,DX,Work,
95: &f,&f_full,&stepsize,&lsflag);
96: CHKERRQ(info);
98:
99: } /* END MAIN LOOP */
101: TaoFunctionReturn(0);
102: }
105: /*------------------------------------------------------------*/
108: static int TaoSetDown_BNLS(TAO_SOLVER tao, void*solver)
109: {
110: TAO_BNLS *bnls = (TAO_BNLS *)solver;
111: int info;
112: /* Free allocated memory in BNLS structure */
113: TaoFunctionBegin;
114:
115: info = TaoVecDestroy(bnls->DX);CHKERRQ(info);bnls->DX=0;
116: info = TaoVecDestroy(bnls->Work);CHKERRQ(info);
117: info = TaoVecDestroy(bnls->DXFree);CHKERRQ(info);
118: info = TaoVecDestroy(bnls->R);CHKERRQ(info);
119: info = TaoVecDestroy(bnls->G);CHKERRQ(info);
120: info = TaoVecDestroy(bnls->PG);CHKERRQ(info);
121: info = TaoVecDestroy(bnls->XL);CHKERRQ(info);
122: info = TaoVecDestroy(bnls->XU);CHKERRQ(info);
123:
124: info = TaoIndexSetDestroy(bnls->FreeVariables);CHKERRQ(info);
125: info = TaoMatDestroy(bnls->Hsub);CHKERRQ(info);
126: info = TaoDestroyLinearSolver(tao);CHKERRQ(info);
128: TaoFunctionReturn(0);
129: }
131: /*------------------------------------------------------------*/
134: static int TaoSetOptions_BNLS(TAO_SOLVER tao, void*solver)
135: {
136: int info;
137: TaoInt ival;
138: TaoTruth flg;
140: TaoFunctionBegin;
142: info = TaoOptionsHead("Newton Line Search Method for bound constrained optimization");CHKERRQ(info);
144: info = TaoOptionInt("-redistribute","Redistribute Free variables (> 1 processors, only)","TaoPetscISType",1,&ival,&flg); CHKERRQ(info);
146: info = TaoOptionName("-submatrixfree","Mask full matrix instead of extract submatrices","TaoPetscISType",&flg); CHKERRQ(info);
148: info = TaoOptionsTail();CHKERRQ(info);
149: info = TaoLineSearchSetFromOptions(tao);CHKERRQ(info);
151: TaoFunctionReturn(0);
152: }
154: /*------------------------------------------------------------*/
157: static int TaoView_BNLS(TAO_SOLVER tao,void*solver)
158: {
159: int info;
161: TaoFunctionBegin;
162: info = TaoLineSearchView(tao);CHKERRQ(info);
163: TaoFunctionReturn(0);
164: }
167: /* ---------------------------------------------------------- */
170: static int TaoSetUp_BNLS(TAO_SOLVER tao, void*solver){
172: int info;
173: TAO_BNLS *bnls = (TAO_BNLS *)solver;
174: TaoVec* X;
175: TaoMat *HH;
177: TaoFunctionBegin;
178: info = TaoGetSolution(tao,&bnls->X);CHKERRQ(info); X=bnls->X;
179: info = TaoGetHessian(tao,&bnls->H);CHKERRQ(info); HH=bnls->H;
181: /* Allocate some arrays */
182: info = X->Clone(&bnls->DX); CHKERRQ(info);
183: info = X->Clone(&bnls->Work); CHKERRQ(info);
184: info = X->Clone(&bnls->DXFree); CHKERRQ(info);
185: info = X->Clone(&bnls->R); CHKERRQ(info);
186: info = X->Clone(&bnls->G); CHKERRQ(info);
187: info = X->Clone(&bnls->PG); CHKERRQ(info);
188: info = X->Clone(&bnls->XL); CHKERRQ(info);
189: info = X->Clone(&bnls->XU); CHKERRQ(info);
191: info = TaoSetLagrangianGradientVector(tao,bnls->PG);CHKERRQ(info);
192: info = TaoSetStepDirectionVector(tao,bnls->DX);CHKERRQ(info);
193: info = TaoSetVariableBounds(tao,bnls->XL,bnls->XU);CHKERRQ(info);
195: info = X->CreateIndexSet(&bnls->FreeVariables); CHKERRQ(info);
196: info = bnls->H->CreateReducedMatrix(bnls->FreeVariables,bnls->FreeVariables,&bnls->Hsub); CHKERRQ(info);
198: info = TaoCreateLinearSolver(tao,HH,100,0); CHKERRQ(info);
200: info = TaoCheckFGH(tao);CHKERRQ(info);
202: TaoFunctionReturn(0);
203: }
205: /*------------------------------------------------------------*/
208: int TaoGetDualVariables_BNLS(TAO_SOLVER tao, TaoVec* DXL, TaoVec* DXU, void* solver)
209: {
211: TAO_BNLS *bnls = (TAO_BNLS *) solver;
212: TaoVec *G=bnls->G,*GP=bnls->Work;
213: int info;
215: TaoFunctionBegin;
217: info = DXL->Waxpby(-1,G,1.0,GP); CHKERRQ(info);
218: info = DXU->SetToZero(); CHKERRQ(info);
219: info = DXL->PointwiseMaximum(DXL,DXU); CHKERRQ(info);
221: info = DXU->Waxpby(-1.0,GP,1.0,G); CHKERRQ(info);
222: info = DXU->Axpy(1.0,DXL); CHKERRQ(info);
224: TaoFunctionReturn(0);
225: }
228: /*------------------------------------------------------------*/
232: int TaoCreate_BNLS(TAO_SOLVER tao)
233: {
234: TAO_BNLS *bnls;
235: int info;
237: TaoFunctionBegin;
239: info = TaoNew(TAO_BNLS,&bnls); CHKERRQ(info);
240: info = PetscLogObjectMemory(tao,sizeof(TAO_BNLS)); CHKERRQ(info);
242: info=TaoSetTaoSolveRoutine(tao,TaoSolve_BNLS,(void*)bnls); CHKERRQ(info);
243: info=TaoSetTaoSetUpDownRoutines(tao,TaoSetUp_BNLS,TaoSetDown_BNLS); CHKERRQ(info);
244: info=TaoSetTaoOptionsRoutine(tao,TaoSetOptions_BNLS); CHKERRQ(info);
245: info=TaoSetTaoViewRoutine(tao,TaoView_BNLS); CHKERRQ(info);
246: info=TaoSetTaoDualVariablesRoutine(tao,TaoGetDualVariables_BNLS); CHKERRQ(info);
248: info = TaoSetMaximumIterates(tao,500); CHKERRQ(info);
249: info = TaoSetTolerances(tao,1e-12,1e-12,0,0); CHKERRQ(info);
251: /* Initialize pointers and variables */
253: bnls->gamma = 0.0;
254: bnls->gamma_factor = 0.01;
255: bnls->DX=0;
256: bnls->DXFree=0;
257: bnls->R=0;
258: bnls->Work=0;
259: bnls->FreeVariables=0;
260: bnls->Hsub=0;
262: info = TaoCreateMoreThuenteBoundLineSearch(tao,0,0.9); CHKERRQ(info);
264: TaoFunctionReturn(0);
265: }