Actual source code: lmvmmat.c
petsc-3.8.4 2018-03-24
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: extern 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;
59:
60: ctx->phi = 0.125;
61:
62: ctx->scalar_history = 1;
63: ctx->rescale_history = 1;
64:
65: ctx->delta_min = 1e-7;
66: ctx->delta_max = 100.0;
68: /* Begin configuration */
69: PetscOptionsBegin(comm,NULL,NULL,NULL);
70: PetscOptionsInt("-tao_lmm_vectors", "vectors to use for approximation", "", ctx->lm, &ctx->lm,NULL);
71: PetscOptionsReal("-tao_lmm_limit_mu", "mu limiting factor", "", ctx->mu, &ctx->mu,NULL);
72: PetscOptionsReal("-tao_lmm_limit_nu", "nu limiting factor", "", ctx->nu, &ctx->nu,NULL);
73: PetscOptionsReal("-tao_lmm_broyden_phi", "phi factor for Broyden scaling", "", ctx->phi, &ctx->phi,NULL);
74: PetscOptionsReal("-tao_lmm_scalar_alpha", "alpha factor for scalar scaling", "",ctx->s_alpha, &ctx->s_alpha,NULL);
75: PetscOptionsReal("-tao_lmm_rescale_alpha", "alpha factor for rescaling diagonal", "", ctx->r_alpha, &ctx->r_alpha,NULL);
76: PetscOptionsReal("-tao_lmm_rescale_beta", "beta factor for rescaling diagonal", "", ctx->r_beta, &ctx->r_beta,NULL);
77: PetscOptionsInt("-tao_lmm_scalar_history", "amount of history for scalar scaling", "", ctx->scalar_history, &ctx->scalar_history,NULL);
78: PetscOptionsInt("-tao_lmm_rescale_history", "amount of history for rescaling diagonal", "", ctx->rescale_history, &ctx->rescale_history,NULL);
79: PetscOptionsReal("-tao_lmm_eps", "rejection tolerance", "", ctx->eps, &ctx->eps,NULL);
80: PetscOptionsEList("-tao_lmm_scale_type", "scale type", "", Scale_Table, MatLMVM_Scale_Types, Scale_Table[ctx->scaleType], &ctx->scaleType,NULL);
81: PetscOptionsEList("-tao_lmm_rescale_type", "rescale type", "", Rescale_Table, MatLMVM_Rescale_Types, Rescale_Table[ctx->rScaleType], &ctx->rScaleType,NULL);
82: PetscOptionsEList("-tao_lmm_limit_type", "limit type", "", Limit_Table, MatLMVM_Limit_Types, Limit_Table[ctx->limitType], &ctx->limitType,NULL);
83: PetscOptionsReal("-tao_lmm_delta_min", "minimum delta value", "", ctx->delta_min, &ctx->delta_min,NULL);
84: PetscOptionsReal("-tao_lmm_delta_max", "maximum delta value", "", ctx->delta_max, &ctx->delta_max,NULL);
85: PetscOptionsEnd();
87: /* Complete configuration */
88: ctx->rescale_history = PetscMin(ctx->rescale_history, ctx->lm);
89: PetscMalloc1(ctx->lm+1,&ctx->rho);
90: PetscMalloc1(ctx->lm+1,&ctx->beta);
92: nhistory = PetscMax(ctx->scalar_history,1);
93: PetscMalloc1(nhistory,&ctx->yy_history);
94: PetscMalloc1(nhistory,&ctx->ys_history);
95: PetscMalloc1(nhistory,&ctx->ss_history);
97: nhistory = PetscMax(ctx->rescale_history,1);
98: PetscMalloc1(nhistory,&ctx->yy_rhistory);
99: PetscMalloc1(nhistory,&ctx->ys_rhistory);
100: PetscMalloc1(nhistory,&ctx->ss_rhistory);
102: /* Finish initializations */
103: ctx->lmnow = 0;
104: ctx->iter = 0;
105: ctx->nupdates = 0;
106: ctx->nrejects = 0;
107: ctx->delta = 1.0;
109: ctx->Gprev = 0;
110: ctx->Xprev = 0;
112: ctx->scale = 0;
113: ctx->useScale = PETSC_FALSE;
115: ctx->H0_mat = 0;
116: ctx->H0_ksp = 0;
117: ctx->H0_norm = 0;
118: ctx->useDefaultH0 = PETSC_TRUE;
120: MatCreateShell(comm, n, n, N, N, ctx, A);
121: MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
122: MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
123: return(0);
124: }
126: extern PetscErrorCode MatLMVMSolve(Mat A, Vec b, Vec x)
127: {
128: PetscReal sq, yq, dd;
129: PetscInt ll;
130: PetscBool scaled;
131: MatLMVMCtx *shell;
138: MatShellGetContext(A,(void**)&shell);
139: if (shell->lmnow < 1) {
140: shell->rho[0] = 1.0;
141: shell->theta = 1.0;
142: }
144: VecCopy(b,x);
145: for (ll = 0; ll < shell->lmnow; ++ll) {
146: VecDot(x,shell->S[ll],&sq);
147: shell->beta[ll] = sq * shell->rho[ll];
148: VecAXPY(x,-shell->beta[ll],shell->Y[ll]);
149: }
151: scaled = PETSC_FALSE;
152: if (!scaled && !shell->useDefaultH0 && shell->H0_mat) {
153: KSPSolve(shell->H0_ksp,x,shell->U);
154: VecScale(shell->U, shell->theta);
155: VecDot(x,shell->U,&dd);
156: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
157: /* Accept Hessian solve */
158: VecCopy(shell->U,x);
159: scaled = PETSC_TRUE;
160: }
161: }
163: if (!scaled && shell->useScale) {
164: VecPointwiseMult(shell->U,x,shell->scale);
165: VecDot(x,shell->U,&dd);
166: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
167: /* Accept scaling */
168: VecCopy(shell->U,x);
169: scaled = PETSC_TRUE;
170: }
171: }
173: if (!scaled) {
174: switch(shell->scaleType) {
175: case MatLMVM_Scale_None:
176: break;
178: case MatLMVM_Scale_Scalar:
179: VecScale(x,shell->sigma);
180: break;
182: case MatLMVM_Scale_Broyden:
183: VecPointwiseMult(x,x,shell->D);
184: break;
185: }
186: }
187: for (ll = shell->lmnow-1; ll >= 0; --ll) {
188: VecDot(x,shell->Y[ll],&yq);
189: VecAXPY(x,shell->beta[ll]-yq*shell->rho[ll],shell->S[ll]);
190: }
191: return(0);
192: }
194: extern PetscErrorCode MatView_LMVM(Mat A, PetscViewer pv)
195: {
196: PetscBool isascii;
198: MatLMVMCtx *lmP;
201: MatShellGetContext(A,(void**)&lmP);
202: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
203: if (isascii) {
204: PetscViewerASCIIPrintf(pv,"LMVM Matrix\n");
205: PetscViewerASCIIPrintf(pv," Number of vectors: %D\n",lmP->lm);
206: PetscViewerASCIIPrintf(pv," scale type: %s\n",Scale_Table[lmP->scaleType]);
207: PetscViewerASCIIPrintf(pv," rescale type: %s\n",Rescale_Table[lmP->rScaleType]);
208: PetscViewerASCIIPrintf(pv," limit type: %s\n",Limit_Table[lmP->limitType]);
209: PetscViewerASCIIPrintf(pv," updates: %D\n",lmP->nupdates);
210: PetscViewerASCIIPrintf(pv," rejects: %D\n",lmP->nrejects);
211: }
212: return(0);
213: }
215: extern PetscErrorCode MatDestroy_LMVM(Mat M)
216: {
217: MatLMVMCtx *ctx;
221: MatShellGetContext(M,(void**)&ctx);
222: if (ctx->allocated) {
223: if (ctx->Xprev) {
224: PetscObjectDereference((PetscObject)ctx->Xprev);
225: }
226: if (ctx->Gprev) {
227: PetscObjectDereference((PetscObject)ctx->Gprev);
228: }
230: VecDestroyVecs(ctx->lm+1,&ctx->S);
231: VecDestroyVecs(ctx->lm+1,&ctx->Y);
232: VecDestroy(&ctx->D);
233: VecDestroy(&ctx->U);
234: VecDestroy(&ctx->V);
235: VecDestroy(&ctx->W);
236: VecDestroy(&ctx->P);
237: VecDestroy(&ctx->Q);
238: VecDestroy(&ctx->H0_norm);
239: if (ctx->scale) {
240: VecDestroy(&ctx->scale);
241: }
242: }
243: if (ctx->H0_mat) {
244: PetscObjectDereference((PetscObject)ctx->H0_mat);
245: KSPDestroy(&ctx->H0_ksp);
246: }
247: PetscFree(ctx->rho);
248: PetscFree(ctx->beta);
249: PetscFree(ctx->yy_history);
250: PetscFree(ctx->ys_history);
251: PetscFree(ctx->ss_history);
252: PetscFree(ctx->yy_rhistory);
253: PetscFree(ctx->ys_rhistory);
254: PetscFree(ctx->ss_rhistory);
255: PetscFree(ctx);
256: return(0);
257: }
259: extern PetscErrorCode MatLMVMReset(Mat M)
260: {
262: MatLMVMCtx *ctx;
263: PetscInt i;
266: MatShellGetContext(M,(void**)&ctx);
267: if (ctx->Gprev) {
268: PetscObjectDereference((PetscObject)ctx->Gprev);
269: }
270: if (ctx->Xprev) {
271: PetscObjectDereference((PetscObject)ctx->Xprev);
272: }
273: ctx->Gprev = ctx->Y[ctx->lm];
274: ctx->Xprev = ctx->S[ctx->lm];
275: PetscObjectReference((PetscObject)ctx->Gprev);
276: PetscObjectReference((PetscObject)ctx->Xprev);
277: for (i=0; i<ctx->lm; ++i) {
278: ctx->rho[i] = 0.0;
279: }
280: ctx->rho[0] = 1.0;
282: /* Set the scaling and diagonal scaling matrix */
283: switch(ctx->scaleType) {
284: case MatLMVM_Scale_None:
285: ctx->sigma = 1.0;
286: break;
287: case MatLMVM_Scale_Scalar:
288: ctx->sigma = ctx->delta;
289: break;
290: case MatLMVM_Scale_Broyden:
291: VecSet(ctx->D,ctx->delta);
292: break;
293: }
295: ctx->iter=0;
296: ctx->nupdates=0;
297: ctx->lmnow=0;
298: return(0);
299: }
301: extern PetscErrorCode MatLMVMUpdate(Mat M, Vec x, Vec g)
302: {
303: MatLMVMCtx *ctx;
304: PetscReal rhotemp, rhotol;
305: PetscReal y0temp, s0temp;
306: PetscReal yDy, yDs, sDs;
307: PetscReal sigmanew, denom;
309: PetscInt i;
310: PetscBool same;
311: PetscReal yy_sum=0.0, ys_sum=0.0, ss_sum=0.0;
316: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
317: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
318: MatShellGetContext(M,(void**)&ctx);
319: if (!ctx->allocated) {
320: MatLMVMAllocateVectors(M, x);
321: }
323: if (0 == ctx->iter) {
324: MatLMVMReset(M);
325: } else {
326: VecAYPX(ctx->Gprev,-1.0,g);
327: VecAYPX(ctx->Xprev,-1.0,x);
329: VecDot(ctx->Gprev,ctx->Xprev,&rhotemp);
330: VecDot(ctx->Gprev,ctx->Gprev,&y0temp);
332: rhotol = ctx->eps * y0temp;
333: if (rhotemp > rhotol) {
334: ++ctx->nupdates;
336: ctx->lmnow = PetscMin(ctx->lmnow+1, ctx->lm);
337: ierr=PetscObjectDereference((PetscObject)ctx->S[ctx->lm]);
338: ierr=PetscObjectDereference((PetscObject)ctx->Y[ctx->lm]);
339: for (i = ctx->lm-1; i >= 0; --i) {
340: ctx->S[i+1] = ctx->S[i];
341: ctx->Y[i+1] = ctx->Y[i];
342: ctx->rho[i+1] = ctx->rho[i];
343: }
344: ctx->S[0] = ctx->Xprev;
345: ctx->Y[0] = ctx->Gprev;
346: PetscObjectReference((PetscObject)ctx->S[0]);
347: PetscObjectReference((PetscObject)ctx->Y[0]);
348: ctx->rho[0] = 1.0 / rhotemp;
350: /* Compute the scaling */
351: switch(ctx->scaleType) {
352: case MatLMVM_Scale_None:
353: break;
355: case MatLMVM_Scale_Scalar:
356: /* Compute s^T s */
357: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
359: /* Scalar is positive; safeguards are not required. */
361: /* Save information for scalar scaling */
362: ctx->yy_history[(ctx->nupdates - 1) % ctx->scalar_history] = y0temp;
363: ctx->ys_history[(ctx->nupdates - 1) % ctx->scalar_history] = rhotemp;
364: ctx->ss_history[(ctx->nupdates - 1) % ctx->scalar_history] = s0temp;
366: /* Compute summations for scalar scaling */
367: yy_sum = 0; /* No safeguard required; y^T y > 0 */
368: ys_sum = 0; /* No safeguard required; y^T s > 0 */
369: ss_sum = 0; /* No safeguard required; s^T s > 0 */
370: for (i = 0; i < PetscMin(ctx->nupdates, ctx->scalar_history); ++i) {
371: yy_sum += ctx->yy_history[i];
372: ys_sum += ctx->ys_history[i];
373: ss_sum += ctx->ss_history[i];
374: }
376: if (0.0 == ctx->s_alpha) {
377: /* Safeguard ys_sum */
378: if (0.0 == ys_sum) {
379: ys_sum = TAO_ZERO_SAFEGUARD;
380: }
382: sigmanew = ss_sum / ys_sum;
383: } else if (1.0 == ctx->s_alpha) {
384: /* Safeguard yy_sum */
385: if (0.0 == yy_sum) {
386: yy_sum = TAO_ZERO_SAFEGUARD;
387: }
389: sigmanew = ys_sum / yy_sum;
390: } else {
391: denom = 2*ctx->s_alpha*yy_sum;
393: /* Safeguard denom */
394: if (0.0 == denom) {
395: denom = TAO_ZERO_SAFEGUARD;
396: }
398: 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;
399: }
401: switch(ctx->limitType) {
402: case MatLMVM_Limit_Average:
403: if (1.0 == ctx->mu) {
404: ctx->sigma = sigmanew;
405: } else if (ctx->mu) {
406: ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
407: }
408: break;
410: case MatLMVM_Limit_Relative:
411: if (ctx->mu) {
412: ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
413: }
414: break;
416: case MatLMVM_Limit_Absolute:
417: if (ctx->nu) {
418: ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
419: }
420: break;
422: default:
423: ctx->sigma = sigmanew;
424: break;
425: }
426: break;
428: case MatLMVM_Scale_Broyden:
429: /* Original version */
430: /* Combine DFP and BFGS */
432: /* This code appears to be numerically unstable. We use the */
433: /* original version because this was used to generate all of */
434: /* the data and because it may be the least unstable of the */
435: /* bunch. */
437: /* P = Q = inv(D); */
438: VecCopy(ctx->D,ctx->P);
439: VecReciprocal(ctx->P);
440: VecCopy(ctx->P,ctx->Q);
442: /* V = y*y */
443: VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev);
445: /* W = inv(D)*s */
446: VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P);
447: VecDot(ctx->W,ctx->Xprev,&sDs);
449: /* Safeguard rhotemp and sDs */
450: if (0.0 == rhotemp) {
451: rhotemp = TAO_ZERO_SAFEGUARD;
452: }
454: if (0.0 == sDs) {
455: sDs = TAO_ZERO_SAFEGUARD;
456: }
458: if (1.0 != ctx->phi) {
459: /* BFGS portion of the update */
460: /* U = (inv(D)*s)*(inv(D)*s) */
461: VecPointwiseMult(ctx->U,ctx->W,ctx->W);
463: /* Assemble */
464: VecAXPY(ctx->P,1.0/rhotemp,ctx->V);
465: VecAXPY(ctx->P,-1.0/sDs,ctx->U);
466: }
468: if (0.0 != ctx->phi) {
469: /* DFP portion of the update */
470: /* U = inv(D)*s*y */
471: VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev);
473: /* Assemble */
474: VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V);
475: VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U);
476: }
478: if (0.0 == ctx->phi) {
479: VecCopy(ctx->P,ctx->U);
480: } else if (1.0 == ctx->phi) {
481: VecCopy(ctx->Q,ctx->U);
482: } else {
483: /* Broyden update U=(1-phi)*P + phi*Q */
484: VecCopy(ctx->Q,ctx->U);
485: VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P);
486: }
488: /* Obtain inverse and ensure positive definite */
489: VecReciprocal(ctx->U);
490: VecAbs(ctx->U);
492: switch(ctx->rScaleType) {
493: case MatLMVM_Rescale_None:
494: break;
496: case MatLMVM_Rescale_Scalar:
497: case MatLMVM_Rescale_GL:
498: if (ctx->rScaleType == MatLMVM_Rescale_GL) {
499: /* Gilbert and Lemarachal use the old diagonal */
500: VecCopy(ctx->D,ctx->P);
501: } else {
502: /* The default version uses the current diagonal */
503: VecCopy(ctx->U,ctx->P);
504: }
506: /* Compute s^T s */
507: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
509: /* Save information for special cases of scalar rescaling */
510: ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
511: ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
512: ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;
514: if (0.5 == ctx->r_beta) {
515: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
516: VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P);
517: VecDot(ctx->V,ctx->Y[0],&yy_sum);
519: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
520: VecDot(ctx->W,ctx->S[0],&ss_sum);
522: ys_sum = ctx->ys_rhistory[0];
523: } else {
524: VecCopy(ctx->P,ctx->Q);
525: VecReciprocal(ctx->Q);
527: /* Compute summations for scalar scaling */
528: yy_sum = 0; /* No safeguard required */
529: ys_sum = 0; /* No safeguard required */
530: ss_sum = 0; /* No safeguard required */
531: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
532: VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P);
533: VecDot(ctx->V,ctx->Y[i],&yDy);
534: yy_sum += yDy;
536: VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q);
537: VecDot(ctx->W,ctx->S[i],&sDs);
538: ss_sum += sDs;
539: ys_sum += ctx->ys_rhistory[i];
540: }
541: }
542: } else if (0.0 == ctx->r_beta) {
543: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
544: /* Compute summations for scalar scaling */
545: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
547: VecDot(ctx->W, ctx->Y[0], &ys_sum);
548: VecDot(ctx->W, ctx->W, &ss_sum);
549: yy_sum += ctx->yy_rhistory[0];
550: } else {
551: VecCopy(ctx->Q, ctx->P);
552: VecReciprocal(ctx->Q);
554: /* Compute summations for scalar scaling */
555: yy_sum = 0; /* No safeguard required */
556: ys_sum = 0; /* No safeguard required */
557: ss_sum = 0; /* No safeguard required */
558: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
559: VecPointwiseMult(ctx->W, ctx->S[i], ctx->Q);
560: VecDot(ctx->W, ctx->Y[i], &yDs);
561: ys_sum += yDs;
563: VecDot(ctx->W, ctx->W, &sDs);
564: ss_sum += sDs;
566: yy_sum += ctx->yy_rhistory[i];
567: }
568: }
569: } else if (1.0 == ctx->r_beta) {
570: /* Compute summations for scalar scaling */
571: yy_sum = 0; /* No safeguard required */
572: ys_sum = 0; /* No safeguard required */
573: ss_sum = 0; /* No safeguard required */
574: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
575: VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P);
576: VecDot(ctx->V, ctx->S[i], &yDs);
577: ys_sum += yDs;
579: VecDot(ctx->V, ctx->V, &yDy);
580: yy_sum += yDy;
582: ss_sum += ctx->ss_rhistory[i];
583: }
584: } else {
585: VecCopy(ctx->Q, ctx->P);
587: VecPow(ctx->P, ctx->r_beta);
588: VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q);
590: /* Compute summations for scalar scaling */
591: yy_sum = 0; /* No safeguard required */
592: ys_sum = 0; /* No safeguard required */
593: ss_sum = 0; /* No safeguard required */
594: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
595: VecPointwiseMult(ctx->V, ctx->P, ctx->Y[i]);
596: VecPointwiseMult(ctx->W, ctx->Q, ctx->S[i]);
598: VecDot(ctx->V, ctx->V, &yDy);
599: VecDot(ctx->V, ctx->W, &yDs);
600: VecDot(ctx->W, ctx->W, &sDs);
602: yy_sum += yDy;
603: ys_sum += yDs;
604: ss_sum += sDs;
605: }
606: }
608: if (0.0 == ctx->r_alpha) {
609: /* Safeguard ys_sum */
610: if (0.0 == ys_sum) {
611: ys_sum = TAO_ZERO_SAFEGUARD;
612: }
614: sigmanew = ss_sum / ys_sum;
615: } else if (1.0 == ctx->r_alpha) {
616: /* Safeguard yy_sum */
617: if (0.0 == yy_sum) {
618: ys_sum = TAO_ZERO_SAFEGUARD;
619: }
621: sigmanew = ys_sum / yy_sum;
622: } else {
623: denom = 2*ctx->r_alpha*yy_sum;
625: /* Safeguard denom */
626: if (0.0 == denom) {
627: denom = TAO_ZERO_SAFEGUARD;
628: }
630: 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;
631: }
633: /* If Q has small values, then Q^(r_beta - 1) */
634: /* can have very large values. Hence, ys_sum */
635: /* and ss_sum can be infinity. In this case, */
636: /* sigmanew can either be not-a-number or infinity. */
638: if (PetscIsInfOrNanReal(sigmanew)) {
639: /* sigmanew is not-a-number; skip rescaling */
640: } else if (!sigmanew) {
641: /* sigmanew is zero; this is a bad case; skip rescaling */
642: } else {
643: /* sigmanew is positive */
644: VecScale(ctx->U, sigmanew);
645: }
646: break;
647: }
649: /* Modify for previous information */
650: switch(ctx->limitType) {
651: case MatLMVM_Limit_Average:
652: if (1.0 == ctx->mu) {
653: VecCopy(ctx->D, ctx->U);
654: } else if (ctx->mu) {
655: VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U);
656: }
657: break;
659: case MatLMVM_Limit_Relative:
660: if (ctx->mu) {
661: /* P = (1-mu) * D */
662: VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D);
663: /* Q = (1+mu) * D */
664: VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D);
665: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D);
666: }
667: break;
669: case MatLMVM_Limit_Absolute:
670: if (ctx->nu) {
671: VecCopy(ctx->P, ctx->D);
672: VecShift(ctx->P, -ctx->nu);
673: VecCopy(ctx->D, ctx->Q);
674: VecShift(ctx->Q, ctx->nu);
675: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P);
676: }
677: break;
679: default:
680: VecCopy(ctx->U, ctx->D);
681: break;
682: }
683: break;
684: }
685: PetscObjectDereference((PetscObject)ctx->Xprev);
686: PetscObjectDereference((PetscObject)ctx->Gprev);
687: ctx->Xprev = ctx->S[ctx->lm];
688: ctx->Gprev = ctx->Y[ctx->lm];
689: PetscObjectReference((PetscObject)ctx->S[ctx->lm]);
690: PetscObjectReference((PetscObject)ctx->Y[ctx->lm]);
692: } else {
693: ++ctx->nrejects;
694: }
695: }
697: ++ctx->iter;
698: VecCopy(x, ctx->Xprev);
699: VecCopy(g, ctx->Gprev);
700: return(0);
701: }
703: extern PetscErrorCode MatLMVMSetDelta(Mat m, PetscReal d)
704: {
705: MatLMVMCtx *ctx;
707: PetscBool same;
711: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
712: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
713: MatShellGetContext(m,(void**)&ctx);
714: ctx->delta = PetscAbsReal(d);
715: ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
716: ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
717: return(0);
718: }
720: extern PetscErrorCode MatLMVMSetScale(Mat m, Vec s)
721: {
722: MatLMVMCtx *ctx;
724: PetscBool same;
728: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
729: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
730: MatShellGetContext(m,(void**)&ctx);
732: VecDestroy(&ctx->scale);
733: if (s) {
734: VecDuplicate(s,&ctx->scale);
735: }
736: return(0);
737: }
739: extern PetscErrorCode MatLMVMGetRejects(Mat m, PetscInt *nrejects)
740: {
741: MatLMVMCtx *ctx;
743: PetscBool same;
747: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
748: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
749: MatShellGetContext(m,(void**)&ctx);
750: *nrejects = ctx->nrejects;
751: return(0);
752: }
754: extern PetscErrorCode MatLMVMSetH0(Mat m, Mat H0)
755: {
756: MatLMVMCtx *ctx;
758: PetscBool same;
762: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
763: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
764: MatShellGetContext(m,(void**)&ctx);
766: ctx->H0_mat = H0;
767: PetscObjectReference((PetscObject)ctx->H0_mat);
769: ctx->useDefaultH0 = PETSC_FALSE;
771: KSPCreate(PetscObjectComm((PetscObject)H0), &ctx->H0_ksp);
772: KSPSetOperators(ctx->H0_ksp, H0, H0);
773: /* its options prefix and setup is handled in TaoSolve_LMVM/TaoSolve_BLMVM */
774: return(0);
775: }
777: extern PetscErrorCode MatLMVMGetH0(Mat m, Mat *H0)
778: {
779: MatLMVMCtx *ctx;
781: PetscBool same;
785: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
786: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
788: MatShellGetContext(m,(void**)&ctx);
789: *H0 = ctx->H0_mat;
790: return(0);
791: }
793: extern PetscErrorCode MatLMVMGetH0KSP(Mat m, KSP *H0ksp)
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");
804: MatShellGetContext(m,(void**)&ctx);
805: *H0ksp = ctx->H0_ksp;
806: return(0);
807: }
809: extern PetscErrorCode MatLMVMGetX0(Mat m, Vec x)
810: {
812: return(0);
813: }
815: extern PetscErrorCode MatLMVMSetPrev(Mat M, Vec x, Vec g)
816: {
817: MatLMVMCtx *ctx;
819: PetscBool same;
824: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
825: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
826: MatShellGetContext(M,(void**)&ctx);
827: if (ctx->nupdates == 0) {
828: MatLMVMUpdate(M,x,g);
829: } else {
830: VecCopy(x,ctx->Xprev);
831: VecCopy(g,ctx->Gprev);
832: /* TODO scaling specific terms */
833: }
834: return(0);
835: }
837: extern PetscErrorCode MatLMVMRefine(Mat coarse, Mat op, Vec fineX, Vec fineG)
838: {
840: PetscBool same;
847: PetscObjectTypeCompare((PetscObject)coarse,MATSHELL,&same);
848: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
849: PetscObjectTypeCompare((PetscObject)op,MATSHELL,&same);
850: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
851: return(0);
852: }
854: extern PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v)
855: {
857: MatLMVMCtx *ctx;
858: PetscBool same;
863: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
864: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
865: MatShellGetContext(m,(void**)&ctx);
867: /* Perform allocations */
868: VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
869: VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
870: VecDuplicate(v,&ctx->D);
871: VecDuplicate(v,&ctx->U);
872: VecDuplicate(v,&ctx->V);
873: VecDuplicate(v,&ctx->W);
874: VecDuplicate(v,&ctx->P);
875: VecDuplicate(v,&ctx->Q);
876: VecDuplicate(v,&ctx->H0_norm);
877: ctx->allocated = PETSC_TRUE;
878: return(0);
879: }