Actual source code: lmvmmat.c
petsc-3.9.4 2018-09-11
1: #include <../src/tao/matrix/lmvmmat.h>
2: #include <petsctao.h>
3: #include <petscksp.h>
4: #include <petsc/private/petscimpl.h>
6: /* This is a vile hack */
7: #if defined(PETSC_USE_COMPLEX)
8: #define VecDot VecDotRealPart
9: #endif
11: #define TaoMid(a,b,c) (((a) < (b)) ? \
12: (((b) < (c)) ? (b) : \
13: (((a) < (c)) ? (c) : (a))) : \
14: (((a) < (c)) ? (a) : \
15: (((b) < (c)) ? (c) : (b))))
17: /* These lists are used for setting options */
18: static const char *Scale_Table[64] = {"none","scalar","broyden"};
20: static const char *Rescale_Table[64] = {"none","scalar","gl"};
22: static const char *Limit_Table[64] = {"none","average","relative","absolute"};
24: /*@C
25: MatCreateLMVM - Creates a limited memory matrix for lmvm algorithms.
27: Collective on A
29: Input Parameters:
30: + comm - MPI Communicator
31: . n - local size of vectors
32: - N - global size of vectors
34: Output Parameters:
35: . A - New LMVM matrix
37: Level: developer
39: @*/
40: PetscErrorCode MatCreateLMVM(MPI_Comm comm, PetscInt n, PetscInt N, Mat *A)
41: {
42: MatLMVMCtx *ctx;
44: PetscInt nhistory;
47: /* create data structure and populate with default values */
48: PetscNew(&ctx);
49: ctx->lm=5;
50: ctx->eps=0.0;
51: ctx->limitType=MatLMVM_Limit_None;
52: ctx->scaleType=MatLMVM_Scale_Broyden;
53: ctx->rScaleType = MatLMVM_Rescale_Scalar;
54: ctx->s_alpha = 1.0;
55: ctx->r_alpha = 1.0;
56: ctx->r_beta = 0.5;
57: ctx->mu = 1.0;
58: ctx->nu = 100.0;
60: ctx->phi = 0.125;
62: ctx->scalar_history = 1;
63: ctx->rescale_history = 1;
65: ctx->delta_min = 1e-7;
66: ctx->delta_max = 100.0;
67:
68: ctx->recycle = PETSC_FALSE;
70: /* Begin configuration */
71: PetscOptionsBegin(comm,NULL,NULL,NULL);
72: PetscOptionsInt("-tao_lmm_vectors", "vectors to use for approximation", "", ctx->lm, &ctx->lm,NULL);
73: PetscOptionsReal("-tao_lmm_limit_mu", "mu limiting factor", "", ctx->mu, &ctx->mu,NULL);
74: PetscOptionsReal("-tao_lmm_limit_nu", "nu limiting factor", "", ctx->nu, &ctx->nu,NULL);
75: PetscOptionsReal("-tao_lmm_broyden_phi", "phi factor for Broyden scaling", "", ctx->phi, &ctx->phi,NULL);
76: PetscOptionsReal("-tao_lmm_scalar_alpha", "alpha factor for scalar scaling", "",ctx->s_alpha, &ctx->s_alpha,NULL);
77: PetscOptionsReal("-tao_lmm_rescale_alpha", "alpha factor for rescaling diagonal", "", ctx->r_alpha, &ctx->r_alpha,NULL);
78: PetscOptionsReal("-tao_lmm_rescale_beta", "beta factor for rescaling diagonal", "", ctx->r_beta, &ctx->r_beta,NULL);
79: PetscOptionsInt("-tao_lmm_scalar_history", "amount of history for scalar scaling", "", ctx->scalar_history, &ctx->scalar_history,NULL);
80: PetscOptionsInt("-tao_lmm_rescale_history", "amount of history for rescaling diagonal", "", ctx->rescale_history, &ctx->rescale_history,NULL);
81: PetscOptionsReal("-tao_lmm_eps", "rejection tolerance", "", ctx->eps, &ctx->eps,NULL);
82: PetscOptionsEList("-tao_lmm_scale_type", "scale type", "", Scale_Table, MatLMVM_Scale_Types, Scale_Table[ctx->scaleType], &ctx->scaleType,NULL);
83: PetscOptionsEList("-tao_lmm_rescale_type", "rescale type", "", Rescale_Table, MatLMVM_Rescale_Types, Rescale_Table[ctx->rScaleType], &ctx->rScaleType,NULL);
84: PetscOptionsEList("-tao_lmm_limit_type", "limit type", "", Limit_Table, MatLMVM_Limit_Types, Limit_Table[ctx->limitType], &ctx->limitType,NULL);
85: PetscOptionsReal("-tao_lmm_delta_min", "minimum delta value", "", ctx->delta_min, &ctx->delta_min,NULL);
86: PetscOptionsReal("-tao_lmm_delta_max", "maximum delta value", "", ctx->delta_max, &ctx->delta_max,NULL);
87: PetscOptionsBool("-tao_lmm_recycle", "enable recycling of accumulated corrections between subsequent TaoSolve() calls", "", ctx->recycle, &ctx->recycle, NULL);
88: PetscOptionsEnd();
90: /* Complete configuration */
91: ctx->rescale_history = PetscMin(ctx->rescale_history, ctx->lm);
92: PetscMalloc1(ctx->lm+1,&ctx->rho);
93: PetscMalloc1(ctx->lm+1,&ctx->beta);
95: nhistory = PetscMax(ctx->scalar_history,1);
96: PetscMalloc1(nhistory,&ctx->yy_history);
97: PetscMalloc1(nhistory,&ctx->ys_history);
98: PetscMalloc1(nhistory,&ctx->ss_history);
100: nhistory = PetscMax(ctx->rescale_history,1);
101: PetscMalloc1(nhistory,&ctx->yy_rhistory);
102: PetscMalloc1(nhistory,&ctx->ys_rhistory);
103: PetscMalloc1(nhistory,&ctx->ss_rhistory);
105: /* Finish initializations */
106: ctx->lmnow = 0;
107: ctx->iter = 0;
108: ctx->nupdates = 0;
109: ctx->nrejects = 0;
110: ctx->delta = 1.0;
112: ctx->Gprev = 0;
113: ctx->Xprev = 0;
115: ctx->scale = 0;
116: ctx->useScale = PETSC_FALSE;
118: ctx->H0_mat = 0;
119: ctx->H0_ksp = 0;
120: ctx->H0_norm = 0;
121: ctx->useDefaultH0 = PETSC_TRUE;
123: MatCreateShell(comm, n, n, N, N, ctx, A);
124: MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
125: MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
126: return(0);
127: }
129: PetscErrorCode MatLMVMSolve(Mat A, Vec b, Vec x)
130: {
131: PetscReal sq, yq, dd;
132: PetscInt ll;
133: PetscBool scaled;
134: MatLMVMCtx *shell;
141: MatShellGetContext(A,(void**)&shell);
142: if (shell->lmnow < 1) {
143: shell->rho[0] = 1.0;
144: shell->theta = 1.0;
145: }
147: VecCopy(b,x);
148: for (ll = 0; ll < shell->lmnow; ++ll) {
149: VecDot(x,shell->S[ll],&sq);
150: shell->beta[ll] = sq * shell->rho[ll];
151: VecAXPY(x,-shell->beta[ll],shell->Y[ll]);
152: }
154: scaled = PETSC_FALSE;
155: if (!scaled && !shell->useDefaultH0 && shell->H0_mat) {
156: KSPSolve(shell->H0_ksp,x,shell->U);
157: VecScale(shell->U, shell->theta);
158: VecDot(x,shell->U,&dd);
159: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
160: /* Accept Hessian solve */
161: VecCopy(shell->U,x);
162: scaled = PETSC_TRUE;
163: }
164: }
166: if (!scaled && shell->useScale) {
167: VecPointwiseMult(shell->U,x,shell->scale);
168: VecDot(x,shell->U,&dd);
169: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
170: /* Accept scaling */
171: VecCopy(shell->U,x);
172: scaled = PETSC_TRUE;
173: }
174: }
176: if (!scaled) {
177: switch(shell->scaleType) {
178: case MatLMVM_Scale_None:
179: break;
181: case MatLMVM_Scale_Scalar:
182: VecScale(x,shell->sigma);
183: break;
185: case MatLMVM_Scale_Broyden:
186: VecPointwiseMult(x,x,shell->D);
187: break;
188: }
189: }
190: for (ll = shell->lmnow-1; ll >= 0; --ll) {
191: VecDot(x,shell->Y[ll],&yq);
192: VecAXPY(x,shell->beta[ll]-yq*shell->rho[ll],shell->S[ll]);
193: }
194: return(0);
195: }
197: PetscErrorCode MatView_LMVM(Mat A, PetscViewer pv)
198: {
199: PetscBool isascii;
201: MatLMVMCtx *lmP;
204: MatShellGetContext(A,(void**)&lmP);
205: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
206: if (isascii) {
207: PetscViewerASCIIPrintf(pv,"LMVM Matrix\n");
208: PetscViewerASCIIPrintf(pv," Number of vectors: %D\n",lmP->lm);
209: PetscViewerASCIIPrintf(pv," scale type: %s\n",Scale_Table[lmP->scaleType]);
210: PetscViewerASCIIPrintf(pv," rescale type: %s\n",Rescale_Table[lmP->rScaleType]);
211: PetscViewerASCIIPrintf(pv," limit type: %s\n",Limit_Table[lmP->limitType]);
212: PetscViewerASCIIPrintf(pv," updates: %D\n",lmP->nupdates);
213: PetscViewerASCIIPrintf(pv," rejects: %D\n",lmP->nrejects);
214: }
215: return(0);
216: }
218: PetscErrorCode MatDestroy_LMVM(Mat M)
219: {
220: MatLMVMCtx *ctx;
224: MatShellGetContext(M,(void**)&ctx);
225: if (ctx->allocated) {
226: if (ctx->Xprev) {
227: PetscObjectDereference((PetscObject)ctx->Xprev);
228: }
229: if (ctx->Gprev) {
230: PetscObjectDereference((PetscObject)ctx->Gprev);
231: }
233: VecDestroyVecs(ctx->lm+1,&ctx->S);
234: VecDestroyVecs(ctx->lm+1,&ctx->Y);
235: VecDestroy(&ctx->D);
236: VecDestroy(&ctx->U);
237: VecDestroy(&ctx->V);
238: VecDestroy(&ctx->W);
239: VecDestroy(&ctx->P);
240: VecDestroy(&ctx->Q);
241: VecDestroy(&ctx->H0_norm);
242: if (ctx->scale) {
243: VecDestroy(&ctx->scale);
244: }
245: }
246: if (ctx->H0_mat) {
247: PetscObjectDereference((PetscObject)ctx->H0_mat);
248: KSPDestroy(&ctx->H0_ksp);
249: }
250: PetscFree(ctx->rho);
251: PetscFree(ctx->beta);
252: PetscFree(ctx->yy_history);
253: PetscFree(ctx->ys_history);
254: PetscFree(ctx->ss_history);
255: PetscFree(ctx->yy_rhistory);
256: PetscFree(ctx->ys_rhistory);
257: PetscFree(ctx->ss_rhistory);
258: PetscFree(ctx);
259: return(0);
260: }
262: PetscErrorCode MatLMVMGetUpdates(Mat M, PetscInt *nupdates)
263: {
265: MatLMVMCtx *ctx;
266:
268: MatShellGetContext(M, (void**)&ctx);
269: *nupdates = ctx->nupdates;
270: return(0);
271: }
273: PetscErrorCode MatLMVMSetRecycleFlag(Mat M, PetscBool flg)
274: {
276: MatLMVMCtx *ctx;
277:
279: MatShellGetContext(M, (void**)&ctx);
280: ctx->recycle = flg;
281: return(0);
282: }
284: PetscErrorCode MatLMVMGetRecycleFlag(Mat M, PetscBool *flg)
285: {
287: MatLMVMCtx *ctx;
288:
290: MatShellGetContext(M, (void**)&ctx);
291: *flg = ctx->recycle;
292: return(0);
293: }
295: PetscErrorCode MatLMVMReset(Mat M)
296: {
298: MatLMVMCtx *ctx;
299: PetscInt i;
302: MatShellGetContext(M,(void**)&ctx);
303: if (ctx->recycle) {
304: PetscInfo(M, "WARNING: Correction vectors are being reset while recycling is enabled!\n");
305: }
306: if (ctx->Gprev) {
307: PetscObjectDereference((PetscObject)ctx->Gprev);
308: }
309: if (ctx->Xprev) {
310: PetscObjectDereference((PetscObject)ctx->Xprev);
311: }
312: ctx->Gprev = ctx->Y[ctx->lm];
313: ctx->Xprev = ctx->S[ctx->lm];
314: PetscObjectReference((PetscObject)ctx->Gprev);
315: PetscObjectReference((PetscObject)ctx->Xprev);
316: for (i=0; i<ctx->lm; ++i) {
317: ctx->rho[i] = 0.0;
318: }
319: ctx->rho[0] = 1.0;
321: /* Set the scaling and diagonal scaling matrix */
322: switch(ctx->scaleType) {
323: case MatLMVM_Scale_None:
324: ctx->sigma = 1.0;
325: break;
326: case MatLMVM_Scale_Scalar:
327: ctx->sigma = ctx->delta;
328: break;
329: case MatLMVM_Scale_Broyden:
330: VecSet(ctx->D,ctx->delta);
331: break;
332: }
334: ctx->iter=0;
335: ctx->nupdates=0;
336: ctx->lmnow=0;
337: return(0);
338: }
340: PetscErrorCode MatLMVMUpdate(Mat M, Vec x, Vec g)
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;
355: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
356: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
357: MatShellGetContext(M,(void**)&ctx);
358: if (!ctx->allocated) {
359: MatLMVMAllocateVectors(M, x);
360: }
362: if (0 == ctx->iter) {
363: MatLMVMReset(M);
364: } else {
365: VecAYPX(ctx->Gprev,-1.0,g);
366: VecAYPX(ctx->Xprev,-1.0,x);
368: VecDot(ctx->Gprev,ctx->Xprev,&rhotemp);
369: VecDot(ctx->Gprev,ctx->Gprev,&y0temp);
371: rhotol = ctx->eps * y0temp;
372: if (rhotemp > rhotol) {
373: ++ctx->nupdates;
375: ctx->lmnow = PetscMin(ctx->lmnow+1, ctx->lm);
376: ierr=PetscObjectDereference((PetscObject)ctx->S[ctx->lm]);
377: ierr=PetscObjectDereference((PetscObject)ctx->Y[ctx->lm]);
378: for (i = ctx->lm-1; i >= 0; --i) {
379: ctx->S[i+1] = ctx->S[i];
380: ctx->Y[i+1] = ctx->Y[i];
381: ctx->rho[i+1] = ctx->rho[i];
382: }
383: ctx->S[0] = ctx->Xprev;
384: ctx->Y[0] = ctx->Gprev;
385: PetscObjectReference((PetscObject)ctx->S[0]);
386: PetscObjectReference((PetscObject)ctx->Y[0]);
387: ctx->rho[0] = 1.0 / rhotemp;
389: /* Compute the scaling */
390: switch(ctx->scaleType) {
391: case MatLMVM_Scale_None:
392: break;
394: case MatLMVM_Scale_Scalar:
395: /* Compute s^T s */
396: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
398: /* Scalar is positive; safeguards are not required. */
400: /* Save information for scalar scaling */
401: ctx->yy_history[(ctx->nupdates - 1) % ctx->scalar_history] = y0temp;
402: ctx->ys_history[(ctx->nupdates - 1) % ctx->scalar_history] = rhotemp;
403: ctx->ss_history[(ctx->nupdates - 1) % ctx->scalar_history] = s0temp;
405: /* Compute summations for scalar scaling */
406: yy_sum = 0; /* No safeguard required; y^T y > 0 */
407: ys_sum = 0; /* No safeguard required; y^T s > 0 */
408: ss_sum = 0; /* No safeguard required; s^T s > 0 */
409: for (i = 0; i < PetscMin(ctx->nupdates, ctx->scalar_history); ++i) {
410: yy_sum += ctx->yy_history[i];
411: ys_sum += ctx->ys_history[i];
412: ss_sum += ctx->ss_history[i];
413: }
415: if (0.0 == ctx->s_alpha) {
416: /* Safeguard ys_sum */
417: if (0.0 == ys_sum) {
418: ys_sum = TAO_ZERO_SAFEGUARD;
419: }
421: sigmanew = ss_sum / ys_sum;
422: } else if (1.0 == ctx->s_alpha) {
423: /* Safeguard yy_sum */
424: if (0.0 == yy_sum) {
425: yy_sum = TAO_ZERO_SAFEGUARD;
426: }
428: sigmanew = ys_sum / yy_sum;
429: } else {
430: denom = 2*ctx->s_alpha*yy_sum;
432: /* Safeguard denom */
433: if (0.0 == denom) {
434: denom = TAO_ZERO_SAFEGUARD;
435: }
437: sigmanew = ((2*ctx->s_alpha-1)*ys_sum + PetscSqrtReal((2*ctx->s_alpha-1)*(2*ctx->s_alpha-1)*ys_sum*ys_sum - 4*(ctx->s_alpha)*(ctx->s_alpha-1)*yy_sum*ss_sum)) / denom;
438: }
440: switch(ctx->limitType) {
441: case MatLMVM_Limit_Average:
442: if (1.0 == ctx->mu) {
443: ctx->sigma = sigmanew;
444: } else if (ctx->mu) {
445: ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
446: }
447: break;
449: case MatLMVM_Limit_Relative:
450: if (ctx->mu) {
451: ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
452: }
453: break;
455: case MatLMVM_Limit_Absolute:
456: if (ctx->nu) {
457: ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
458: }
459: break;
461: default:
462: ctx->sigma = sigmanew;
463: break;
464: }
465: break;
467: case MatLMVM_Scale_Broyden:
468: /* Original version */
469: /* Combine DFP and BFGS */
471: /* This code appears to be numerically unstable. We use the */
472: /* original version because this was used to generate all of */
473: /* the data and because it may be the least unstable of the */
474: /* bunch. */
476: /* P = Q = inv(D); */
477: VecCopy(ctx->D,ctx->P);
478: VecReciprocal(ctx->P);
479: VecCopy(ctx->P,ctx->Q);
481: /* V = y*y */
482: VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev);
484: /* W = inv(D)*s */
485: VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P);
486: VecDot(ctx->W,ctx->Xprev,&sDs);
488: /* Safeguard rhotemp and sDs */
489: if (0.0 == rhotemp) {
490: rhotemp = TAO_ZERO_SAFEGUARD;
491: }
493: if (0.0 == sDs) {
494: sDs = TAO_ZERO_SAFEGUARD;
495: }
497: if (1.0 != ctx->phi) {
498: /* BFGS portion of the update */
499: /* U = (inv(D)*s)*(inv(D)*s) */
500: VecPointwiseMult(ctx->U,ctx->W,ctx->W);
502: /* Assemble */
503: VecAXPY(ctx->P,1.0/rhotemp,ctx->V);
504: VecAXPY(ctx->P,-1.0/sDs,ctx->U);
505: }
507: if (0.0 != ctx->phi) {
508: /* DFP portion of the update */
509: /* U = inv(D)*s*y */
510: VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev);
512: /* Assemble */
513: VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V);
514: VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U);
515: }
517: if (0.0 == ctx->phi) {
518: VecCopy(ctx->P,ctx->U);
519: } else if (1.0 == ctx->phi) {
520: VecCopy(ctx->Q,ctx->U);
521: } else {
522: /* Broyden update U=(1-phi)*P + phi*Q */
523: VecCopy(ctx->Q,ctx->U);
524: VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P);
525: }
527: /* Obtain inverse and ensure positive definite */
528: VecReciprocal(ctx->U);
529: VecAbs(ctx->U);
531: switch(ctx->rScaleType) {
532: case MatLMVM_Rescale_None:
533: break;
535: case MatLMVM_Rescale_Scalar:
536: case MatLMVM_Rescale_GL:
537: if (ctx->rScaleType == MatLMVM_Rescale_GL) {
538: /* Gilbert and Lemarachal use the old diagonal */
539: VecCopy(ctx->D,ctx->P);
540: } else {
541: /* The default version uses the current diagonal */
542: VecCopy(ctx->U,ctx->P);
543: }
545: /* Compute s^T s */
546: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
548: /* Save information for special cases of scalar rescaling */
549: ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
550: ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
551: ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;
553: if (0.5 == ctx->r_beta) {
554: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
555: VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P);
556: VecDot(ctx->V,ctx->Y[0],&yy_sum);
558: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
559: VecDot(ctx->W,ctx->S[0],&ss_sum);
561: ys_sum = ctx->ys_rhistory[0];
562: } else {
563: VecCopy(ctx->P,ctx->Q);
564: VecReciprocal(ctx->Q);
566: /* Compute summations for scalar scaling */
567: yy_sum = 0; /* No safeguard required */
568: ys_sum = 0; /* No safeguard required */
569: ss_sum = 0; /* No safeguard required */
570: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
571: VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P);
572: VecDot(ctx->V,ctx->Y[i],&yDy);
573: yy_sum += yDy;
575: VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q);
576: VecDot(ctx->W,ctx->S[i],&sDs);
577: ss_sum += sDs;
578: ys_sum += ctx->ys_rhistory[i];
579: }
580: }
581: } else if (0.0 == ctx->r_beta) {
582: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
583: /* Compute summations for scalar scaling */
584: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
586: VecDot(ctx->W, ctx->Y[0], &ys_sum);
587: VecDot(ctx->W, ctx->W, &ss_sum);
588: yy_sum += ctx->yy_rhistory[0];
589: } else {
590: VecCopy(ctx->Q, ctx->P);
591: VecReciprocal(ctx->Q);
593: /* Compute summations for scalar scaling */
594: yy_sum = 0; /* No safeguard required */
595: ys_sum = 0; /* No safeguard required */
596: ss_sum = 0; /* No safeguard required */
597: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
598: VecPointwiseMult(ctx->W, ctx->S[i], ctx->Q);
599: VecDot(ctx->W, ctx->Y[i], &yDs);
600: ys_sum += yDs;
602: VecDot(ctx->W, ctx->W, &sDs);
603: ss_sum += sDs;
605: yy_sum += ctx->yy_rhistory[i];
606: }
607: }
608: } else if (1.0 == ctx->r_beta) {
609: /* Compute summations for scalar scaling */
610: yy_sum = 0; /* No safeguard required */
611: ys_sum = 0; /* No safeguard required */
612: ss_sum = 0; /* No safeguard required */
613: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
614: VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P);
615: VecDot(ctx->V, ctx->S[i], &yDs);
616: ys_sum += yDs;
618: VecDot(ctx->V, ctx->V, &yDy);
619: yy_sum += yDy;
621: ss_sum += ctx->ss_rhistory[i];
622: }
623: } else {
624: VecCopy(ctx->Q, ctx->P);
626: VecPow(ctx->P, ctx->r_beta);
627: VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q);
629: /* Compute summations for scalar scaling */
630: yy_sum = 0; /* No safeguard required */
631: ys_sum = 0; /* No safeguard required */
632: ss_sum = 0; /* No safeguard required */
633: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
634: VecPointwiseMult(ctx->V, ctx->P, ctx->Y[i]);
635: VecPointwiseMult(ctx->W, ctx->Q, ctx->S[i]);
637: VecDot(ctx->V, ctx->V, &yDy);
638: VecDot(ctx->V, ctx->W, &yDs);
639: VecDot(ctx->W, ctx->W, &sDs);
641: yy_sum += yDy;
642: ys_sum += yDs;
643: ss_sum += sDs;
644: }
645: }
647: if (0.0 == ctx->r_alpha) {
648: /* Safeguard ys_sum */
649: if (0.0 == ys_sum) {
650: ys_sum = TAO_ZERO_SAFEGUARD;
651: }
653: sigmanew = ss_sum / ys_sum;
654: } else if (1.0 == ctx->r_alpha) {
655: /* Safeguard yy_sum */
656: if (0.0 == yy_sum) {
657: ys_sum = TAO_ZERO_SAFEGUARD;
658: }
660: sigmanew = ys_sum / yy_sum;
661: } else {
662: denom = 2*ctx->r_alpha*yy_sum;
664: /* Safeguard denom */
665: if (0.0 == denom) {
666: denom = TAO_ZERO_SAFEGUARD;
667: }
669: sigmanew = ((2*ctx->r_alpha-1)*ys_sum + PetscSqrtReal((2*ctx->r_alpha-1)*(2*ctx->r_alpha-1)*ys_sum*ys_sum - 4*ctx->r_alpha*(ctx->r_alpha-1)*yy_sum*ss_sum)) / denom;
670: }
672: /* If Q has small values, then Q^(r_beta - 1) */
673: /* can have very large values. Hence, ys_sum */
674: /* and ss_sum can be infinity. In this case, */
675: /* sigmanew can either be not-a-number or infinity. */
677: if (PetscIsInfOrNanReal(sigmanew)) {
678: /* sigmanew is not-a-number; skip rescaling */
679: } else if (!sigmanew) {
680: /* sigmanew is zero; this is a bad case; skip rescaling */
681: } else {
682: /* sigmanew is positive */
683: VecScale(ctx->U, sigmanew);
684: }
685: break;
686: }
688: /* Modify for previous information */
689: switch(ctx->limitType) {
690: case MatLMVM_Limit_Average:
691: if (1.0 == ctx->mu) {
692: VecCopy(ctx->D, ctx->U);
693: } else if (ctx->mu) {
694: VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U);
695: }
696: break;
698: case MatLMVM_Limit_Relative:
699: if (ctx->mu) {
700: /* P = (1-mu) * D */
701: VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D);
702: /* Q = (1+mu) * D */
703: VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D);
704: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D);
705: }
706: break;
708: case MatLMVM_Limit_Absolute:
709: if (ctx->nu) {
710: VecCopy(ctx->P, ctx->D);
711: VecShift(ctx->P, -ctx->nu);
712: VecCopy(ctx->D, ctx->Q);
713: VecShift(ctx->Q, ctx->nu);
714: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P);
715: }
716: break;
718: default:
719: VecCopy(ctx->U, ctx->D);
720: break;
721: }
722: break;
723: }
724: PetscObjectDereference((PetscObject)ctx->Xprev);
725: PetscObjectDereference((PetscObject)ctx->Gprev);
726: ctx->Xprev = ctx->S[ctx->lm];
727: ctx->Gprev = ctx->Y[ctx->lm];
728: PetscObjectReference((PetscObject)ctx->S[ctx->lm]);
729: PetscObjectReference((PetscObject)ctx->Y[ctx->lm]);
731: } else {
732: ++ctx->nrejects;
733: }
734: }
736: ++ctx->iter;
737: VecCopy(x, ctx->Xprev);
738: VecCopy(g, ctx->Gprev);
739: return(0);
740: }
742: PetscErrorCode MatLMVMSetDelta(Mat m, PetscReal d)
743: {
744: MatLMVMCtx *ctx;
746: PetscBool same;
750: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
751: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
752: MatShellGetContext(m,(void**)&ctx);
753: ctx->delta = PetscAbsReal(d);
754: ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
755: ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
756: return(0);
757: }
759: PetscErrorCode MatLMVMSetScale(Mat m, Vec s)
760: {
761: MatLMVMCtx *ctx;
763: PetscBool same;
767: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
768: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
769: MatShellGetContext(m,(void**)&ctx);
771: VecDestroy(&ctx->scale);
772: if (s) {
773: VecDuplicate(s,&ctx->scale);
774: }
775: return(0);
776: }
778: PetscErrorCode MatLMVMGetRejects(Mat m, PetscInt *nrejects)
779: {
780: MatLMVMCtx *ctx;
782: PetscBool same;
786: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
787: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
788: MatShellGetContext(m,(void**)&ctx);
789: *nrejects = ctx->nrejects;
790: return(0);
791: }
793: PetscErrorCode MatLMVMSetH0(Mat m, Mat H0)
794: {
795: MatLMVMCtx *ctx;
797: PetscBool same;
801: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
802: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
803: MatShellGetContext(m,(void**)&ctx);
805: ctx->H0_mat = H0;
806: PetscObjectReference((PetscObject)ctx->H0_mat);
808: ctx->useDefaultH0 = PETSC_FALSE;
810: KSPCreate(PetscObjectComm((PetscObject)H0), &ctx->H0_ksp);
811: PetscObjectIncrementTabLevel((PetscObject)ctx->H0_ksp, (PetscObject)ctx, 1);
812: KSPSetOperators(ctx->H0_ksp, H0, H0);
813: /* its options prefix and setup is handled in TaoSolve_LMVM/TaoSolve_BLMVM */
814: return(0);
815: }
817: PetscErrorCode MatLMVMGetH0(Mat m, Mat *H0)
818: {
819: MatLMVMCtx *ctx;
821: PetscBool same;
825: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
826: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
828: MatShellGetContext(m,(void**)&ctx);
829: *H0 = ctx->H0_mat;
830: return(0);
831: }
833: PetscErrorCode MatLMVMGetH0KSP(Mat m, KSP *H0ksp)
834: {
835: MatLMVMCtx *ctx;
837: PetscBool same;
841: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
842: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
844: MatShellGetContext(m,(void**)&ctx);
845: *H0ksp = ctx->H0_ksp;
846: return(0);
847: }
849: PetscErrorCode MatLMVMGetX0(Mat m, Vec x)
850: {
852: return(0);
853: }
855: PetscErrorCode MatLMVMSetPrev(Mat M, Vec x, Vec g)
856: {
857: MatLMVMCtx *ctx;
859: PetscBool same;
864: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
865: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
866: MatShellGetContext(M,(void**)&ctx);
867: if (ctx->nupdates == 0) {
868: MatLMVMUpdate(M,x,g);
869: } else {
870: VecCopy(x,ctx->Xprev);
871: VecCopy(g,ctx->Gprev);
872: /* TODO scaling specific terms */
873: }
874: return(0);
875: }
877: PetscErrorCode MatLMVMRefine(Mat coarse, Mat op, Vec fineX, Vec fineG)
878: {
880: PetscBool same;
887: PetscObjectTypeCompare((PetscObject)coarse,MATSHELL,&same);
888: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
889: PetscObjectTypeCompare((PetscObject)op,MATSHELL,&same);
890: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
891: return(0);
892: }
894: PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v)
895: {
897: MatLMVMCtx *ctx;
898: PetscBool same;
903: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
904: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
905: MatShellGetContext(m,(void**)&ctx);
907: /* Perform allocations */
908: VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
909: VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
910: VecDuplicate(v,&ctx->D);
911: VecDuplicate(v,&ctx->U);
912: VecDuplicate(v,&ctx->V);
913: VecDuplicate(v,&ctx->W);
914: VecDuplicate(v,&ctx->P);
915: VecDuplicate(v,&ctx->Q);
916: VecDuplicate(v,&ctx->H0_norm);
917: ctx->allocated = PETSC_TRUE;
918: return(0);
919: }