Actual source code: blmvm.c
petsc-3.7.5 2017-01-01
1: #include <petsctaolinesearch.h>
2: #include <../src/tao/matrix/lmvmmat.h>
3: #include <../src/tao/unconstrained/impls/lmvm/lmvm.h>
4: #include <../src/tao/bound/impls/blmvm/blmvm.h>
6: /*------------------------------------------------------------*/
9: static PetscErrorCode TaoSolve_BLMVM(Tao tao)
10: {
11: PetscErrorCode ierr;
12: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
13: TaoConvergedReason reason = TAO_CONTINUE_ITERATING;
14: TaoLineSearchConvergedReason ls_status = TAOLINESEARCH_CONTINUE_ITERATING;
15: PetscReal f, fold, gdx, gnorm;
16: PetscReal stepsize = 1.0,delta;
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);
25: /* Check convergence criteria */
26: TaoComputeObjectiveAndGradient(tao, tao->solution,&f,blmP->unprojected_gradient);
27: VecBoundGradientProjection(blmP->unprojected_gradient,tao->solution, tao->XL,tao->XU,tao->gradient);
29: TaoGradientNorm(tao, tao->gradient,NORM_2,&gnorm);
30: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Inf pr NaN");
32: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize, &reason);
33: if (reason != TAO_CONTINUE_ITERATING) return(0);
35: /* Set initial scaling for the function */
36: if (f != 0.0) {
37: delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm);
38: } else {
39: delta = 2.0 / (gnorm*gnorm);
40: }
41: MatLMVMSetDelta(blmP->M,delta);
43: /* Set counter for gradient/reset steps */
44: blmP->grad = 0;
45: blmP->reset = 0;
47: /* Have not converged; continue with Newton method */
48: while (reason == TAO_CONTINUE_ITERATING) {
49: /* Compute direction */
50: MatLMVMUpdate(blmP->M, tao->solution, tao->gradient);
51: MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
52: VecBoundGradientProjection(tao->stepdirection,tao->solution,tao->XL,tao->XU,tao->gradient);
54: /* Check for success (descent direction) */
55: VecDot(blmP->unprojected_gradient, tao->gradient, &gdx);
56: if (gdx <= 0) {
57: /* Step is not descent or solve was not successful
58: Use steepest descent direction (scaled) */
59: ++blmP->grad;
61: if (f != 0.0) {
62: delta = 2.0*PetscAbsScalar(f) / (gnorm*gnorm);
63: } else {
64: delta = 2.0 / (gnorm*gnorm);
65: }
66: MatLMVMSetDelta(blmP->M,delta);
67: MatLMVMReset(blmP->M);
68: MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
69: MatLMVMSolve(blmP->M,blmP->unprojected_gradient, tao->stepdirection);
70: }
71: VecScale(tao->stepdirection,-1.0);
73: /* Perform the linesearch */
74: fold = f;
75: VecCopy(tao->solution, blmP->Xold);
76: VecCopy(blmP->unprojected_gradient, blmP->Gold);
77: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
78: TaoLineSearchApply(tao->linesearch, tao->solution, &f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);
79: TaoAddLineSearchCounts(tao);
81: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
82: /* Linesearch failed
83: Reset factors and use scaled (projected) gradient step */
84: ++blmP->reset;
86: f = fold;
87: VecCopy(blmP->Xold, tao->solution);
88: VecCopy(blmP->Gold, blmP->unprojected_gradient);
90: if (f != 0.0) {
91: delta = 2.0* PetscAbsScalar(f) / (gnorm*gnorm);
92: } else {
93: delta = 2.0/ (gnorm*gnorm);
94: }
95: MatLMVMSetDelta(blmP->M,delta);
96: MatLMVMReset(blmP->M);
97: MatLMVMUpdate(blmP->M, tao->solution, blmP->unprojected_gradient);
98: MatLMVMSolve(blmP->M, blmP->unprojected_gradient, tao->stepdirection);
99: VecScale(tao->stepdirection, -1.0);
101: /* This may be incorrect; linesearch has values fo stepmax and stepmin
102: that should be reset. */
103: TaoLineSearchSetInitialStepLength(tao->linesearch,1.0);
104: TaoLineSearchApply(tao->linesearch,tao->solution,&f, blmP->unprojected_gradient, tao->stepdirection, &stepsize, &ls_status);
105: TaoAddLineSearchCounts(tao);
107: if (ls_status != TAOLINESEARCH_SUCCESS && ls_status != TAOLINESEARCH_SUCCESS_USER) {
108: tao->reason = TAO_DIVERGED_LS_FAILURE;
109: break;
110: }
111: }
113: /* Check for converged */
114: VecBoundGradientProjection(blmP->unprojected_gradient, tao->solution, tao->XL, tao->XU, tao->gradient);
115: TaoGradientNorm(tao, tao->gradient, NORM_2, &gnorm);
118: if (PetscIsInfOrNanReal(f) || PetscIsInfOrNanReal(gnorm)) SETERRQ(PETSC_COMM_SELF,1, "User provided compute function generated Not-a-Number");
119: tao->niter++;
120: TaoMonitor(tao, tao->niter, f, gnorm, 0.0, stepsize, &reason);
121: }
122: return(0);
123: }
127: static PetscErrorCode TaoSetup_BLMVM(Tao tao)
128: {
129: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
130: PetscInt n,N;
132: KSP H0ksp;
135: /* Existence of tao->solution checked in TaoSetup() */
136: VecDuplicate(tao->solution,&blmP->Xold);
137: VecDuplicate(tao->solution,&blmP->Gold);
138: VecDuplicate(tao->solution, &blmP->unprojected_gradient);
140: if (!tao->stepdirection) {
141: VecDuplicate(tao->solution, &tao->stepdirection);
142: }
143: if (!tao->gradient) {
144: VecDuplicate(tao->solution,&tao->gradient);
145: }
146: if (!tao->XL) {
147: VecDuplicate(tao->solution,&tao->XL);
148: VecSet(tao->XL,PETSC_NINFINITY);
149: }
150: if (!tao->XU) {
151: VecDuplicate(tao->solution,&tao->XU);
152: VecSet(tao->XU,PETSC_INFINITY);
153: }
154: /* Create matrix for the limited memory approximation */
155: VecGetLocalSize(tao->solution,&n);
156: VecGetSize(tao->solution,&N);
157: MatCreateLMVM(((PetscObject)tao)->comm,n,N,&blmP->M);
158: MatLMVMAllocateVectors(blmP->M,tao->solution);
160: /* If the user has set a matrix to solve as the initial H0, set the options prefix here, and set up the KSP */
161: if (blmP->H0) {
162: const char *prefix;
163: PC H0pc;
165: MatLMVMSetH0(blmP->M, blmP->H0);
166: MatLMVMGetH0KSP(blmP->M, &H0ksp);
168: TaoGetOptionsPrefix(tao, &prefix);
169: KSPSetOptionsPrefix(H0ksp, prefix);
170: PetscObjectAppendOptionsPrefix((PetscObject)H0ksp, "tao_h0_");
171: KSPGetPC(H0ksp, &H0pc);
172: PetscObjectAppendOptionsPrefix((PetscObject)H0pc, "tao_h0_");
174: KSPSetFromOptions(H0ksp);
175: KSPSetUp(H0ksp);
176: }
178: return(0);
179: }
181: /* ---------------------------------------------------------- */
184: static PetscErrorCode TaoDestroy_BLMVM(Tao tao)
185: {
186: TAO_BLMVM *blmP = (TAO_BLMVM *)tao->data;
190: if (tao->setupcalled) {
191: MatDestroy(&blmP->M);
192: VecDestroy(&blmP->unprojected_gradient);
193: VecDestroy(&blmP->Xold);
194: VecDestroy(&blmP->Gold);
195: }
197: if (blmP->H0) {
198: PetscObjectDereference((PetscObject)blmP->H0);
199: }
201: PetscFree(tao->data);
202: return(0);
203: }
205: /*------------------------------------------------------------*/
208: static PetscErrorCode TaoSetFromOptions_BLMVM(PetscOptionItems* PetscOptionsObject,Tao tao)
209: {
213: PetscOptionsHead(PetscOptionsObject,"Limited-memory variable-metric method for bound constrained optimization");
214: TaoLineSearchSetFromOptions(tao->linesearch);
215: PetscOptionsTail();
216: return(0);
217: }
220: /*------------------------------------------------------------*/
223: static int TaoView_BLMVM(Tao tao, PetscViewer viewer)
224: {
225: TAO_BLMVM *lmP = (TAO_BLMVM *)tao->data;
226: PetscBool isascii;
230: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii);
231: if (isascii) {
232: PetscViewerASCIIPushTab(viewer);
233: PetscViewerASCIIPrintf(viewer, "Gradient steps: %D\n", lmP->grad);
234: PetscViewerASCIIPopTab(viewer);
235: }
236: return(0);
237: }
241: static PetscErrorCode TaoComputeDual_BLMVM(Tao tao, Vec DXL, Vec DXU)
242: {
243: TAO_BLMVM *blm = (TAO_BLMVM *) tao->data;
250: if (!tao->gradient || !blm->unprojected_gradient) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Dual variables don't exist yet or no longer exist.\n");
252: VecCopy(tao->gradient,DXL);
253: VecAXPY(DXL,-1.0,blm->unprojected_gradient);
254: VecSet(DXU,0.0);
255: VecPointwiseMax(DXL,DXL,DXU);
257: VecCopy(blm->unprojected_gradient,DXU);
258: VecAXPY(DXU,-1.0,tao->gradient);
259: VecAXPY(DXU,1.0,DXL);
260: return(0);
261: }
263: /* ---------------------------------------------------------- */
264: /*MC
265: TAOBLMVM - Bounded limited memory variable metric is a quasi-Newton method
266: for nonlinear minimization with bound constraints. It is an extension
267: of TAOLMVM
269: Options Database Keys:
270: + -tao_lmm_vectors - number of vectors to use for approximation
271: . -tao_lmm_scale_type - "none","scalar","broyden"
272: . -tao_lmm_limit_type - "none","average","relative","absolute"
273: . -tao_lmm_rescale_type - "none","scalar","gl"
274: . -tao_lmm_limit_mu - mu limiting factor
275: . -tao_lmm_limit_nu - nu limiting factor
276: . -tao_lmm_delta_min - minimum delta value
277: . -tao_lmm_delta_max - maximum delta value
278: . -tao_lmm_broyden_phi - phi factor for Broyden scaling
279: . -tao_lmm_scalar_alpha - alpha factor for scalar scaling
280: . -tao_lmm_rescale_alpha - alpha factor for rescaling diagonal
281: . -tao_lmm_rescale_beta - beta factor for rescaling diagonal
282: . -tao_lmm_scalar_history - amount of history for scalar scaling
283: . -tao_lmm_rescale_history - amount of history for rescaling diagonal
284: - -tao_lmm_eps - rejection tolerance
286: Level: beginner
287: M*/
290: PETSC_EXTERN PetscErrorCode TaoCreate_BLMVM(Tao tao)
291: {
292: TAO_BLMVM *blmP;
293: const char *morethuente_type = TAOLINESEARCHMT;
297: tao->ops->setup = TaoSetup_BLMVM;
298: tao->ops->solve = TaoSolve_BLMVM;
299: tao->ops->view = TaoView_BLMVM;
300: tao->ops->setfromoptions = TaoSetFromOptions_BLMVM;
301: tao->ops->destroy = TaoDestroy_BLMVM;
302: tao->ops->computedual = TaoComputeDual_BLMVM;
304: PetscNewLog(tao,&blmP);
305: blmP->H0 = NULL;
306: tao->data = (void*)blmP;
308: /* Override default settings (unless already changed) */
309: if (!tao->max_it_changed) tao->max_it = 2000;
310: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
312: TaoLineSearchCreate(((PetscObject)tao)->comm, &tao->linesearch);
313: TaoLineSearchSetType(tao->linesearch, morethuente_type);
314: TaoLineSearchUseTaoRoutines(tao->linesearch,tao);
315: TaoLineSearchSetOptionsPrefix(tao->linesearch,tao->hdr.prefix);
316: return(0);
317: }
321: PETSC_EXTERN PetscErrorCode TaoLMVMSetH0(Tao tao, Mat H0)
322: {
323: TAO_LMVM *lmP;
324: TAO_BLMVM *blmP;
325: const TaoType type;
326: PetscBool is_lmvm, is_blmvm;
330: TaoGetType(tao, &type);
331: PetscStrcmp(type, TAOLMVM, &is_lmvm);
332: PetscStrcmp(type, TAOBLMVM, &is_blmvm);
334: if (is_lmvm) {
335: lmP = (TAO_LMVM *)tao->data;
336: PetscObjectReference((PetscObject)H0);
337: lmP->H0 = H0;
338: } else if (is_blmvm) {
339: blmP = (TAO_BLMVM *)tao->data;
340: PetscObjectReference((PetscObject)H0);
341: blmP->H0 = H0;
342: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");
344: return(0);
345: }
349: PETSC_EXTERN PetscErrorCode TaoLMVMGetH0(Tao tao, Mat *H0)
350: {
351: TAO_LMVM *lmP;
352: TAO_BLMVM *blmP;
353: const TaoType type;
354: PetscBool is_lmvm, is_blmvm;
355: Mat M;
359: TaoGetType(tao, &type);
360: PetscStrcmp(type, TAOLMVM, &is_lmvm);
361: PetscStrcmp(type, TAOBLMVM, &is_blmvm);
363: if (is_lmvm) {
364: lmP = (TAO_LMVM *)tao->data;
365: M = lmP->M;
366: } else if (is_blmvm) {
367: blmP = (TAO_BLMVM *)tao->data;
368: M = blmP->M;
369: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");
371: MatLMVMGetH0(M, H0);
372: return(0);
373: }
377: PETSC_EXTERN PetscErrorCode TaoLMVMGetH0KSP(Tao tao, KSP *ksp)
378: {
379: TAO_LMVM *lmP;
380: TAO_BLMVM *blmP;
381: const TaoType type;
382: PetscBool is_lmvm, is_blmvm;
383: Mat M;
386: TaoGetType(tao, &type);
387: PetscStrcmp(type, TAOLMVM, &is_lmvm);
388: PetscStrcmp(type, TAOBLMVM, &is_blmvm);
390: if (is_lmvm) {
391: lmP = (TAO_LMVM *)tao->data;
392: M = lmP->M;
393: } else if (is_blmvm) {
394: blmP = (TAO_BLMVM *)tao->data;
395: M = blmP->M;
396: } else SETERRQ(PetscObjectComm((PetscObject)tao), PETSC_ERR_ARG_WRONGSTATE, "This routine applies to TAO_LMVM and TAO_BLMVM.");
398: MatLMVMGetH0KSP(M, ksp);
399: return(0);
400: }