Actual source code: blmvm.c
petsc-3.5.4 2015-05-23
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/matrix/lmvmmat.h>
3: #include <../src/tao/bound/impls/blmvm/blmvm.h>
5: /*------------------------------------------------------------*/
8: static PetscErrorCode TaoSolve_BLMVM(Tao tao)
9: {
10: PetscErrorCode ierr;
11: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
12: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
13: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
14: PetscReal f, fold, gdx, gnorm;
15: PetscReal stepsize = 1.0,delta;
16: PetscInt iter = 0;
19: /* Project initial point onto bounds */
20: TaoComputeVariableBounds(tao);
21: VecMedian(tao->XL,tao->solution,tao->XU,tao->solution);
22: TaoLineSearchSetVariableBounds(tao->linesearch,tao->XL,tao->XU);
24: /* Check convergence criteria */
25: TaoComputeObjectiveAndGradient(tao, tao->solution,&f,blmP->unprojected_gradient);
26: VecBoundGradientProjection(blmP->unprojected_gradient,tao->solution, tao->XL,tao->XU,tao->gradient);
28: VecNorm(tao->gradient,NORM_2,&gnorm);
29: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN");
31: TaoMonitor(tao, iter, f, gnorm, 0.0, stepsize, &reason);
32: if (reason != TAO_CONTINUE_ITERATING) return(0);
34: /* Set initial scaling for the function */
35: if (f != 0.0) {
36: delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm);
37: } else {
38: delta = 2.0 / (gnorm*gnorm);
39: }
40: MatLMVMSetDelta(blmP->M,delta);
42: /* Set counter for gradient/reset steps */
43: blmP->grad = 0;
44: blmP->reset = 0;
46: /* Have not converged; continue with Newton method */
47: while (reason == TAO_CONTINUE_ITERATING) {
48: /* Compute direction */
49: MatLMVMUpdate(blmP->M, tao->solution, tao->gradient);
50: MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
51: VecBoundGradientProjection(tao->stepdirection,tao->solution,tao->XL,tao->XU,tao->gradient);
53: /* Check for success (descent direction) */
54: VecDot(blmP->unprojected_gradient, tao->gradient, &gdx);
55: if (gdx <= 0) {
56: /* Step is not descent or solve was not successful
57: Use steepest descent direction (scaled) */
58: ++blmP->grad;
60: if (f != 0.0) {
61: delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm);
62: } else {
63: delta = 2.0 / (gnorm*gnorm);
64: }
65: MatLMVMSetDelta(blmP->M,delta);
66: MatLMVMReset(blmP->M);
67: MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
68: MatLMVMSolve(blmP->M,blmP->unprojected_gradient, tao->stepdirection);
69: }
70: VecScale(tao->stepdirection,-1.0);
72: /* Perform the linesearch */
73: fold = f;
74: VecCopy(tao->solution, blmP->Xold);
75: VecCopy(blmP->unprojected_gradient, blmP->Gold);
76: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
77: TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);
78: TaoAddLineSearchCounts(tao);
80: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
81: /* Linesearch failed
82: Reset factors and use scaled (projected) gradient step */
83: ++blmP->reset;
85: f = fold;
86: VecCopy(blmP->Xold, tao->solution);
87: VecCopy(blmP->Gold, blmP->unprojected_gradient);
89: if (f != 0.0) {
90: delta = 2.0* PetscAbsScalar(f) / (gnorm*gnorm);
91: } else {
92: delta = 2.0/ (gnorm*gnorm);
93: }
94: MatLMVMSetDelta(blmP->M,delta);
95: MatLMVMReset(blmP->M);
96: MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
97: MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
98: VecScale(tao->stepdirection, -1.0);
100: /* This may be incorrect; linesearch has values fo stepmax and stepmin
101: that should be reset. */
102: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
103: TaoLineSearchApply(tao->linesearch,tao->solution,&f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);
104: TaoAddLineSearchCounts(tao);
106: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
107: tao->reason = TAO_DIVERGED_LS_FAILURE;
108: break;
109: }
110: }
112: /* Check for converged */
113: VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);
114: VecNorm(tao->gradient, NORM_2, &gnorm);
117: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Not-a-Number");
118: iter++;
119: TaoMonitor(tao, iter, f, gnorm, 0.0, stepsize, &reason);
120: }
121: return(0);
122: }
126: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
127: {
128: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
129: PetscInt n,N;
133: /* Existence of tao->solution checked in TaoSetup() */
134: VecDuplicate(tao->solution,&blmP->Xold);
135: VecDuplicate(tao->solution,&blmP->Gold);
136: VecDuplicate(tao->solution, &blmP->unprojected_gradient);
138: if (!tao->stepdirection) {
139: VecDuplicate(tao->solution, &tao->stepdirection);
140: }
141: if (!tao->gradient) {
142: VecDuplicate(tao->solution,&tao->gradient);
143: }
144: if (!tao->XL) {
145: VecDuplicate(tao->solution,&tao->XL);
146: VecSet(tao->XL,PETSC_NINFINITY);
147: }
148: if (!tao->XU) {
149: VecDuplicate(tao->solution,&tao->XU);
150: VecSet(tao->XU,PETSC_INFINITY);
151: }
152: /* Create matrix for the limited memory approximation */
153: VecGetLocalSize(tao->solution,&n);
154: VecGetSize(tao->solution,&N);
155: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&blmP->M);
156: MatLMVMAllocateVectors(blmP->M,tao->solution);
157: return(0);
158: }
160: /* ---------------------------------------------------------- */
163: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
164: {
165: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
169: if (tao->setupcalled) {
170: MatDestroy(&blmP->M);
171: VecDestroy(&blmP->unprojected_gradient);
172: VecDestroy(&blmP->Xold);
173: VecDestroy(&blmP->Gold);
174: }
175: PetscFree(tao->data);
176: return(0);
177: }
179: /*------------------------------------------------------------*/
182: static PetscErrorCode TaoSetFromOptions_BLMVM(Tao tao)
183: {
187: PetscOptionsHead("Limited-memory variable-metric method for bound constrained optimization");
188: TaoLineSearchSetFromOptions(tao->linesearch);
189: PetscOptionsTail();
190: return(0);
191: }
194: /*------------------------------------------------------------*/
197: static int TaoView_BLMVM(Tao tao, PetscViewer viewer)
198: {
199: TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
200: PetscBool isascii;
204: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
205: if (isascii) {
206: PetscViewerASCIIPushTab(viewer);
207: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lmP->grad);
208: PetscViewerASCIIPopTab(viewer);
209: }
210: return(0);
211: }
215: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
216: {
217: TAO_BLMVM *blm = (TAO_BLMVM *) tao->data;
224: if (!tao->gradient || !blm->unprojected_gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");
226: VecCopy(tao->gradient,DXL);
227: VecAXPY(DXL,-1.0,blm->unprojected_gradient);
228: VecSet(DXU,0.0);
229: VecPointwiseMax(DXL,DXL,DXU);
231: VecCopy(blm->unprojected_gradient,DXU);
232: VecAXPY(DXU,-1.0,tao->gradient);
233: VecAXPY(DXU,1.0,DXL);
234: return(0);
235: }
237: /* ---------------------------------------------------------- */
238: /*MC
239: TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
240: for nonlinear minimization with bound constraints. It is an extension
241: of TAOLMVM
243: Options Database Keys:
244: + -tao_lmm_vectors - number of vectors to use for approximation
245: . -tao_lmm_scale_type - "none","scalar","broyden"
246: . -tao_lmm_limit_type - "none","average","relative","absolute"
247: . -tao_lmm_rescale_type - "none","scalar","gl"
248: . -tao_lmm_limit_mu - mu limiting factor
249: . -tao_lmm_limit_nu - nu limiting factor
250: . -tao_lmm_delta_min - minimum delta value
251: . -tao_lmm_delta_max - maximum delta value
252: . -tao_lmm_broyden_phi - phi factor for Broyden scaling
253: . -tao_lmm_scalar_alpha - alpha factor for scalar scaling
254: . -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal
255: . -tao_lmm_rescale_beta - beta factor for rescaling diagonal
256: . -tao_lmm_scalar_history - amount of history for scalar scaling
257: . -tao_lmm_rescale_history - amount of history for rescaling diagonal
258: - -tao_lmm_eps - rejection tolerance
260: Level: beginner
261: M*/
262: EXTERN_C_BEGIN
265: PetscErrorCode TaoCreate_BLMVM(Tao tao)
266: {
267: TAO_BLMVM *blmP;
268: const char *morethuente_type = TAOLINESEARCHMT;
272: tao->ops->setup = TaoSetup_BLMVM;
273: tao->ops->solve = TaoSolve_BLMVM;
274: tao->ops->view = TaoView_BLMVM;
275: tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
276: tao->ops->destroy = TaoDestroy_BLMVM;
277: tao->ops->computedual = TaoComputeDual_BLMVM;
279: PetscNewLog(tao,&blmP);
280: tao->data = (void*)blmP;
281: tao->max_it = 2000;
282: tao->max_funcs = 4000;
283: tao->fatol = 1e-4;
284: tao->frtol = 1e-4;
286: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
287: TaoLineSearchSetType(tao->linesearch, morethuente_type);
288: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
289: return(0);
290: }
291: EXTERN_C_END