Actual source code: lmvmmat.c
petsc-3.7.5 2017-01-01
1: #include <../src/tao/matrix/lmvmmat.h> /*I "lmvmmat.h" */
2: #include <petsctao.h> /*I "petsctao.h" */
3: #include <petscksp.h>
5: /* This is a vile hack */
6: #if defined(PETSC_USE_COMPLEX)
7: #define VecDot VecDotRealPart
8: #endif
10: #define TaoMid(a,b,c) (((a) < (b)) ? \
11: (((b) < (c)) ? (b) : \
12: (((a) < (c)) ? (c) : (a))) : \
13: (((a) < (c)) ? (a) : \
14: (((b) < (c)) ? (c) : (b))))
16: /* These lists are used for setting options */
17: static const char *Scale_Table[64] = {"none","scalar","broyden"};
19: static const char *Rescale_Table[64] = {"none","scalar","gl"};
21: static const char *Limit_Table[64] = {"none","average","relative","absolute"};
25: /*@C
26: MatCreateLMVM - Creates a limited memory matrix for lmvm algorithms.
28: Collective on A
30: Input Parameters:
31: + comm - MPI Communicator
32: . n - local size of vectors
33: - N - global size of vectors
35: Output Parameters:
36: . A - New LMVM matrix
38: Level: developer
40: @*/
41: extern PetscErrorCode MatCreateLMVM(MPI_Comm comm, PetscInt n, PetscInt N, Mat *A)
42: {
43: MatLMVMCtx *ctx;
45: PetscInt nhistory;
48: /* create data structure and populate with default values */
49: PetscNew(&ctx);
50: ctx->lm=5;
51: ctx->eps=0.0;
52: ctx->limitType=MatLMVM_Limit_None;
53: ctx->scaleType=MatLMVM_Scale_Broyden;
54: ctx->rScaleType = MatLMVM_Rescale_Scalar;
55: ctx->s_alpha = 1.0;
56: ctx->r_alpha = 1.0;
57: ctx->r_beta = 0.5;
58: ctx->mu = 1.0;
59: ctx->nu = 100.0;
60:
61: ctx->phi = 0.125;
62:
63: ctx->scalar_history = 1;
64: ctx->rescale_history = 1;
65:
66: ctx->delta_min = 1e-7;
67: ctx->delta_max = 100.0;
69: /* Begin configuration */
70: PetscOptionsBegin(comm,NULL,NULL,NULL);
71: PetscOptionsInt("-tao_lmm_vectors", "vectors to use for approximation", "", ctx->lm, &ctx->lm,NULL);
72: PetscOptionsReal("-tao_lmm_limit_mu", "mu limiting factor", "", ctx->mu, &ctx->mu,NULL);
73: PetscOptionsReal("-tao_lmm_limit_nu", "nu limiting factor", "", ctx->nu, &ctx->nu,NULL);
74: PetscOptionsReal("-tao_lmm_broyden_phi", "phi factor for Broyden scaling", "", ctx->phi, &ctx->phi,NULL);
75: PetscOptionsReal("-tao_lmm_scalar_alpha", "alpha factor for scalar scaling", "",ctx->s_alpha, &ctx->s_alpha,NULL);
76: PetscOptionsReal("-tao_lmm_rescale_alpha", "alpha factor for rescaling diagonal", "", ctx->r_alpha, &ctx->r_alpha,NULL);
77: PetscOptionsReal("-tao_lmm_rescale_beta", "beta factor for rescaling diagonal", "", ctx->r_beta, &ctx->r_beta,NULL);
78: PetscOptionsInt("-tao_lmm_scalar_history", "amount of history for scalar scaling", "", ctx->scalar_history, &ctx->scalar_history,NULL);
79: PetscOptionsInt("-tao_lmm_rescale_history", "amount of history for rescaling diagonal", "", ctx->rescale_history, &ctx->rescale_history,NULL);
80: PetscOptionsReal("-tao_lmm_eps", "rejection tolerance", "", ctx->eps, &ctx->eps,NULL);
81: PetscOptionsEList("-tao_lmm_scale_type", "scale type", "", Scale_Table, MatLMVM_Scale_Types, Scale_Table[ctx->scaleType], &ctx->scaleType,NULL);
82: PetscOptionsEList("-tao_lmm_rescale_type", "rescale type", "", Rescale_Table, MatLMVM_Rescale_Types, Rescale_Table[ctx->rScaleType], &ctx->rScaleType,NULL);
83: PetscOptionsEList("-tao_lmm_limit_type", "limit type", "", Limit_Table, MatLMVM_Limit_Types, Limit_Table[ctx->limitType], &ctx->limitType,NULL);
84: PetscOptionsReal("-tao_lmm_delta_min", "minimum delta value", "", ctx->delta_min, &ctx->delta_min,NULL);
85: PetscOptionsReal("-tao_lmm_delta_max", "maximum delta value", "", ctx->delta_max, &ctx->delta_max,NULL);
86: PetscOptionsEnd();
88: /* Complete configuration */
89: ctx->rescale_history = PetscMin(ctx->rescale_history, ctx->lm);
90: PetscMalloc1(ctx->lm+1,&ctx->rho);
91: PetscMalloc1(ctx->lm+1,&ctx->beta);
93: nhistory = PetscMax(ctx->scalar_history,1);
94: PetscMalloc1(nhistory,&ctx->yy_history);
95: PetscMalloc1(nhistory,&ctx->ys_history);
96: PetscMalloc1(nhistory,&ctx->ss_history);
98: nhistory = PetscMax(ctx->rescale_history,1);
99: PetscMalloc1(nhistory,&ctx->yy_rhistory);
100: PetscMalloc1(nhistory,&ctx->ys_rhistory);
101: PetscMalloc1(nhistory,&ctx->ss_rhistory);
103: /* Finish initializations */
104: ctx->lmnow = 0;
105: ctx->iter = 0;
106: ctx->nupdates = 0;
107: ctx->nrejects = 0;
108: ctx->delta = 1.0;
110: ctx->Gprev = 0;
111: ctx->Xprev = 0;
113: ctx->scale = 0;
114: ctx->useScale = PETSC_FALSE;
116: ctx->H0_mat = 0;
117: ctx->H0_ksp = 0;
118: ctx->H0_norm = 0;
119: ctx->useDefaultH0 = PETSC_TRUE;
121: MatCreateShell(comm, n, n, N, N, ctx, A);
122: MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
123: MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
124: return(0);
125: }
129: extern 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: }
199: extern PetscErrorCode MatView_LMVM(Mat A, PetscViewer pv)
200: {
201: PetscBool isascii;
203: MatLMVMCtx *lmP;
206: MatShellGetContext(A,(void**)&lmP);
207: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
208: if (isascii) {
209: PetscViewerASCIIPrintf(pv,"LMVM Matrix\n");
210: PetscViewerASCIIPrintf(pv," Number of vectors: %D\n",lmP->lm);
211: PetscViewerASCIIPrintf(pv," scale type: %s\n",Scale_Table[lmP->scaleType]);
212: PetscViewerASCIIPrintf(pv," rescale type: %s\n",Rescale_Table[lmP->rScaleType]);
213: PetscViewerASCIIPrintf(pv," limit type: %s\n",Limit_Table[lmP->limitType]);
214: PetscViewerASCIIPrintf(pv," updates: %D\n",lmP->nupdates);
215: PetscViewerASCIIPrintf(pv," rejects: %D\n",lmP->nrejects);
216: }
217: return(0);
218: }
222: extern PetscErrorCode MatDestroy_LMVM(Mat M)
223: {
224: MatLMVMCtx *ctx;
228: MatShellGetContext(M,(void**)&ctx);
229: if (ctx->allocated) {
230: if (ctx->Xprev) {
231: PetscObjectDereference((PetscObject)ctx->Xprev);
232: }
233: if (ctx->Gprev) {
234: PetscObjectDereference((PetscObject)ctx->Gprev);
235: }
237: VecDestroyVecs(ctx->lm+1,&ctx->S);
238: VecDestroyVecs(ctx->lm+1,&ctx->Y);
239: VecDestroy(&ctx->D);
240: VecDestroy(&ctx->U);
241: VecDestroy(&ctx->V);
242: VecDestroy(&ctx->W);
243: VecDestroy(&ctx->P);
244: VecDestroy(&ctx->Q);
245: VecDestroy(&ctx->H0_norm);
246: if (ctx->scale) {
247: VecDestroy(&ctx->scale);
248: }
249: }
250: if (ctx->H0_mat) {
251: PetscObjectDereference((PetscObject)ctx->H0_mat);
252: KSPDestroy(&ctx->H0_ksp);
253: }
254: PetscFree(ctx->rho);
255: PetscFree(ctx->beta);
256: PetscFree(ctx->yy_history);
257: PetscFree(ctx->ys_history);
258: PetscFree(ctx->ss_history);
259: PetscFree(ctx->yy_rhistory);
260: PetscFree(ctx->ys_rhistory);
261: PetscFree(ctx->ss_rhistory);
262: PetscFree(ctx);
263: return(0);
264: }
268: extern PetscErrorCode MatLMVMReset(Mat M)
269: {
271: MatLMVMCtx *ctx;
272: PetscInt i;
275: MatShellGetContext(M,(void**)&ctx);
276: if (ctx->Gprev) {
277: PetscObjectDereference((PetscObject)ctx->Gprev);
278: }
279: if (ctx->Xprev) {
280: PetscObjectDereference((PetscObject)ctx->Xprev);
281: }
282: ctx->Gprev = ctx->Y[ctx->lm];
283: ctx->Xprev = ctx->S[ctx->lm];
284: PetscObjectReference((PetscObject)ctx->Gprev);
285: PetscObjectReference((PetscObject)ctx->Xprev);
286: for (i=0; i<ctx->lm; ++i) {
287: ctx->rho[i] = 0.0;
288: }
289: ctx->rho[0] = 1.0;
291: /* Set the scaling and diagonal scaling matrix */
292: switch(ctx->scaleType) {
293: case MatLMVM_Scale_None:
294: ctx->sigma = 1.0;
295: break;
296: case MatLMVM_Scale_Scalar:
297: ctx->sigma = ctx->delta;
298: break;
299: case MatLMVM_Scale_Broyden:
300: VecSet(ctx->D,ctx->delta);
301: break;
302: }
304: ctx->iter=0;
305: ctx->nupdates=0;
306: ctx->lmnow=0;
307: return(0);
308: }
312: extern PetscErrorCode MatLMVMUpdate(Mat M, Vec x, Vec g)
313: {
314: MatLMVMCtx *ctx;
315: PetscReal rhotemp, rhotol;
316: PetscReal y0temp, s0temp;
317: PetscReal yDy, yDs, sDs;
318: PetscReal sigmanew, denom;
320: PetscInt i;
321: PetscBool same;
322: PetscReal yy_sum=0.0, ys_sum=0.0, ss_sum=0.0;
327: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
328: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
329: MatShellGetContext(M,(void**)&ctx);
330: if (!ctx->allocated) {
331: MatLMVMAllocateVectors(M, x);
332: }
334: if (0 == ctx->iter) {
335: MatLMVMReset(M);
336: } else {
337: VecAYPX(ctx->Gprev,-1.0,g);
338: VecAYPX(ctx->Xprev,-1.0,x);
340: VecDot(ctx->Gprev,ctx->Xprev,&rhotemp);
341: VecDot(ctx->Gprev,ctx->Gprev,&y0temp);
343: rhotol = ctx->eps * y0temp;
344: if (rhotemp > rhotol) {
345: ++ctx->nupdates;
347: ctx->lmnow = PetscMin(ctx->lmnow+1, ctx->lm);
348: ierr=PetscObjectDereference((PetscObject)ctx->S[ctx->lm]);
349: ierr=PetscObjectDereference((PetscObject)ctx->Y[ctx->lm]);
350: for (i = ctx->lm-1; i >= 0; --i) {
351: ctx->S[i+1] = ctx->S[i];
352: ctx->Y[i+1] = ctx->Y[i];
353: ctx->rho[i+1] = ctx->rho[i];
354: }
355: ctx->S[0] = ctx->Xprev;
356: ctx->Y[0] = ctx->Gprev;
357: PetscObjectReference((PetscObject)ctx->S[0]);
358: PetscObjectReference((PetscObject)ctx->Y[0]);
359: ctx->rho[0] = 1.0 / rhotemp;
361: /* Compute the scaling */
362: switch(ctx->scaleType) {
363: case MatLMVM_Scale_None:
364: break;
366: case MatLMVM_Scale_Scalar:
367: /* Compute s^T s */
368: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
370: /* Scalar is positive; safeguards are not required. */
372: /* Save information for scalar scaling */
373: ctx->yy_history[(ctx->nupdates - 1) % ctx->scalar_history] = y0temp;
374: ctx->ys_history[(ctx->nupdates - 1) % ctx->scalar_history] = rhotemp;
375: ctx->ss_history[(ctx->nupdates - 1) % ctx->scalar_history] = s0temp;
377: /* Compute summations for scalar scaling */
378: yy_sum = 0; /* No safeguard required; y^T y > 0 */
379: ys_sum = 0; /* No safeguard required; y^T s > 0 */
380: ss_sum = 0; /* No safeguard required; s^T s > 0 */
381: for (i = 0; i < PetscMin(ctx->nupdates, ctx->scalar_history); ++i) {
382: yy_sum += ctx->yy_history[i];
383: ys_sum += ctx->ys_history[i];
384: ss_sum += ctx->ss_history[i];
385: }
387: if (0.0 == ctx->s_alpha) {
388: /* Safeguard ys_sum */
389: if (0.0 == ys_sum) {
390: ys_sum = TAO_ZERO_SAFEGUARD;
391: }
393: sigmanew = ss_sum / ys_sum;
394: } else if (1.0 == ctx->s_alpha) {
395: /* Safeguard yy_sum */
396: if (0.0 == yy_sum) {
397: yy_sum = TAO_ZERO_SAFEGUARD;
398: }
400: sigmanew = ys_sum / yy_sum;
401: } else {
402: denom = 2*ctx->s_alpha*yy_sum;
404: /* Safeguard denom */
405: if (0.0 == denom) {
406: denom = TAO_ZERO_SAFEGUARD;
407: }
409: 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;
410: }
412: switch(ctx->limitType) {
413: case MatLMVM_Limit_Average:
414: if (1.0 == ctx->mu) {
415: ctx->sigma = sigmanew;
416: } else if (ctx->mu) {
417: ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
418: }
419: break;
421: case MatLMVM_Limit_Relative:
422: if (ctx->mu) {
423: ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
424: }
425: break;
427: case MatLMVM_Limit_Absolute:
428: if (ctx->nu) {
429: ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
430: }
431: break;
433: default:
434: ctx->sigma = sigmanew;
435: break;
436: }
437: break;
439: case MatLMVM_Scale_Broyden:
440: /* Original version */
441: /* Combine DFP and BFGS */
443: /* This code appears to be numerically unstable. We use the */
444: /* original version because this was used to generate all of */
445: /* the data and because it may be the least unstable of the */
446: /* bunch. */
448: /* P = Q = inv(D); */
449: VecCopy(ctx->D,ctx->P);
450: VecReciprocal(ctx->P);
451: VecCopy(ctx->P,ctx->Q);
453: /* V = y*y */
454: VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev);
456: /* W = inv(D)*s */
457: VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P);
458: VecDot(ctx->W,ctx->Xprev,&sDs);
460: /* Safeguard rhotemp and sDs */
461: if (0.0 == rhotemp) {
462: rhotemp = TAO_ZERO_SAFEGUARD;
463: }
465: if (0.0 == sDs) {
466: sDs = TAO_ZERO_SAFEGUARD;
467: }
469: if (1.0 != ctx->phi) {
470: /* BFGS portion of the update */
471: /* U = (inv(D)*s)*(inv(D)*s) */
472: VecPointwiseMult(ctx->U,ctx->W,ctx->W);
474: /* Assemble */
475: VecAXPY(ctx->P,1.0/rhotemp,ctx->V);
476: VecAXPY(ctx->P,-1.0/sDs,ctx->U);
477: }
479: if (0.0 != ctx->phi) {
480: /* DFP portion of the update */
481: /* U = inv(D)*s*y */
482: VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev);
484: /* Assemble */
485: VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V);
486: VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U);
487: }
489: if (0.0 == ctx->phi) {
490: VecCopy(ctx->P,ctx->U);
491: } else if (1.0 == ctx->phi) {
492: VecCopy(ctx->Q,ctx->U);
493: } else {
494: /* Broyden update U=(1-phi)*P + phi*Q */
495: VecCopy(ctx->Q,ctx->U);
496: VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P);
497: }
499: /* Obtain inverse and ensure positive definite */
500: VecReciprocal(ctx->U);
501: VecAbs(ctx->U);
503: switch(ctx->rScaleType) {
504: case MatLMVM_Rescale_None:
505: break;
507: case MatLMVM_Rescale_Scalar:
508: case MatLMVM_Rescale_GL:
509: if (ctx->rScaleType == MatLMVM_Rescale_GL) {
510: /* Gilbert and Lemarachal use the old diagonal */
511: VecCopy(ctx->D,ctx->P);
512: } else {
513: /* The default version uses the current diagonal */
514: VecCopy(ctx->U,ctx->P);
515: }
517: /* Compute s^T s */
518: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
520: /* Save information for special cases of scalar rescaling */
521: ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
522: ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
523: ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;
525: if (0.5 == ctx->r_beta) {
526: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
527: VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P);
528: VecDot(ctx->V,ctx->Y[0],&yy_sum);
530: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
531: VecDot(ctx->W,ctx->S[0],&ss_sum);
533: ys_sum = ctx->ys_rhistory[0];
534: } else {
535: VecCopy(ctx->P,ctx->Q);
536: VecReciprocal(ctx->Q);
538: /* Compute summations for scalar scaling */
539: yy_sum = 0; /* No safeguard required */
540: ys_sum = 0; /* No safeguard required */
541: ss_sum = 0; /* No safeguard required */
542: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
543: VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P);
544: VecDot(ctx->V,ctx->Y[i],&yDy);
545: yy_sum += yDy;
547: VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q);
548: VecDot(ctx->W,ctx->S[i],&sDs);
549: ss_sum += sDs;
550: ys_sum += ctx->ys_rhistory[i];
551: }
552: }
553: } else if (0.0 == ctx->r_beta) {
554: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
555: /* Compute summations for scalar scaling */
556: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
558: VecDot(ctx->W, ctx->Y[0], &ys_sum);
559: VecDot(ctx->W, ctx->W, &ss_sum);
560: yy_sum += ctx->yy_rhistory[0];
561: } else {
562: VecCopy(ctx->Q, ctx->P);
563: VecReciprocal(ctx->Q);
565: /* Compute summations for scalar scaling */
566: yy_sum = 0; /* No safeguard required */
567: ys_sum = 0; /* No safeguard required */
568: ss_sum = 0; /* No safeguard required */
569: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
570: VecPointwiseMult(ctx->W, ctx->S[i], ctx->Q);
571: VecDot(ctx->W, ctx->Y[i], &yDs);
572: ys_sum += yDs;
574: VecDot(ctx->W, ctx->W, &sDs);
575: ss_sum += sDs;
577: yy_sum += ctx->yy_rhistory[i];
578: }
579: }
580: } else if (1.0 == ctx->r_beta) {
581: /* Compute summations for scalar scaling */
582: yy_sum = 0; /* No safeguard required */
583: ys_sum = 0; /* No safeguard required */
584: ss_sum = 0; /* No safeguard required */
585: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
586: VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P);
587: VecDot(ctx->V, ctx->S[i], &yDs);
588: ys_sum += yDs;
590: VecDot(ctx->V, ctx->V, &yDy);
591: yy_sum += yDy;
593: ss_sum += ctx->ss_rhistory[i];
594: }
595: } else {
596: VecCopy(ctx->Q, ctx->P);
598: VecPow(ctx->P, ctx->r_beta);
599: VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q);
601: /* Compute summations for scalar scaling */
602: yy_sum = 0; /* No safeguard required */
603: ys_sum = 0; /* No safeguard required */
604: ss_sum = 0; /* No safeguard required */
605: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
606: VecPointwiseMult(ctx->V, ctx->P, ctx->Y[i]);
607: VecPointwiseMult(ctx->W, ctx->Q, ctx->S[i]);
609: VecDot(ctx->V, ctx->V, &yDy);
610: VecDot(ctx->V, ctx->W, &yDs);
611: VecDot(ctx->W, ctx->W, &sDs);
613: yy_sum += yDy;
614: ys_sum += yDs;
615: ss_sum += sDs;
616: }
617: }
619: if (0.0 == ctx->r_alpha) {
620: /* Safeguard ys_sum */
621: if (0.0 == ys_sum) {
622: ys_sum = TAO_ZERO_SAFEGUARD;
623: }
625: sigmanew = ss_sum / ys_sum;
626: } else if (1.0 == ctx->r_alpha) {
627: /* Safeguard yy_sum */
628: if (0.0 == yy_sum) {
629: ys_sum = TAO_ZERO_SAFEGUARD;
630: }
632: sigmanew = ys_sum / yy_sum;
633: } else {
634: denom = 2*ctx->r_alpha*yy_sum;
636: /* Safeguard denom */
637: if (0.0 == denom) {
638: denom = TAO_ZERO_SAFEGUARD;
639: }
641: 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;
642: }
644: /* If Q has small values, then Q^(r_beta - 1) */
645: /* can have very large values. Hence, ys_sum */
646: /* and ss_sum can be infinity. In this case, */
647: /* sigmanew can either be not-a-number or infinity. */
649: if (PetscIsInfOrNanReal(sigmanew)) {
650: /* sigmanew is not-a-number; skip rescaling */
651: } else if (!sigmanew) {
652: /* sigmanew is zero; this is a bad case; skip rescaling */
653: } else {
654: /* sigmanew is positive */
655: VecScale(ctx->U, sigmanew);
656: }
657: break;
658: }
660: /* Modify for previous information */
661: switch(ctx->limitType) {
662: case MatLMVM_Limit_Average:
663: if (1.0 == ctx->mu) {
664: VecCopy(ctx->D, ctx->U);
665: } else if (ctx->mu) {
666: VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U);
667: }
668: break;
670: case MatLMVM_Limit_Relative:
671: if (ctx->mu) {
672: /* P = (1-mu) * D */
673: VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D);
674: /* Q = (1+mu) * D */
675: VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D);
676: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D);
677: }
678: break;
680: case MatLMVM_Limit_Absolute:
681: if (ctx->nu) {
682: VecCopy(ctx->P, ctx->D);
683: VecShift(ctx->P, -ctx->nu);
684: VecCopy(ctx->D, ctx->Q);
685: VecShift(ctx->Q, ctx->nu);
686: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P);
687: }
688: break;
690: default:
691: VecCopy(ctx->U, ctx->D);
692: break;
693: }
694: break;
695: }
696: PetscObjectDereference((PetscObject)ctx->Xprev);
697: PetscObjectDereference((PetscObject)ctx->Gprev);
698: ctx->Xprev = ctx->S[ctx->lm];
699: ctx->Gprev = ctx->Y[ctx->lm];
700: PetscObjectReference((PetscObject)ctx->S[ctx->lm]);
701: PetscObjectReference((PetscObject)ctx->Y[ctx->lm]);
703: } else {
704: ++ctx->nrejects;
705: }
706: }
708: ++ctx->iter;
709: VecCopy(x, ctx->Xprev);
710: VecCopy(g, ctx->Gprev);
711: return(0);
712: }
716: extern PetscErrorCode MatLMVMSetDelta(Mat m, PetscReal d)
717: {
718: MatLMVMCtx *ctx;
720: PetscBool same;
724: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
725: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
726: MatShellGetContext(m,(void**)&ctx);
727: ctx->delta = PetscAbsReal(d);
728: ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
729: ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
730: return(0);
731: }
735: extern PetscErrorCode MatLMVMSetScale(Mat m, Vec s)
736: {
737: MatLMVMCtx *ctx;
739: PetscBool same;
743: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
744: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
745: MatShellGetContext(m,(void**)&ctx);
747: VecDestroy(&ctx->scale);
748: if (s) {
749: VecDuplicate(s,&ctx->scale);
750: }
751: return(0);
752: }
756: extern PetscErrorCode MatLMVMGetRejects(Mat m, PetscInt *nrejects)
757: {
758: MatLMVMCtx *ctx;
760: PetscBool same;
764: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
765: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
766: MatShellGetContext(m,(void**)&ctx);
767: *nrejects = ctx->nrejects;
768: return(0);
769: }
773: extern PetscErrorCode MatLMVMSetH0(Mat m, Mat H0)
774: {
775: MatLMVMCtx *ctx;
777: PetscBool same;
781: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
782: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
783: MatShellGetContext(m,(void**)&ctx);
785: ctx->H0_mat = H0;
786: PetscObjectReference((PetscObject)ctx->H0_mat);
788: ctx->useDefaultH0 = PETSC_FALSE;
790: KSPCreate(PetscObjectComm((PetscObject)H0), &ctx->H0_ksp);
791: KSPSetOperators(ctx->H0_ksp, H0, H0);
792: /* its options prefix and setup is handled in TaoSolve_LMVM/TaoSolve_BLMVM */
793: return(0);
794: }
798: extern PetscErrorCode MatLMVMGetH0(Mat m, Mat *H0)
799: {
800: MatLMVMCtx *ctx;
802: PetscBool same;
806: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
807: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
809: MatShellGetContext(m,(void**)&ctx);
810: *H0 = ctx->H0_mat;
811: return(0);
812: }
816: extern PetscErrorCode MatLMVMGetH0KSP(Mat m, KSP *H0ksp)
817: {
818: MatLMVMCtx *ctx;
820: PetscBool same;
824: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
825: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
827: MatShellGetContext(m,(void**)&ctx);
828: *H0ksp = ctx->H0_ksp;
829: return(0);
830: }
834: extern PetscErrorCode MatLMVMGetX0(Mat m, Vec x)
835: {
837: return(0);
838: }
842: extern PetscErrorCode MatLMVMSetPrev(Mat M, Vec x, Vec g)
843: {
844: MatLMVMCtx *ctx;
846: PetscBool same;
851: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
852: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
853: MatShellGetContext(M,(void**)&ctx);
854: if (ctx->nupdates == 0) {
855: MatLMVMUpdate(M,x,g);
856: } else {
857: VecCopy(x,ctx->Xprev);
858: VecCopy(g,ctx->Gprev);
859: /* TODO scaling specific terms */
860: }
861: return(0);
862: }
866: extern PetscErrorCode MatLMVMRefine(Mat coarse, Mat op, Vec fineX, Vec fineG)
867: {
869: PetscBool same;
876: PetscObjectTypeCompare((PetscObject)coarse,MATSHELL,&same);
877: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
878: PetscObjectTypeCompare((PetscObject)op,MATSHELL,&same);
879: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
880: return(0);
881: }
885: extern PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v)
886: {
888: MatLMVMCtx *ctx;
889: PetscBool same;
894: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
895: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
896: MatShellGetContext(m,(void**)&ctx);
898: /* Perform allocations */
899: VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
900: VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
901: VecDuplicate(v,&ctx->D);
902: VecDuplicate(v,&ctx->U);
903: VecDuplicate(v,&ctx->V);
904: VecDuplicate(v,&ctx->W);
905: VecDuplicate(v,&ctx->P);
906: VecDuplicate(v,&ctx->Q);
907: VecDuplicate(v,&ctx->H0_norm);
908: ctx->allocated = PETSC_TRUE;
909: return(0);
910: }