Actual source code: lmvmmat.c
1: #include "lmvmmat.h" /*I "lmvmmat.h" */
2: #include "tao_util.h" /*I "tao_util.h" */
3: #include "petscksp.h"
5: #define TaoMid(a,b,c) (((a) < (b)) ? \
6: (((b) < (c)) ? (b) : \
7: (((a) < (c)) ? (c) : (a))) : \
8: (((a) < (c)) ? (a) : \
9: (((b) < (c)) ? (c) : (b))))
11: /* These lists are used for setting options */
12: static const char *Scale_Table[64] = {
13: "none","scalar","broyden"
14: };
16: static const char *Rescale_Table[64] = {
17: "none","scalar","gl"
18: };
20: static const char *Limit_Table[64] = {
21: "none","average","relative","absolute"
22: };
27: /*@C
28: MatCreateLMVM - Creates a limited memory matrix for lmvm algorithms.
30: Collective on A
32: Input Parameters:
33: + comm - MPI Communicator
34: . n - local size of vectors
35: - N - global size of vectors
37: Output Parameters:
38: . A - New LMVM matrix
40: Level: developer
42: @*/
44: {
45: MatLMVMCtx *ctx;
47: PetscInt nhistory;
51: /* create data structure and populate with default values */
52: PetscMalloc(sizeof(MatLMVMCtx),(void**)&ctx);
53: ctx->lm=5;
54: ctx->eps=0.0;
55: ctx->limitType=MatLMVM_Limit_None;
56: ctx->scaleType=MatLMVM_Scale_Broyden;
57: ctx->rScaleType = MatLMVM_Rescale_Scalar;
58: ctx->s_alpha = 1.0;
59: ctx->r_alpha = 1.0;
60: ctx->r_beta = 0.5;
61: ctx->mu = 1.0;
62: ctx->nu = 100.0;
64: ctx->phi = 0.125;
66: ctx->scalar_history = 1;
67: ctx->rescale_history = 1;
69: ctx->delta_min = 1e-7;
70: ctx->delta_max = 100.0;
72: /* Begin configuration */
73: PetscOptionsInt("-tao_lmm_vectors", "vectors to use for approximation", "", ctx->lm, &ctx->lm, 0);
74: PetscOptionsReal("-tao_lmm_limit_mu", "mu limiting factor", "", ctx->mu, &ctx->mu, 0);
75: PetscOptionsReal("-tao_lmm_limit_nu", "nu limiting factor", "", ctx->nu, &ctx->nu, 0);
76: PetscOptionsReal("-tao_lmm_broyden_phi", "phi factor for Broyden scaling", "", ctx->phi, &ctx->phi, 0);
77: PetscOptionsReal("-tao_lmm_scalar_alpha", "alpha factor for scalar scaling", "",ctx->s_alpha, &ctx->s_alpha, 0);
78: PetscOptionsReal("-tao_lmm_rescale_alpha", "alpha factor for rescaling diagonal", "", ctx->r_alpha, &ctx->r_alpha, 0);
79: PetscOptionsReal("-tao_lmm_rescale_beta", "beta factor for rescaling diagonal", "", ctx->r_beta, &ctx->r_beta, 0);
80: PetscOptionsInt("-tao_lmm_scalar_history", "amount of history for scalar scaling", "", ctx->scalar_history, &ctx->scalar_history, 0);
81: PetscOptionsInt("-tao_lmm_rescale_history", "amount of history for rescaling diagonal", "", ctx->rescale_history, &ctx->rescale_history, 0);
82: PetscOptionsReal("-tao_lmm_eps", "rejection tolerance", "", ctx->eps, &ctx->eps, 0);
83: PetscOptionsEList("-tao_lmm_scale_type", "scale type", "", Scale_Table, MatLMVM_Scale_Types, Scale_Table[ctx->scaleType], &ctx->scaleType, 0);
84: PetscOptionsEList("-tao_lmm_rescale_type", "rescale type", "", Rescale_Table, MatLMVM_Rescale_Types, Rescale_Table[ctx->rScaleType], &ctx->rScaleType, 0);
85: PetscOptionsEList("-tao_lmm_limit_type", "limit type", "", Limit_Table, MatLMVM_Limit_Types, Limit_Table[ctx->limitType], &ctx->limitType, 0);
86: PetscOptionsReal("-tao_lmm_delta_min", "minimum delta value", "", ctx->delta_min, &ctx->delta_min, 0);
87: PetscOptionsReal("-tao_lmm_delta_max", "maximum delta value", "", ctx->delta_max, &ctx->delta_max, 0);
89: /* Complete configuration */
90: ctx->rescale_history = PetscMin(ctx->rescale_history, ctx->lm);
93: PetscMalloc((ctx->lm+1)*sizeof(PetscReal),(void**)&ctx->rho);
94:
95: PetscMalloc((ctx->lm+1)*sizeof(PetscReal),(void**)&ctx->beta);
96:
98: nhistory = PetscMax(ctx->scalar_history,1);
99: PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->yy_history);
100:
101: PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->ys_history);
102:
103: PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->ss_history);
104:
106: nhistory = PetscMax(ctx->rescale_history,1);
107: PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->yy_rhistory);
108:
109: PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->ys_rhistory);
110:
111: PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->ss_rhistory);
112:
115: /* Finish initializations */
116: ctx->lmnow = 0;
117: ctx->iter = 0;
118: ctx->nupdates = 0;
119: ctx->nrejects = 0;
120: ctx->delta = 1.0;
122: ctx->Gprev = 0;
123: ctx->Xprev = 0;
125: ctx->scale = 0;
126: ctx->useScale = PETSC_FALSE;
128: ctx->H0 = 0;
129: ctx->useDefaultH0=PETSC_TRUE;
130:
131: MatCreateShell(comm, n, n, N, N, ctx, A);
132: MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
133:
134: MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
135:
137: return(0);
138: }
142:
143:
147: {
148: PetscReal sq, yq, dd;
149: PetscInt ll;
150: PetscBool scaled;
151: MatLMVMCtx *shell;
158: MatShellGetContext(A,(void**)&shell);
159: if (shell->lmnow < 1) {
160: shell->rho[0] = 1.0;
161: }
163: VecCopy(b,x);
164: for (ll = 0; ll < shell->lmnow; ++ll) {
165: VecDot(x,shell->S[ll],&sq);
166: shell->beta[ll] = sq * shell->rho[ll];
167: VecAXPY(x,-shell->beta[ll],shell->Y[ll]);
168: }
170: scaled = PETSC_FALSE;
171: if (!scaled && !shell->useDefaultH0 && shell->H0) {
172: MatSolve(shell->H0,x,shell->U);
173: VecDot(x,shell->U,&dd);
174: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
175: /* Accept Hessian solve */
176: VecCopy(shell->U,x);
177: scaled = PETSC_TRUE;
178: }
179: }
181: if (!scaled && shell->useScale) {
182: VecPointwiseMult(shell->U,x,shell->scale);
183: VecDot(x,shell->U,&dd);
184: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
185: /* Accept scaling */
186: VecCopy(shell->U,x);
187: scaled = PETSC_TRUE;
188: }
189: }
190:
191: if (!scaled) {
192: switch(shell->scaleType) {
193: case MatLMVM_Scale_None:
194: break;
196: case MatLMVM_Scale_Scalar:
197: VecScale(x,shell->sigma);
198: break;
199:
200: case MatLMVM_Scale_Broyden:
201: VecPointwiseMult(x,x,shell->D);
202: break;
203: }
204: }
206: for (ll = shell->lmnow-1; ll >= 0; --ll) {
207: VecDot(x,shell->Y[ll],&yq);
208: VecAXPY(x,shell->beta[ll]-yq*shell->rho[ll],shell->S[ll]);
209:
210: }
211: return(0);
212: }
216:
220: {
221: PetscBool isascii;
223: MatLMVMCtx *lmP;
226: MatShellGetContext(A,(void**)&lmP);
228: PetscTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
229: if (isascii) {
230: PetscViewerASCIIPrintf(pv,"LMVM Matrix\n");
231: PetscViewerASCIIPrintf(pv," Number of vectors: %D\n",lmP->lm);
232: PetscViewerASCIIPrintf(pv," scale type: %s\n",Scale_Table[lmP->scaleType]);
233: PetscViewerASCIIPrintf(pv," rescale type: %s\n",Rescale_Table[lmP->rScaleType]);
234: PetscViewerASCIIPrintf(pv," limit type: %s\n",Limit_Table[lmP->limitType]);
235: PetscViewerASCIIPrintf(pv," updates: %D\n",lmP->nupdates);
236: PetscViewerASCIIPrintf(pv," rejects: %D\n",lmP->nrejects);
237:
238: }
239: else {
240: SETERRQ(PETSC_COMM_SELF,1,"non-ascii viewers not implemented for MatLMVM\n");
241: }
243: return(0);
244:
245: }
250: {
251: MatLMVMCtx *ctx;
255: MatShellGetContext(M,(void**)&ctx);
256: if (ctx->allocated) {
257: if (ctx->Xprev) {
258: PetscObjectDereference((PetscObject)ctx->Xprev);
259: }
260: if (ctx->Gprev) {
261: PetscObjectDereference((PetscObject)ctx->Gprev);
262: }
264: VecDestroyVecs(ctx->lm+1,&ctx->S);
265: VecDestroyVecs(ctx->lm+1,&ctx->Y);
266: VecDestroy(&ctx->D);
267: VecDestroy(&ctx->U);
268: VecDestroy(&ctx->V);
269: VecDestroy(&ctx->W);
270: VecDestroy(&ctx->P);
271: VecDestroy(&ctx->Q);
272: if (ctx->scale) {
273: VecDestroy(&ctx->scale);
274: }
275: }
276: PetscFree(ctx->rho);
277: PetscFree(ctx->beta);
278: PetscFree(ctx->yy_history);
279: PetscFree(ctx->ys_history);
280: PetscFree(ctx->ss_history);
281: PetscFree(ctx->yy_rhistory);
282: PetscFree(ctx->ys_rhistory);
283: PetscFree(ctx->ss_rhistory);
285: PetscFree(ctx);
287: return(0);
288:
289: }
295: {
297: MatLMVMCtx *ctx;
298: PetscInt i;
300: MatShellGetContext(M,(void**)&ctx);
301: if (ctx->Gprev) {
302: PetscObjectDereference((PetscObject)ctx->Gprev);
303: }
304: if (ctx->Xprev) {
305: PetscObjectDereference((PetscObject)ctx->Xprev);
306: }
307: ctx->Gprev = ctx->Y[ctx->lm];
308: ctx->Xprev = ctx->S[ctx->lm];
309: PetscObjectReference((PetscObject)ctx->Gprev);
310: PetscObjectReference((PetscObject)ctx->Xprev);
311: for (i=0; i<ctx->lm; ++i) {
312: ctx->rho[i] = 0.0;
313: }
314: ctx->rho[0] = 1.0;
315:
316: /* Set the scaling and diagonal scaling matrix */
317: switch(ctx->scaleType) {
318: case MatLMVM_Scale_None:
319: ctx->sigma = 1.0;
320: break;
321: case MatLMVM_Scale_Scalar:
322: ctx->sigma = ctx->delta;
323: break;
324: case MatLMVM_Scale_Broyden:
325: VecSet(ctx->D,ctx->delta);
326: break;
327: }
329: ctx->iter=0;
330: ctx->nupdates=0;
331: ctx->lmnow=0;
333:
334: return(0);
335: }
341: {
342: MatLMVMCtx *ctx;
343: PetscReal rhotemp, rhotol;
344: PetscReal y0temp, s0temp;
345: PetscReal yDy, yDs, sDs;
346: PetscReal sigmanew, denom;
348: PetscInt i;
349: PetscBool same;
350: PetscReal yy_sum=0.0, ys_sum=0.0, ss_sum=0.0;
356: PetscTypeCompare((PetscObject)M,MATSHELL,&same);
357: if (!same) {SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");}
358: MatShellGetContext(M,(void**)&ctx);
359: if (!ctx->allocated) {
360: MatLMVMAllocateVectors(M, x);
361: }
363: if (0 == ctx->iter) {
364: MatLMVMReset(M);
365: }
366: else {
367: VecAYPX(ctx->Gprev,-1.0,g);
368: VecAYPX(ctx->Xprev,-1.0,x);
370: VecDot(ctx->Gprev,ctx->Xprev,&rhotemp);
371: VecDot(ctx->Gprev,ctx->Gprev,&y0temp);
373: rhotol = ctx->eps * y0temp;
374: if (rhotemp > rhotol) {
375: ++ctx->nupdates;
377: ctx->lmnow = PetscMin(ctx->lmnow+1, ctx->lm);
378: ierr=PetscObjectDereference((PetscObject)ctx->S[ctx->lm]);
379: ierr=PetscObjectDereference((PetscObject)ctx->Y[ctx->lm]);
380: for (i = ctx->lm-1; i >= 0; --i) {
381: ctx->S[i+1] = ctx->S[i];
382: ctx->Y[i+1] = ctx->Y[i];
383: ctx->rho[i+1] = ctx->rho[i];
384: }
385: ctx->S[0] = ctx->Xprev;
386: ctx->Y[0] = ctx->Gprev;
387: PetscObjectReference((PetscObject)ctx->S[0]);
388: PetscObjectReference((PetscObject)ctx->Y[0]);
389: ctx->rho[0] = 1.0 / rhotemp;
391: /* Compute the scaling */
392: switch(ctx->scaleType) {
393: case MatLMVM_Scale_None:
394: break;
396: case MatLMVM_Scale_Scalar:
397: /* Compute s^T s */
398: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
400: /* Scalar is positive; safeguards are not required. */
402: /* Save information for scalar scaling */
403: ctx->yy_history[(ctx->nupdates - 1) % ctx->scalar_history] = y0temp;
404: ctx->ys_history[(ctx->nupdates - 1) % ctx->scalar_history] = rhotemp;
405: ctx->ss_history[(ctx->nupdates - 1) % ctx->scalar_history] = s0temp;
407: /* Compute summations for scalar scaling */
408: yy_sum = 0; /* No safeguard required; y^T y > 0 */
409: ys_sum = 0; /* No safeguard required; y^T s > 0 */
410: ss_sum = 0; /* No safeguard required; s^T s > 0 */
411: for (i = 0; i < PetscMin(ctx->nupdates, ctx->scalar_history); ++i) {
412: yy_sum += ctx->yy_history[i];
413: ys_sum += ctx->ys_history[i];
414: ss_sum += ctx->ss_history[i];
415: }
417: if (0.0 == ctx->s_alpha) {
418: /* Safeguard ys_sum */
419: if (0.0 == ys_sum) {
420: ys_sum = TAO_ZERO_SAFEGUARD;
421: }
423: sigmanew = ss_sum / ys_sum;
424: }
425: else if (1.0 == ctx->s_alpha) {
426: /* Safeguard yy_sum */
427: if (0.0 == yy_sum) {
428: yy_sum = TAO_ZERO_SAFEGUARD;
429: }
431: sigmanew = ys_sum / yy_sum;
432: }
433: else {
434: denom = 2*ctx->s_alpha*yy_sum;
436: /* Safeguard denom */
437: if (0.0 == denom) {
438: denom = TAO_ZERO_SAFEGUARD;
439: }
441: sigmanew = ((2*ctx->s_alpha-1)*ys_sum +
442: PetscSqrtScalar((2*ctx->s_alpha-1)*(2*ctx->s_alpha-1)*ys_sum*ys_sum -
443: 4*(ctx->s_alpha)*(ctx->s_alpha-1)*yy_sum*ss_sum)) / denom;
444: }
446: switch(ctx->limitType) {
447: case MatLMVM_Limit_Average:
448: if (1.0 == ctx->mu) {
449: ctx->sigma = sigmanew;
450: }
451: else if (ctx->mu) {
452: ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
453: }
454: break;
456: case MatLMVM_Limit_Relative:
457: if (ctx->mu) {
458: ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
459: }
460: break;
462: case MatLMVM_Limit_Absolute:
463: if (ctx->nu) {
464: ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
465: }
466: break;
468: default:
469: ctx->sigma = sigmanew;
470: break;
471: }
472: break;
474: case MatLMVM_Scale_Broyden:
475: /* Original version */
476: /* Combine DFP and BFGS */
478: /* This code appears to be numerically unstable. We use the */
479: /* original version because this was used to generate all of */
480: /* the data and because it may be the least unstable of the */
481: /* bunch. */
483: /* P = Q = inv(D); */
484: VecCopy(ctx->D,ctx->P);
485: VecReciprocal(ctx->P);
486: VecCopy(ctx->P,ctx->Q);
488: /* V = y*y */
489: VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev);
491: /* W = inv(D)*s */
492: VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P);
493: VecDot(ctx->W,ctx->Xprev,&sDs);
495: /* Safeguard rhotemp and sDs */
496: if (0.0 == rhotemp) {
497: rhotemp = TAO_ZERO_SAFEGUARD;
498: }
500: if (0.0 == sDs) {
501: sDs = TAO_ZERO_SAFEGUARD;
502: }
504: if (1.0 != ctx->phi) {
505: /* BFGS portion of the update */
506: /* U = (inv(D)*s)*(inv(D)*s) */
507: VecPointwiseMult(ctx->U,ctx->W,ctx->W);
510: /* Assemble */
511: VecAXPY(ctx->P,1.0/rhotemp,ctx->V);
512: VecAXPY(ctx->P,-1.0/sDs,ctx->U);
513: }
515: if (0.0 != ctx->phi) {
516: /* DFP portion of the update */
517: /* U = inv(D)*s*y */
518: VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev);
520: /* Assemble */
521: VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V);
522: VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U);
523: }
525: if (0.0 == ctx->phi) {
526: VecCopy(ctx->P,ctx->U);
527: }
528: else if (1.0 == ctx->phi) {
529: VecCopy(ctx->Q,ctx->U);
530: }
531: else {
532: /* Broyden update U=(1-phi)*P + phi*Q */
533: VecCopy(ctx->Q,ctx->U);
534: VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P);
535: }
537: /* Obtain inverse and ensure positive definite */
538: VecReciprocal(ctx->U);
539: VecAbs(ctx->U);
542: switch(ctx->rScaleType) {
543: case MatLMVM_Rescale_None:
544: break;
546: case MatLMVM_Rescale_Scalar:
547: case MatLMVM_Rescale_GL:
548: if (ctx->rScaleType == MatLMVM_Rescale_GL) {
549: /* Gilbert and Lemarachal use the old diagonal */
550: VecCopy(ctx->D,ctx->P);
551: }
552: else {
553: /* The default version uses the current diagonal */
554: VecCopy(ctx->U,ctx->P);
555: }
557: /* Compute s^T s */
558: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
560: /* Save information for special cases of scalar rescaling */
561: ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
562: ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
563: ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;
565: if (0.5 == ctx->r_beta) {
566: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
567: VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P);
568: VecDot(ctx->V,ctx->Y[0],&yy_sum);
569:
570: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
571: VecDot(ctx->W,ctx->S[0],&ss_sum);
573: ys_sum = ctx->ys_rhistory[0];
574: }
575: else {
576: VecCopy(ctx->P,ctx->Q);
577: VecReciprocal(ctx->Q);
579: /* Compute summations for scalar scaling */
580: yy_sum = 0; /* No safeguard required */
581: ys_sum = 0; /* No safeguard required */
582: ss_sum = 0; /* No safeguard required */
583: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
584: VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P);
585: VecDot(ctx->V,ctx->Y[i],&yDy);
586: yy_sum += yDy;
588: VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q);
589: VecDot(ctx->W,ctx->S[i],&sDs);
590: ss_sum += sDs;
591: ys_sum += ctx->ys_rhistory[i];
592: }
593: }
594: }
595: else if (0.0 == ctx->r_beta) {
596: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
597: /* Compute summations for scalar scaling */
598: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
600: VecDot(ctx->W, ctx->Y[0], &ys_sum);
601: VecDot(ctx->W, ctx->W, &ss_sum);
602: yy_sum += ctx->yy_rhistory[0];
603: }
604: else {
605: VecCopy(ctx->Q, ctx->P);
606: VecReciprocal(ctx->Q);
608: /* Compute summations for scalar scaling */
609: yy_sum = 0; /* No safeguard required */
610: ys_sum = 0; /* No safeguard required */
611: ss_sum = 0; /* No safeguard required */
612: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
613: VecPointwiseMult(ctx->W, ctx->S[i], ctx->Q);
614: VecDot(ctx->W, ctx->Y[i], &yDs);
615: ys_sum += yDs;
617: VecDot(ctx->W, ctx->W, &sDs);
618: ss_sum += sDs;
619:
620: yy_sum += ctx->yy_rhistory[i];
621: }
622: }
623: }
624: else if (1.0 == ctx->r_beta) {
625: /* Compute summations for scalar scaling */
626: yy_sum = 0; /* No safeguard required */
627: ys_sum = 0; /* No safeguard required */
628: ss_sum = 0; /* No safeguard required */
629: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
630: VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P);
631: VecDot(ctx->V, ctx->S[i], &yDs);
632: ys_sum += yDs;
634: VecDot(ctx->V, ctx->V, &yDy);
635: yy_sum += yDy;
637: ss_sum += ctx->ss_rhistory[i];
638: }
639: }
640: else {
641: VecCopy(ctx->Q, ctx->P);
643: VecPow(ctx->P, ctx->r_beta);
644: VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q);
646: /* Compute summations for scalar scaling */
647: yy_sum = 0; /* No safeguard required */
648: ys_sum = 0; /* No safeguard required */
649: ss_sum = 0; /* No safeguard required */
650: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
651: VecPointwiseMult(ctx->V, ctx->P, ctx->Y[i]);
652: VecPointwiseMult(ctx->W, ctx->Q, ctx->S[i]);
654: VecDot(ctx->V, ctx->V, &yDy);
655: VecDot(ctx->V, ctx->W, &yDs);
656: VecDot(ctx->W, ctx->W, &sDs);
658: yy_sum += yDy;
659: ys_sum += yDs;
660: ss_sum += sDs;
661: }
662: }
664: if (0.0 == ctx->r_alpha) {
665: /* Safeguard ys_sum */
666: if (0.0 == ys_sum) {
667: ys_sum = TAO_ZERO_SAFEGUARD;
668: }
670: sigmanew = ss_sum / ys_sum;
671: }
672: else if (1.0 == ctx->r_alpha) {
673: /* Safeguard yy_sum */
674: if (0.0 == yy_sum) {
675: ys_sum = TAO_ZERO_SAFEGUARD;
676: }
678: sigmanew = ys_sum / yy_sum;
679: }
680: else {
681: denom = 2*ctx->r_alpha*yy_sum;
683: /* Safeguard denom */
684: if (0.0 == denom) {
685: denom = TAO_ZERO_SAFEGUARD;
686: }
688: sigmanew = ((2*ctx->r_alpha-1)*ys_sum +
689: PetscSqrtScalar((2*ctx->r_alpha-1)*(2*ctx->r_alpha-1)*ys_sum*ys_sum -
690: 4*ctx->r_alpha*(ctx->r_alpha-1)*yy_sum*ss_sum)) / denom;
691: }
693: /* If Q has small values, then Q^(r_beta - 1) */
694: /* can have very large values. Hence, ys_sum */
695: /* and ss_sum can be infinity. In this case, */
696: /* sigmanew can either be not-a-number or infinity. */
698: if (PetscIsInfOrNanReal(sigmanew)) {
699: /* sigmanew is not-a-number; skip rescaling */
700: }
701: else if (!sigmanew) {
702: /* sigmanew is zero; this is a bad case; skip rescaling */
703: }
704: else {
705: /* sigmanew is positive */
706: VecScale(ctx->U, sigmanew);
707: }
708: break;
709: }
711: /* Modify for previous information */
712: switch(ctx->limitType) {
713: case MatLMVM_Limit_Average:
714: if (1.0 == ctx->mu) {
715: VecCopy(ctx->D, ctx->U);
716: }
717: else if (ctx->mu) {
718: VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U);
719: }
720: break;
721:
722: case MatLMVM_Limit_Relative:
723: if (ctx->mu) {
724: /* P = (1-mu) * D */
725: VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D);
726: /* Q = (1+mu) * D */
727: VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D);
728: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D);
729: }
730: break;
732: case MatLMVM_Limit_Absolute:
733: if (ctx->nu) {
734: VecCopy(ctx->P, ctx->D);
735: VecShift(ctx->P, -ctx->nu);
736: VecCopy(ctx->D, ctx->Q);
737: VecShift(ctx->Q, ctx->nu);
738: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P);
739: }
740: break;
742: default:
743: VecCopy(ctx->U, ctx->D);
744: break;
745: }
746: break;
747: }
748: PetscObjectDereference((PetscObject)ctx->Xprev);
749: PetscObjectDereference((PetscObject)ctx->Gprev);
750: ctx->Xprev = ctx->S[ctx->lm];
751: ctx->Gprev = ctx->Y[ctx->lm];
752: PetscObjectReference((PetscObject)ctx->S[ctx->lm]);
753: PetscObjectReference((PetscObject)ctx->Y[ctx->lm]);
755: }
756: else {
757: ++ctx->nrejects;
758: }
759: }
760:
761: ++ctx->iter;
762: VecCopy(x, ctx->Xprev);
763: VecCopy(g, ctx->Gprev);
764: return(0);
765: }
770: {
771: MatLMVMCtx *ctx;
773: PetscBool same;
776: PetscTypeCompare((PetscObject)m,MATSHELL,&same);
777: if (!same) {
778: SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
779: }
780: MatShellGetContext(m,(void**)&ctx);
781: ctx->delta = PetscAbsReal(d);
782: ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
783: ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
784: return(0);
785:
786: }
791: {
792: MatLMVMCtx *ctx;
794: PetscBool same;
797: PetscTypeCompare((PetscObject)m,MATSHELL,&same);
798: if (!same) {
799: SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
800: }
801: MatShellGetContext(m,(void**)&ctx);
802:
803: if (ctx->scale) {
804: VecDestroy(&ctx->scale);
805: }
806: if (s) {
807: VecDuplicate(s,&ctx->scale);
808: } else {
809: ctx->scale = PETSC_NULL;
810: }
811: return(0);
812: }
817: {
818: MatLMVMCtx *ctx;
820: PetscBool same;
823: PetscTypeCompare((PetscObject)m,MATSHELL,&same);
824: if (!same) {
825: SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
826: }
827: MatShellGetContext(m,(void**)&ctx);
828: *nrejects = ctx->nrejects;
829: return(0);
830: }
835: {
837: return(0);
838: }
843: {
845: return(0);
846: }
851: {
852: MatLMVMCtx *ctx;
854: PetscBool same;
860: PetscTypeCompare((PetscObject)M,MATSHELL,&same);
861: if (!same) {SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");}
862: MatShellGetContext(M,(void**)&ctx);
863: if (ctx->nupdates == 0) {
864: MatLMVMUpdate(M,x,g);
865: } else {
866: VecCopy(x,ctx->Xprev);
867: VecCopy(g,ctx->Gprev);
868: /* TODO scaling specific terms */
869: }
870: return(0);
871:
872: }
876: {
878: PetscBool same;
884: PetscTypeCompare((PetscObject)coarse,MATSHELL,&same);
885: if (!same) {
886: SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
887: }
888: PetscTypeCompare((PetscObject)op,MATSHELL,&same);
889: if (!same) {
890: SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
891: }
893: return(0);
894: }
899: {
901: MatLMVMCtx *ctx;
902: PetscBool same;
907: PetscTypeCompare((PetscObject)m,MATSHELL,&same);
908: if (!same) {
909: SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
910: }
911: MatShellGetContext(m,(void**)&ctx);
912:
914: /* Perform allocations */
915: VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
916: VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
917: VecDuplicate(v,&ctx->D);
918: VecDuplicate(v,&ctx->U);
919: VecDuplicate(v,&ctx->V);
920: VecDuplicate(v,&ctx->W);
921: VecDuplicate(v,&ctx->P);
922: VecDuplicate(v,&ctx->Q);
923: ctx->allocated = PETSC_TRUE;
925: return(0);
926: }