Actual source code: bmrm.c
petsc-3.6.1 2015-08-06
1: #include <../src/tao/unconstrained/impls/bmrm/bmrm.h>
3: static PetscErrorCode init_df_solver(TAO_DF*);
4: static PetscErrorCode ensure_df_space(PetscInt, TAO_DF*);
5: static PetscErrorCode destroy_df_solver(TAO_DF*);
6: static PetscReal phi(PetscReal*,PetscInt,PetscReal,PetscReal*,PetscReal,PetscReal*,PetscReal*,PetscReal*);
7: static PetscInt project(PetscInt,PetscReal*,PetscReal,PetscReal*,PetscReal*,PetscReal*,PetscReal*,PetscReal*,TAO_DF*);
8: static PetscErrorCode solve(TAO_DF*);
11: /*------------------------------------------------------------*/
12: /* The main solver function
14: f = Remp(W) This is what the user provides us from the application layer
15: So the ComputeGradient function for instance should get us back the subgradient of Remp(W)
17: Regularizer assumed to be L2 norm = lambda*0.5*W'W ()
18: */
22: static PetscErrorCode make_grad_node(Vec X, Vec_Chain **p)
23: {
27: PetscNew(p);
28: VecDuplicate(X, &(*p)->V);
29: VecCopy(X, (*p)->V);
30: (*p)->next = NULL;
31: return(0);
32: }
36: static PetscErrorCode destroy_grad_list(Vec_Chain *head)
37: {
39: Vec_Chain *p = head->next, *q;
42: while(p) {
43: q = p->next;
44: VecDestroy(&p->V);
45: PetscFree(p);
46: p = q;
47: }
48: head->next = NULL;
49: return(0);
50: }
55: static PetscErrorCode TaoSolve_BMRM(Tao tao)
56: {
57: PetscErrorCode ierr;
58: TaoConvergedReason reason;
59: TAO_DF df;
60: TAO_BMRM *bmrm = (TAO_BMRM*)tao->data;
62: /* Values and pointers to parts of the optimization problem */
63: PetscReal f = 0.0;
64: Vec W = tao->solution;
65: Vec G = tao->gradient;
66: PetscReal lambda;
67: PetscReal bt;
68: Vec_Chain grad_list, *tail_glist, *pgrad;
69: PetscInt i;
70: PetscMPIInt rank;
72: /* Used in converged criteria check */
73: PetscReal reg;
74: PetscReal jtwt, max_jtwt, pre_epsilon, epsilon, jw, min_jw;
75: PetscReal innerSolverTol;
76: MPI_Comm comm;
79: PetscObjectGetComm((PetscObject)tao,&comm);
80: MPI_Comm_rank(comm, &rank);
81: lambda = bmrm->lambda;
83: /* Check Stopping Condition */
84: tao->step = 1.0;
85: max_jtwt = -BMRM_INFTY;
86: min_jw = BMRM_INFTY;
87: innerSolverTol = 1.0;
88: epsilon = 0.0;
90: if (!rank) {
91: init_df_solver(&df);
92: grad_list.next = NULL;
93: tail_glist = &grad_list;
94: }
96: df.tol = 1e-6;
97: reason = TAO_CONTINUE_ITERATING;
99: /*-----------------Algorithm Begins------------------------*/
100: /* make the scatter */
101: VecScatterCreateToZero(W, &bmrm->scatter, &bmrm->local_w);
102: VecAssemblyBegin(bmrm->local_w);
103: VecAssemblyEnd(bmrm->local_w);
105: /* NOTE: In application pass the sub-gradient of Remp(W) */
106: TaoComputeObjectiveAndGradient(tao, W, &f, G);
107: TaoMonitor(tao,tao->niter,f,1.0,0.0,tao->step,&reason);
108: while (reason == TAO_CONTINUE_ITERATING) {
109: /* compute bt = Remp(Wt-1) - <Wt-1, At> */
110: VecDot(W, G, &bt);
111: bt = f - bt;
113: /* First gather the gradient to the master node */
114: VecScatterBegin(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
115: VecScatterEnd(bmrm->scatter, G, bmrm->local_w, INSERT_VALUES, SCATTER_FORWARD);
117: /* Bring up the inner solver */
118: if (!rank) {
119: ensure_df_space(tao->niter+1, &df);
120: make_grad_node(bmrm->local_w, &pgrad);
121: tail_glist->next = pgrad;
122: tail_glist = pgrad;
124: df.a[tao->niter] = 1.0;
125: df.f[tao->niter] = -bt;
126: df.u[tao->niter] = 1.0;
127: df.l[tao->niter] = 0.0;
129: /* set up the Q */
130: pgrad = grad_list.next;
131: for (i=0; i<=tao->niter; i++) {
132: VecDot(pgrad->V, bmrm->local_w, ®);
133: df.Q[i][tao->niter] = df.Q[tao->niter][i] = reg / lambda;
134: pgrad = pgrad->next;
135: }
137: if (tao->niter > 0) {
138: df.x[tao->niter] = 0.0;
139: solve(&df);
140: } else
141: df.x[0] = 1.0;
143: /* now computing Jt*(alpha_t) which should be = Jt(wt) to check convergence */
144: jtwt = 0.0;
145: VecSet(bmrm->local_w, 0.0);
146: pgrad = grad_list.next;
147: for (i=0; i<=tao->niter; i++) {
148: jtwt -= df.x[i] * df.f[i];
149: VecAXPY(bmrm->local_w, -df.x[i] / lambda, pgrad->V);
150: pgrad = pgrad->next;
151: }
153: VecNorm(bmrm->local_w, NORM_2, ®);
154: reg = 0.5*lambda*reg*reg;
155: jtwt -= reg;
156: } /* end if rank == 0 */
158: /* scatter the new W to all nodes */
159: VecScatterBegin(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);
160: VecScatterEnd(bmrm->scatter,bmrm->local_w,W,INSERT_VALUES,SCATTER_REVERSE);
162: TaoComputeObjectiveAndGradient(tao, W, &f, G);
164: MPI_Bcast(&jtwt,1,MPIU_REAL,0,comm);
165: MPI_Bcast(®,1,MPIU_REAL,0,comm);
167: jw = reg + f; /* J(w) = regularizer + Remp(w) */
168: if (jw < min_jw) min_jw = jw;
169: if (jtwt > max_jtwt) max_jtwt = jtwt;
171: pre_epsilon = epsilon;
172: epsilon = min_jw - jtwt;
174: if (!rank) {
175: if (innerSolverTol > epsilon) innerSolverTol = epsilon;
176: else if (innerSolverTol < 1e-7) innerSolverTol = 1e-7;
178: /* if the annealing doesn't work well, lower the inner solver tolerance */
179: if(pre_epsilon < epsilon) innerSolverTol *= 0.2;
181: df.tol = innerSolverTol*0.5;
182: }
184: tao->niter++;
185: TaoMonitor(tao,tao->niter,min_jw,epsilon,0.0,tao->step,&reason);
186: }
188: /* free all the memory */
189: if (!rank) {
190: destroy_grad_list(&grad_list);
191: destroy_df_solver(&df);
192: }
194: VecDestroy(&bmrm->local_w);
195: VecScatterDestroy(&bmrm->scatter);
196: return(0);
197: }
200: /* ---------------------------------------------------------- */
204: static PetscErrorCode TaoSetup_BMRM(Tao tao)
205: {
210: /* Allocate some arrays */
211: if (!tao->gradient) {
212: VecDuplicate(tao->solution, &tao->gradient);
213: }
214: return(0);
215: }
217: /*------------------------------------------------------------*/
220: static PetscErrorCode TaoDestroy_BMRM(Tao tao)
221: {
225: PetscFree(tao->data);
226: return(0);
227: }
231: static PetscErrorCode TaoSetFromOptions_BMRM(PetscOptions *PetscOptionsObject,Tao tao)
232: {
234: TAO_BMRM* bmrm = (TAO_BMRM*)tao->data;
237: PetscOptionsHead(PetscOptionsObject,"BMRM for regularized risk minimization");
238: PetscOptionsReal("-tao_bmrm_lambda", "regulariser weight","", 100,&bmrm->lambda,NULL);
239: PetscOptionsTail();
240: return(0);
241: }
243: /*------------------------------------------------------------*/
246: static PetscErrorCode TaoView_BMRM(Tao tao, PetscViewer viewer)
247: {
248: PetscBool isascii;
252: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
253: if (isascii) {
254: PetscViewerASCIIPushTab(viewer);
255: PetscViewerASCIIPopTab(viewer);
256: }
257: return(0);
258: }
260: /*------------------------------------------------------------*/
261: /*MC
262: TAOBMRM - bundle method for regularized risk minimization
264: Options Database Keys:
265: . - tao_bmrm_lambda - regulariser weight
267: Level: beginner
268: M*/
272: PETSC_EXTERN PetscErrorCode TaoCreate_BMRM(Tao tao)
273: {
274: TAO_BMRM *bmrm;
278: tao->ops->setup = TaoSetup_BMRM;
279: tao->ops->solve = TaoSolve_BMRM;
280: tao->ops->view = TaoView_BMRM;
281: tao->ops->setfromoptions = TaoSetFromOptions_BMRM;
282: tao->ops->destroy = TaoDestroy_BMRM;
284: PetscNewLog(tao,&bmrm);
285: bmrm->lambda = 1.0;
286: tao->data = (void*)bmrm;
288: /* Override default settings (unless already changed) */
289: if (!tao->max_it_changed) tao->max_it = 2000;
290: if (!tao->max_funcs_changed) tao->max_funcs = 4000;
291: if (!tao->fatol_changed) tao->fatol = 1.0e-12;
292: if (!tao->frtol_changed) tao->frtol = 1.0e-12;
293: if (!tao->gatol_changed) tao->gatol = 1.0e-12;
294: if (!tao->grtol_changed) tao->grtol = 1.0e-12;
296: return(0);
297: }
301: PetscErrorCode init_df_solver(TAO_DF *df)
302: {
303: PetscInt i, n = INCRE_DIM;
307: /* default values */
308: df->maxProjIter = 200;
309: df->maxPGMIter = 300000;
310: df->b = 1.0;
312: /* memory space required by Dai-Fletcher */
313: df->cur_num_cp = n;
314: PetscMalloc1(n, &df->f);
315: PetscMalloc1(n, &df->a);
316: PetscMalloc1(n, &df->l);
317: PetscMalloc1(n, &df->u);
318: PetscMalloc1(n, &df->x);
319: PetscMalloc1(n, &df->Q);
321: for (i = 0; i < n; i ++) {
322: PetscMalloc1(n, &df->Q[i]);
323: }
325: PetscMalloc1(n, &df->g);
326: PetscMalloc1(n, &df->y);
327: PetscMalloc1(n, &df->tempv);
328: PetscMalloc1(n, &df->d);
329: PetscMalloc1(n, &df->Qd);
330: PetscMalloc1(n, &df->t);
331: PetscMalloc1(n, &df->xplus);
332: PetscMalloc1(n, &df->tplus);
333: PetscMalloc1(n, &df->sk);
334: PetscMalloc1(n, &df->yk);
336: PetscMalloc1(n, &df->ipt);
337: PetscMalloc1(n, &df->ipt2);
338: PetscMalloc1(n, &df->uv);
339: return(0);
340: }
344: PetscErrorCode ensure_df_space(PetscInt dim, TAO_DF *df)
345: {
347: PetscReal *tmp, **tmp_Q;
348: PetscInt i, n, old_n;
351: df->dim = dim;
352: if (dim <= df->cur_num_cp) return(0);
354: old_n = df->cur_num_cp;
355: df->cur_num_cp += INCRE_DIM;
356: n = df->cur_num_cp;
358: /* memory space required by dai-fletcher */
359: PetscMalloc1(n, &tmp);
360: PetscMemcpy(tmp, df->f, sizeof(PetscReal)*old_n);
361: PetscFree(df->f);
362: df->f = tmp;
364: PetscMalloc1(n, &tmp);
365: PetscMemcpy(tmp, df->a, sizeof(PetscReal)*old_n);
366: PetscFree(df->a);
367: df->a = tmp;
369: PetscMalloc1(n, &tmp);
370: PetscMemcpy(tmp, df->l, sizeof(PetscReal)*old_n);
371: PetscFree(df->l);
372: df->l = tmp;
374: PetscMalloc1(n, &tmp);
375: PetscMemcpy(tmp, df->u, sizeof(PetscReal)*old_n);
376: PetscFree(df->u);
377: df->u = tmp;
379: PetscMalloc1(n, &tmp);
380: PetscMemcpy(tmp, df->x, sizeof(PetscReal)*old_n);
381: PetscFree(df->x);
382: df->x = tmp;
384: PetscMalloc1(n, &tmp_Q);
385: for (i = 0; i < n; i ++) {
386: PetscMalloc1(n, &tmp_Q[i]);
387: if (i < old_n) {
388: PetscMemcpy(tmp_Q[i], df->Q[i], sizeof(PetscReal)*old_n);
389: PetscFree(df->Q[i]);
390: }
391: }
393: PetscFree(df->Q);
394: df->Q = tmp_Q;
396: PetscFree(df->g);
397: PetscMalloc1(n, &df->g);
399: PetscFree(df->y);
400: PetscMalloc1(n, &df->y);
402: PetscFree(df->tempv);
403: PetscMalloc1(n, &df->tempv);
405: PetscFree(df->d);
406: PetscMalloc1(n, &df->d);
408: PetscFree(df->Qd);
409: PetscMalloc1(n, &df->Qd);
411: PetscFree(df->t);
412: PetscMalloc1(n, &df->t);
414: PetscFree(df->xplus);
415: PetscMalloc1(n, &df->xplus);
417: PetscFree(df->tplus);
418: PetscMalloc1(n, &df->tplus);
420: PetscFree(df->sk);
421: PetscMalloc1(n, &df->sk);
423: PetscFree(df->yk);
424: PetscMalloc1(n, &df->yk);
426: PetscFree(df->ipt);
427: PetscMalloc1(n, &df->ipt);
429: PetscFree(df->ipt2);
430: PetscMalloc1(n, &df->ipt2);
432: PetscFree(df->uv);
433: PetscMalloc1(n, &df->uv);
434: return(0);
435: }
439: PetscErrorCode destroy_df_solver(TAO_DF *df)
440: {
442: PetscInt i;
445: PetscFree(df->f);
446: PetscFree(df->a);
447: PetscFree(df->l);
448: PetscFree(df->u);
449: PetscFree(df->x);
451: for (i = 0; i < df->cur_num_cp; i ++) {
452: PetscFree(df->Q[i]);
453: }
454: PetscFree(df->Q);
455: PetscFree(df->ipt);
456: PetscFree(df->ipt2);
457: PetscFree(df->uv);
458: PetscFree(df->g);
459: PetscFree(df->y);
460: PetscFree(df->tempv);
461: PetscFree(df->d);
462: PetscFree(df->Qd);
463: PetscFree(df->t);
464: PetscFree(df->xplus);
465: PetscFree(df->tplus);
466: PetscFree(df->sk);
467: PetscFree(df->yk);
468: return(0);
469: }
471: /* Piecewise linear monotone target function for the Dai-Fletcher projector */
474: PetscReal phi(PetscReal *x,PetscInt n,PetscReal lambda,PetscReal *a,PetscReal b,PetscReal *c,PetscReal *l,PetscReal *u)
475: {
476: PetscReal r = 0.0;
477: PetscInt i;
479: for (i = 0; i < n; i++){
480: x[i] = -c[i] + lambda*a[i];
481: if (x[i] > u[i]) x[i] = u[i];
482: else if(x[i] < l[i]) x[i] = l[i];
483: r += a[i]*x[i];
484: }
485: return r - b;
486: }
488: /** Modified Dai-Fletcher QP projector solves the problem:
489: *
490: * minimise 0.5*x'*x - c'*x
491: * subj to a'*x = b
492: * l \leq x \leq u
493: *
494: * \param c The point to be projected onto feasible set
495: */
498: PetscInt project(PetscInt n,PetscReal *a,PetscReal b,PetscReal *c,PetscReal *l,PetscReal *u,PetscReal *x,PetscReal *lam_ext,TAO_DF *df)
499: {
500: PetscReal lambda, lambdal, lambdau, dlambda, lambda_new;
501: PetscReal r, rl, ru, s;
502: PetscInt innerIter;
503: PetscBool nonNegativeSlack = PETSC_FALSE;
506: *lam_ext = 0;
507: lambda = 0;
508: dlambda = 0.5;
509: innerIter = 1;
511: /* \phi(x;lambda) := 0.5*x'*x + c'*x - lambda*(a'*x-b)
512: *
513: * Optimality conditions for \phi:
514: *
515: * 1. lambda <= 0
516: * 2. r <= 0
517: * 3. r*lambda == 0
518: */
520: /* Bracketing Phase */
521: r = phi(x, n, lambda, a, b, c, l, u);
523: if(nonNegativeSlack) {
524: /* inequality constraint, i.e., with \xi >= 0 constraint */
525: if (r < TOL_R) return 0;
526: } else {
527: /* equality constraint ,i.e., without \xi >= 0 constraint */
528: if (fabs(r) < TOL_R) return 0;
529: }
531: if (r < 0.0){
532: lambdal = lambda;
533: rl = r;
534: lambda = lambda + dlambda;
535: r = phi(x, n, lambda, a, b, c, l, u);
536: while (r < 0.0 && dlambda < BMRM_INFTY) {
537: lambdal = lambda;
538: s = rl/r - 1.0;
539: if (s < 0.1) s = 0.1;
540: dlambda = dlambda + dlambda/s;
541: lambda = lambda + dlambda;
542: rl = r;
543: r = phi(x, n, lambda, a, b, c, l, u);
544: }
545: lambdau = lambda;
546: ru = r;
547: } else {
548: lambdau = lambda;
549: ru = r;
550: lambda = lambda - dlambda;
551: r = phi(x, n, lambda, a, b, c, l, u);
552: while (r > 0.0 && dlambda > -BMRM_INFTY) {
553: lambdau = lambda;
554: s = ru/r - 1.0;
555: if (s < 0.1) s = 0.1;
556: dlambda = dlambda + dlambda/s;
557: lambda = lambda - dlambda;
558: ru = r;
559: r = phi(x, n, lambda, a, b, c, l, u);
560: }
561: lambdal = lambda;
562: rl = r;
563: }
565: if(fabs(dlambda) > BMRM_INFTY) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"L2N2_DaiFletcherPGM detected Infeasible QP problem!");
567: if(ru == 0){
568: lambda = lambdau;
569: r = phi(x, n, lambda, a, b, c, l, u);
570: return innerIter;
571: }
573: /* Secant Phase */
574: s = 1.0 - rl/ru;
575: dlambda = dlambda/s;
576: lambda = lambdau - dlambda;
577: r = phi(x, n, lambda, a, b, c, l, u);
579: while (fabs(r) > TOL_R
580: && dlambda > TOL_LAM * (1.0 + fabs(lambda))
581: && innerIter < df->maxProjIter){
582: innerIter++;
583: if (r > 0.0){
584: if (s <= 2.0){
585: lambdau = lambda;
586: ru = r;
587: s = 1.0 - rl/ru;
588: dlambda = (lambdau - lambdal) / s;
589: lambda = lambdau - dlambda;
590: } else {
591: s = ru/r-1.0;
592: if (s < 0.1) s = 0.1;
593: dlambda = (lambdau - lambda) / s;
594: lambda_new = 0.75*lambdal + 0.25*lambda;
595: if (lambda_new < (lambda - dlambda))
596: lambda_new = lambda - dlambda;
597: lambdau = lambda;
598: ru = r;
599: lambda = lambda_new;
600: s = (lambdau - lambdal) / (lambdau - lambda);
601: }
602: } else {
603: if (s >= 2.0){
604: lambdal = lambda;
605: rl = r;
606: s = 1.0 - rl/ru;
607: dlambda = (lambdau - lambdal) / s;
608: lambda = lambdau - dlambda;
609: } else {
610: s = rl/r - 1.0;
611: if (s < 0.1) s = 0.1;
612: dlambda = (lambda-lambdal) / s;
613: lambda_new = 0.75*lambdau + 0.25*lambda;
614: if (lambda_new > (lambda + dlambda))
615: lambda_new = lambda + dlambda;
616: lambdal = lambda;
617: rl = r;
618: lambda = lambda_new;
619: s = (lambdau - lambdal) / (lambdau-lambda);
620: }
621: }
622: r = phi(x, n, lambda, a, b, c, l, u);
623: }
625: *lam_ext = lambda;
626: if(innerIter >= df->maxProjIter) {
627: PetscPrintf(PETSC_COMM_SELF, "WARNING: DaiFletcher max iterations\n");
628: }
629: return innerIter;
630: }
635: PetscErrorCode solve(TAO_DF *df)
636: {
638: PetscInt i, j, innerIter, it, it2, luv, info, lscount = 0, projcount = 0;
639: PetscReal gd, max, ak, bk, akold, bkold, lamnew, alpha, kktlam=0.0, lam_ext;
640: PetscReal DELTAsv, ProdDELTAsv;
641: PetscReal c, *tempQ;
642: PetscReal *x = df->x, *a = df->a, b = df->b, *l = df->l, *u = df->u, tol = df->tol;
643: PetscReal *tempv = df->tempv, *y = df->y, *g = df->g, *d = df->d, *Qd = df->Qd;
644: PetscReal *xplus = df->xplus, *tplus = df->tplus, *sk = df->sk, *yk = df->yk;
645: PetscReal **Q = df->Q, *f = df->f, *t = df->t;
646: PetscInt dim = df->dim, *ipt = df->ipt, *ipt2 = df->ipt2, *uv = df->uv;
648: /*** variables for the adaptive nonmonotone linesearch ***/
649: PetscInt L, llast;
650: PetscReal fr, fbest, fv, fc, fv0;
652: c = BMRM_INFTY;
654: DELTAsv = EPS_SV;
655: if (tol <= 1.0e-5 || dim <= 20) ProdDELTAsv = 0.0F;
656: else ProdDELTAsv = EPS_SV;
658: for (i = 0; i < dim; i++) tempv[i] = -x[i];
660: lam_ext = 0.0;
662: /* Project the initial solution */
663: projcount += project(dim, a, b, tempv, l, u, x, &lam_ext, df);
665: /* Compute gradient
666: g = Q*x + f; */
668: it = 0;
669: for (i = 0; i < dim; i++) {
670: if (fabs(x[i]) > ProdDELTAsv) ipt[it++] = i;
671: }
673: PetscMemzero(t, dim*sizeof(PetscReal));
674: for (i = 0; i < it; i++){
675: tempQ = Q[ipt[i]];
676: for (j = 0; j < dim; j++) t[j] += (tempQ[j]*x[ipt[i]]);
677: }
678: for (i = 0; i < dim; i++){
679: g[i] = t[i] + f[i];
680: }
683: /* y = -(x_{k} - g_{k}) */
684: for (i = 0; i < dim; i++){
685: y[i] = g[i] - x[i];
686: }
688: /* Project x_{k} - g_{k} */
689: projcount += project(dim, a, b, y, l, u, tempv, &lam_ext, df);
691: /* y = P(x_{k} - g_{k}) - x_{k} */
692: max = ALPHA_MIN;
693: for (i = 0; i < dim; i++){
694: y[i] = tempv[i] - x[i];
695: if (fabs(y[i]) > max) max = fabs(y[i]);
696: }
698: if (max < tol*1e-3){
699: lscount = 0;
700: innerIter = 0;
701: return 0;
702: }
704: alpha = 1.0 / max;
706: /* fv0 = f(x_{0}). Recall t = Q x_{k} */
707: fv0 = 0.0;
708: for (i = 0; i < dim; i++) fv0 += x[i] * (0.5*t[i] + f[i]);
710: /*** adaptive nonmonotone linesearch ***/
711: L = 2;
712: fr = ALPHA_MAX;
713: fbest = fv0;
714: fc = fv0;
715: llast = 0;
716: akold = bkold = 0.0;
718: /*** Iterator begins ***/
719: for (innerIter = 1; innerIter <= df->maxPGMIter; innerIter++) {
721: /* tempv = -(x_{k} - alpha*g_{k}) */
722: for (i = 0; i < dim; i++) tempv[i] = alpha*g[i] - x[i];
724: /* Project x_{k} - alpha*g_{k} */
725: projcount += project(dim, a, b, tempv, l, u, y, &lam_ext, df);
728: /* gd = \inner{d_{k}}{g_{k}}
729: d = P(x_{k} - alpha*g_{k}) - x_{k}
730: */
731: gd = 0.0;
732: for (i = 0; i < dim; i++){
733: d[i] = y[i] - x[i];
734: gd += d[i] * g[i];
735: }
737: /* Gradient computation */
739: /* compute Qd = Q*d or Qd = Q*y - t depending on their sparsity */
741: it = it2 = 0;
742: for (i = 0; i < dim; i++){
743: if (fabs(d[i]) > (ProdDELTAsv*1.0e-2)) ipt[it++] = i;
744: }
745: for (i = 0; i < dim; i++) {
746: if (fabs(y[i]) > ProdDELTAsv) ipt2[it2++] = i;
747: }
749: PetscMemzero(Qd, dim*sizeof(PetscReal));
750: /* compute Qd = Q*d */
751: if (it < it2){
752: for (i = 0; i < it; i++){
753: tempQ = Q[ipt[i]];
754: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * d[ipt[i]]);
755: }
756: } else { /* compute Qd = Q*y-t */
757: for (i = 0; i < it2; i++){
758: tempQ = Q[ipt2[i]];
759: for (j = 0; j < dim; j++) Qd[j] += (tempQ[j] * y[ipt2[i]]);
760: }
761: for (j = 0; j < dim; j++) Qd[j] -= t[j];
762: }
764: /* ak = inner{d_{k}}{d_{k}} */
765: ak = 0.0;
766: for (i = 0; i < dim; i++) ak += d[i] * d[i];
768: bk = 0.0;
769: for (i = 0; i < dim; i++) bk += d[i]*Qd[i];
771: if (bk > EPS*ak && gd < 0.0) lamnew = -gd/bk;
772: else lamnew = 1.0;
774: /* fv is computing f(x_{k} + d_{k}) */
775: fv = 0.0;
776: for (i = 0; i < dim; i++){
777: xplus[i] = x[i] + d[i];
778: tplus[i] = t[i] + Qd[i];
779: fv += xplus[i] * (0.5*tplus[i] + f[i]);
780: }
782: /* fr is fref */
783: if ((innerIter == 1 && fv >= fv0) || (innerIter > 1 && fv >= fr)){
784: lscount++;
785: fv = 0.0;
786: for (i = 0; i < dim; i++){
787: xplus[i] = x[i] + lamnew*d[i];
788: tplus[i] = t[i] + lamnew*Qd[i];
789: fv += xplus[i] * (0.5*tplus[i] + f[i]);
790: }
791: }
793: for (i = 0; i < dim; i++){
794: sk[i] = xplus[i] - x[i];
795: yk[i] = tplus[i] - t[i];
796: x[i] = xplus[i];
797: t[i] = tplus[i];
798: g[i] = t[i] + f[i];
799: }
801: /* update the line search control parameters */
802: if (fv < fbest){
803: fbest = fv;
804: fc = fv;
805: llast = 0;
806: } else {
807: fc = (fc > fv ? fc : fv);
808: llast++;
809: if (llast == L){
810: fr = fc;
811: fc = fv;
812: llast = 0;
813: }
814: }
816: ak = bk = 0.0;
817: for (i = 0; i < dim; i++){
818: ak += sk[i] * sk[i];
819: bk += sk[i] * yk[i];
820: }
822: if (bk <= EPS*ak) alpha = ALPHA_MAX;
823: else {
824: if (bkold < EPS*akold) alpha = ak/bk;
825: else alpha = (akold+ak)/(bkold+bk);
827: if (alpha > ALPHA_MAX) alpha = ALPHA_MAX;
828: else if (alpha < ALPHA_MIN) alpha = ALPHA_MIN;
829: }
831: akold = ak;
832: bkold = bk;
834: /*** stopping criterion based on KKT conditions ***/
835: /* at optimal, gradient of lagrangian w.r.t. x is zero */
837: bk = 0.0;
838: for (i = 0; i < dim; i++) bk += x[i] * x[i];
840: if (PetscSqrtReal(ak) < tol*10 * PetscSqrtReal(bk)){
841: it = 0;
842: luv = 0;
843: kktlam = 0.0;
844: for (i = 0; i < dim; i++){
845: /* x[i] is active hence lagrange multipliers for box constraints
846: are zero. The lagrange multiplier for ineq. const. is then
847: defined as below
848: */
849: if ((x[i] > DELTAsv) && (x[i] < c-DELTAsv)){
850: ipt[it++] = i;
851: kktlam = kktlam - a[i]*g[i];
852: } else uv[luv++] = i;
853: }
855: if (it == 0 && PetscSqrtReal(ak) < tol*0.5 * PetscSqrtReal(bk)) return 0;
856: else {
857: kktlam = kktlam/it;
858: info = 1;
859: for (i = 0; i < it; i++) {
860: if (fabs(a[ipt[i]] * g[ipt[i]] + kktlam) > tol) {
861: info = 0;
862: break;
863: }
864: }
865: if (info == 1) {
866: for (i = 0; i < luv; i++) {
867: if (x[uv[i]] <= DELTAsv){
868: /* x[i] == lower bound, hence, lagrange multiplier (say, beta) for lower bound may
869: not be zero. So, the gradient without beta is > 0
870: */
871: if (g[uv[i]] + kktlam*a[uv[i]] < -tol){
872: info = 0;
873: break;
874: }
875: } else {
876: /* x[i] == upper bound, hence, lagrange multiplier (say, eta) for upper bound may
877: not be zero. So, the gradient without eta is < 0
878: */
879: if (g[uv[i]] + kktlam*a[uv[i]] > tol) {
880: info = 0;
881: break;
882: }
883: }
884: }
885: }
887: if (info == 1) return 0;
888: }
889: }
890: }
891: return 0;
892: }