Actual source code: lmvmmat.c
petsc-3.6.1 2015-08-06
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 = 0;
117: ctx->useDefaultH0=PETSC_TRUE;
119: MatCreateShell(comm, n, n, N, N, ctx, A);
120: MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
121: MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
122: return(0);
123: }
127: extern PetscErrorCode MatLMVMSolve(Mat A, Vec b, Vec x)
128: {
129: PetscReal sq, yq, dd;
130: PetscInt ll;
131: PetscBool scaled;
132: MatLMVMCtx *shell;
139: MatShellGetContext(A,(void**)&shell);
140: if (shell->lmnow < 1) {
141: shell->rho[0] = 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) {
153: MatSolve(shell->H0,x,shell->U);
154: VecDot(x,shell->U,&dd);
155: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
156: /* Accept Hessian solve */
157: VecCopy(shell->U,x);
158: scaled = PETSC_TRUE;
159: }
160: }
162: if (!scaled && shell->useScale) {
163: VecPointwiseMult(shell->U,x,shell->scale);
164: VecDot(x,shell->U,&dd);
165: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
166: /* Accept scaling */
167: VecCopy(shell->U,x);
168: scaled = PETSC_TRUE;
169: }
170: }
172: if (!scaled) {
173: switch(shell->scaleType) {
174: case MatLMVM_Scale_None:
175: break;
177: case MatLMVM_Scale_Scalar:
178: VecScale(x,shell->sigma);
179: break;
181: case MatLMVM_Scale_Broyden:
182: VecPointwiseMult(x,x,shell->D);
183: break;
184: }
185: }
186: for (ll = shell->lmnow-1; ll >= 0; --ll) {
187: VecDot(x,shell->Y[ll],&yq);
188: VecAXPY(x,shell->beta[ll]-yq*shell->rho[ll],shell->S[ll]);
189: }
190: return(0);
191: }
195: extern PetscErrorCode MatView_LMVM(Mat A, PetscViewer pv)
196: {
197: PetscBool isascii;
199: MatLMVMCtx *lmP;
202: MatShellGetContext(A,(void**)&lmP);
203: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
204: if (isascii) {
205: PetscViewerASCIIPrintf(pv,"LMVM Matrix\n");
206: PetscViewerASCIIPrintf(pv," Number of vectors: %D\n",lmP->lm);
207: PetscViewerASCIIPrintf(pv," scale type: %s\n",Scale_Table[lmP->scaleType]);
208: PetscViewerASCIIPrintf(pv," rescale type: %s\n",Rescale_Table[lmP->rScaleType]);
209: PetscViewerASCIIPrintf(pv," limit type: %s\n",Limit_Table[lmP->limitType]);
210: PetscViewerASCIIPrintf(pv," updates: %D\n",lmP->nupdates);
211: PetscViewerASCIIPrintf(pv," rejects: %D\n",lmP->nrejects);
212: }
213: return(0);
214: }
218: extern 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: if (ctx->scale) {
242: VecDestroy(&ctx->scale);
243: }
244: }
245: PetscFree(ctx->rho);
246: PetscFree(ctx->beta);
247: PetscFree(ctx->yy_history);
248: PetscFree(ctx->ys_history);
249: PetscFree(ctx->ss_history);
250: PetscFree(ctx->yy_rhistory);
251: PetscFree(ctx->ys_rhistory);
252: PetscFree(ctx->ss_rhistory);
253: PetscFree(ctx);
254: return(0);
255: }
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: }
303: extern PetscErrorCode MatLMVMUpdate(Mat M, Vec x, Vec g)
304: {
305: MatLMVMCtx *ctx;
306: PetscReal rhotemp, rhotol;
307: PetscReal y0temp, s0temp;
308: PetscReal yDy, yDs, sDs;
309: PetscReal sigmanew, denom;
311: PetscInt i;
312: PetscBool same;
313: PetscReal yy_sum=0.0, ys_sum=0.0, ss_sum=0.0;
318: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
319: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
320: MatShellGetContext(M,(void**)&ctx);
321: if (!ctx->allocated) {
322: MatLMVMAllocateVectors(M, x);
323: }
325: if (0 == ctx->iter) {
326: MatLMVMReset(M);
327: } else {
328: VecAYPX(ctx->Gprev,-1.0,g);
329: VecAYPX(ctx->Xprev,-1.0,x);
331: VecDot(ctx->Gprev,ctx->Xprev,&rhotemp);
332: VecDot(ctx->Gprev,ctx->Gprev,&y0temp);
334: rhotol = ctx->eps * y0temp;
335: if (rhotemp > rhotol) {
336: ++ctx->nupdates;
338: ctx->lmnow = PetscMin(ctx->lmnow+1, ctx->lm);
339: ierr=PetscObjectDereference((PetscObject)ctx->S[ctx->lm]);
340: ierr=PetscObjectDereference((PetscObject)ctx->Y[ctx->lm]);
341: for (i = ctx->lm-1; i >= 0; --i) {
342: ctx->S[i+1] = ctx->S[i];
343: ctx->Y[i+1] = ctx->Y[i];
344: ctx->rho[i+1] = ctx->rho[i];
345: }
346: ctx->S[0] = ctx->Xprev;
347: ctx->Y[0] = ctx->Gprev;
348: PetscObjectReference((PetscObject)ctx->S[0]);
349: PetscObjectReference((PetscObject)ctx->Y[0]);
350: ctx->rho[0] = 1.0 / rhotemp;
352: /* Compute the scaling */
353: switch(ctx->scaleType) {
354: case MatLMVM_Scale_None:
355: break;
357: case MatLMVM_Scale_Scalar:
358: /* Compute s^T s */
359: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
361: /* Scalar is positive; safeguards are not required. */
363: /* Save information for scalar scaling */
364: ctx->yy_history[(ctx->nupdates - 1) % ctx->scalar_history] = y0temp;
365: ctx->ys_history[(ctx->nupdates - 1) % ctx->scalar_history] = rhotemp;
366: ctx->ss_history[(ctx->nupdates - 1) % ctx->scalar_history] = s0temp;
368: /* Compute summations for scalar scaling */
369: yy_sum = 0; /* No safeguard required; y^T y > 0 */
370: ys_sum = 0; /* No safeguard required; y^T s > 0 */
371: ss_sum = 0; /* No safeguard required; s^T s > 0 */
372: for (i = 0; i < PetscMin(ctx->nupdates, ctx->scalar_history); ++i) {
373: yy_sum += ctx->yy_history[i];
374: ys_sum += ctx->ys_history[i];
375: ss_sum += ctx->ss_history[i];
376: }
378: if (0.0 == ctx->s_alpha) {
379: /* Safeguard ys_sum */
380: if (0.0 == ys_sum) {
381: ys_sum = TAO_ZERO_SAFEGUARD;
382: }
384: sigmanew = ss_sum / ys_sum;
385: } else if (1.0 == ctx->s_alpha) {
386: /* Safeguard yy_sum */
387: if (0.0 == yy_sum) {
388: yy_sum = TAO_ZERO_SAFEGUARD;
389: }
391: sigmanew = ys_sum / yy_sum;
392: } else {
393: denom = 2*ctx->s_alpha*yy_sum;
395: /* Safeguard denom */
396: if (0.0 == denom) {
397: denom = TAO_ZERO_SAFEGUARD;
398: }
400: 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;
401: }
403: switch(ctx->limitType) {
404: case MatLMVM_Limit_Average:
405: if (1.0 == ctx->mu) {
406: ctx->sigma = sigmanew;
407: } else if (ctx->mu) {
408: ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
409: }
410: break;
412: case MatLMVM_Limit_Relative:
413: if (ctx->mu) {
414: ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
415: }
416: break;
418: case MatLMVM_Limit_Absolute:
419: if (ctx->nu) {
420: ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
421: }
422: break;
424: default:
425: ctx->sigma = sigmanew;
426: break;
427: }
428: break;
430: case MatLMVM_Scale_Broyden:
431: /* Original version */
432: /* Combine DFP and BFGS */
434: /* This code appears to be numerically unstable. We use the */
435: /* original version because this was used to generate all of */
436: /* the data and because it may be the least unstable of the */
437: /* bunch. */
439: /* P = Q = inv(D); */
440: VecCopy(ctx->D,ctx->P);
441: VecReciprocal(ctx->P);
442: VecCopy(ctx->P,ctx->Q);
444: /* V = y*y */
445: VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev);
447: /* W = inv(D)*s */
448: VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P);
449: VecDot(ctx->W,ctx->Xprev,&sDs);
451: /* Safeguard rhotemp and sDs */
452: if (0.0 == rhotemp) {
453: rhotemp = TAO_ZERO_SAFEGUARD;
454: }
456: if (0.0 == sDs) {
457: sDs = TAO_ZERO_SAFEGUARD;
458: }
460: if (1.0 != ctx->phi) {
461: /* BFGS portion of the update */
462: /* U = (inv(D)*s)*(inv(D)*s) */
463: VecPointwiseMult(ctx->U,ctx->W,ctx->W);
465: /* Assemble */
466: VecAXPY(ctx->P,1.0/rhotemp,ctx->V);
467: VecAXPY(ctx->P,-1.0/sDs,ctx->U);
468: }
470: if (0.0 != ctx->phi) {
471: /* DFP portion of the update */
472: /* U = inv(D)*s*y */
473: VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev);
475: /* Assemble */
476: VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V);
477: VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U);
478: }
480: if (0.0 == ctx->phi) {
481: VecCopy(ctx->P,ctx->U);
482: } else if (1.0 == ctx->phi) {
483: VecCopy(ctx->Q,ctx->U);
484: } else {
485: /* Broyden update U=(1-phi)*P + phi*Q */
486: VecCopy(ctx->Q,ctx->U);
487: VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P);
488: }
490: /* Obtain inverse and ensure positive definite */
491: VecReciprocal(ctx->U);
492: VecAbs(ctx->U);
494: switch(ctx->rScaleType) {
495: case MatLMVM_Rescale_None:
496: break;
498: case MatLMVM_Rescale_Scalar:
499: case MatLMVM_Rescale_GL:
500: if (ctx->rScaleType == MatLMVM_Rescale_GL) {
501: /* Gilbert and Lemarachal use the old diagonal */
502: VecCopy(ctx->D,ctx->P);
503: } else {
504: /* The default version uses the current diagonal */
505: VecCopy(ctx->U,ctx->P);
506: }
508: /* Compute s^T s */
509: VecDot(ctx->Xprev,ctx->Xprev,&s0temp);
511: /* Save information for special cases of scalar rescaling */
512: ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
513: ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
514: ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;
516: if (0.5 == ctx->r_beta) {
517: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
518: VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P);
519: VecDot(ctx->V,ctx->Y[0],&yy_sum);
521: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
522: VecDot(ctx->W,ctx->S[0],&ss_sum);
524: ys_sum = ctx->ys_rhistory[0];
525: } else {
526: VecCopy(ctx->P,ctx->Q);
527: VecReciprocal(ctx->Q);
529: /* Compute summations for scalar scaling */
530: yy_sum = 0; /* No safeguard required */
531: ys_sum = 0; /* No safeguard required */
532: ss_sum = 0; /* No safeguard required */
533: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
534: VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P);
535: VecDot(ctx->V,ctx->Y[i],&yDy);
536: yy_sum += yDy;
538: VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q);
539: VecDot(ctx->W,ctx->S[i],&sDs);
540: ss_sum += sDs;
541: ys_sum += ctx->ys_rhistory[i];
542: }
543: }
544: } else if (0.0 == ctx->r_beta) {
545: if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
546: /* Compute summations for scalar scaling */
547: VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
549: VecDot(ctx->W, ctx->Y[0], &ys_sum);
550: VecDot(ctx->W, ctx->W, &ss_sum);
551: yy_sum += ctx->yy_rhistory[0];
552: } else {
553: VecCopy(ctx->Q, ctx->P);
554: VecReciprocal(ctx->Q);
556: /* Compute summations for scalar scaling */
557: yy_sum = 0; /* No safeguard required */
558: ys_sum = 0; /* No safeguard required */
559: ss_sum = 0; /* No safeguard required */
560: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
561: VecPointwiseMult(ctx->W, ctx->S[i], ctx->Q);
562: VecDot(ctx->W, ctx->Y[i], &yDs);
563: ys_sum += yDs;
565: VecDot(ctx->W, ctx->W, &sDs);
566: ss_sum += sDs;
568: yy_sum += ctx->yy_rhistory[i];
569: }
570: }
571: } else if (1.0 == ctx->r_beta) {
572: /* Compute summations for scalar scaling */
573: yy_sum = 0; /* No safeguard required */
574: ys_sum = 0; /* No safeguard required */
575: ss_sum = 0; /* No safeguard required */
576: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
577: VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P);
578: VecDot(ctx->V, ctx->S[i], &yDs);
579: ys_sum += yDs;
581: VecDot(ctx->V, ctx->V, &yDy);
582: yy_sum += yDy;
584: ss_sum += ctx->ss_rhistory[i];
585: }
586: } else {
587: VecCopy(ctx->Q, ctx->P);
589: VecPow(ctx->P, ctx->r_beta);
590: VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q);
592: /* Compute summations for scalar scaling */
593: yy_sum = 0; /* No safeguard required */
594: ys_sum = 0; /* No safeguard required */
595: ss_sum = 0; /* No safeguard required */
596: for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
597: VecPointwiseMult(ctx->V, ctx->P, ctx->Y[i]);
598: VecPointwiseMult(ctx->W, ctx->Q, ctx->S[i]);
600: VecDot(ctx->V, ctx->V, &yDy);
601: VecDot(ctx->V, ctx->W, &yDs);
602: VecDot(ctx->W, ctx->W, &sDs);
604: yy_sum += yDy;
605: ys_sum += yDs;
606: ss_sum += sDs;
607: }
608: }
610: if (0.0 == ctx->r_alpha) {
611: /* Safeguard ys_sum */
612: if (0.0 == ys_sum) {
613: ys_sum = TAO_ZERO_SAFEGUARD;
614: }
616: sigmanew = ss_sum / ys_sum;
617: } else if (1.0 == ctx->r_alpha) {
618: /* Safeguard yy_sum */
619: if (0.0 == yy_sum) {
620: ys_sum = TAO_ZERO_SAFEGUARD;
621: }
623: sigmanew = ys_sum / yy_sum;
624: } else {
625: denom = 2*ctx->r_alpha*yy_sum;
627: /* Safeguard denom */
628: if (0.0 == denom) {
629: denom = TAO_ZERO_SAFEGUARD;
630: }
632: 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;
633: }
635: /* If Q has small values, then Q^(r_beta - 1) */
636: /* can have very large values. Hence, ys_sum */
637: /* and ss_sum can be infinity. In this case, */
638: /* sigmanew can either be not-a-number or infinity. */
640: if (PetscIsInfOrNanReal(sigmanew)) {
641: /* sigmanew is not-a-number; skip rescaling */
642: } else if (!sigmanew) {
643: /* sigmanew is zero; this is a bad case; skip rescaling */
644: } else {
645: /* sigmanew is positive */
646: VecScale(ctx->U, sigmanew);
647: }
648: break;
649: }
651: /* Modify for previous information */
652: switch(ctx->limitType) {
653: case MatLMVM_Limit_Average:
654: if (1.0 == ctx->mu) {
655: VecCopy(ctx->D, ctx->U);
656: } else if (ctx->mu) {
657: VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U);
658: }
659: break;
661: case MatLMVM_Limit_Relative:
662: if (ctx->mu) {
663: /* P = (1-mu) * D */
664: VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D);
665: /* Q = (1+mu) * D */
666: VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D);
667: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D);
668: }
669: break;
671: case MatLMVM_Limit_Absolute:
672: if (ctx->nu) {
673: VecCopy(ctx->P, ctx->D);
674: VecShift(ctx->P, -ctx->nu);
675: VecCopy(ctx->D, ctx->Q);
676: VecShift(ctx->Q, ctx->nu);
677: VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P);
678: }
679: break;
681: default:
682: VecCopy(ctx->U, ctx->D);
683: break;
684: }
685: break;
686: }
687: PetscObjectDereference((PetscObject)ctx->Xprev);
688: PetscObjectDereference((PetscObject)ctx->Gprev);
689: ctx->Xprev = ctx->S[ctx->lm];
690: ctx->Gprev = ctx->Y[ctx->lm];
691: PetscObjectReference((PetscObject)ctx->S[ctx->lm]);
692: PetscObjectReference((PetscObject)ctx->Y[ctx->lm]);
694: } else {
695: ++ctx->nrejects;
696: }
697: }
699: ++ctx->iter;
700: VecCopy(x, ctx->Xprev);
701: VecCopy(g, ctx->Gprev);
702: return(0);
703: }
707: extern PetscErrorCode MatLMVMSetDelta(Mat m, PetscReal d)
708: {
709: MatLMVMCtx *ctx;
711: PetscBool same;
715: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
716: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
717: MatShellGetContext(m,(void**)&ctx);
718: ctx->delta = PetscAbsReal(d);
719: ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
720: ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
721: return(0);
722: }
726: extern PetscErrorCode MatLMVMSetScale(Mat m, Vec s)
727: {
728: MatLMVMCtx *ctx;
730: PetscBool same;
734: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
735: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
736: MatShellGetContext(m,(void**)&ctx);
738: VecDestroy(&ctx->scale);
739: if (s) {
740: VecDuplicate(s,&ctx->scale);
741: }
742: return(0);
743: }
747: extern PetscErrorCode MatLMVMGetRejects(Mat m, PetscInt *nrejects)
748: {
749: MatLMVMCtx *ctx;
751: PetscBool same;
755: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
756: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
757: MatShellGetContext(m,(void**)&ctx);
758: *nrejects = ctx->nrejects;
759: return(0);
760: }
764: extern PetscErrorCode MatLMVMSetH0(Mat m, Mat A)
765: {
767: return(0);
768: }
772: extern PetscErrorCode MatLMVMGetX0(Mat m, Vec x)
773: {
775: return(0);
776: }
780: extern PetscErrorCode MatLMVMSetPrev(Mat M, Vec x, Vec g)
781: {
782: MatLMVMCtx *ctx;
784: PetscBool same;
789: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
790: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
791: MatShellGetContext(M,(void**)&ctx);
792: if (ctx->nupdates == 0) {
793: MatLMVMUpdate(M,x,g);
794: } else {
795: VecCopy(x,ctx->Xprev);
796: VecCopy(g,ctx->Gprev);
797: /* TODO scaling specific terms */
798: }
799: return(0);
800: }
804: extern PetscErrorCode MatLMVMRefine(Mat coarse, Mat op, Vec fineX, Vec fineG)
805: {
807: PetscBool same;
814: PetscObjectTypeCompare((PetscObject)coarse,MATSHELL,&same);
815: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
816: PetscObjectTypeCompare((PetscObject)op,MATSHELL,&same);
817: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
818: return(0);
819: }
823: extern PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v)
824: {
826: MatLMVMCtx *ctx;
827: PetscBool same;
832: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
833: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
834: MatShellGetContext(m,(void**)&ctx);
836: /* Perform allocations */
837: VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
838: VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
839: VecDuplicate(v,&ctx->D);
840: VecDuplicate(v,&ctx->U);
841: VecDuplicate(v,&ctx->V);
842: VecDuplicate(v,&ctx->W);
843: VecDuplicate(v,&ctx->P);
844: VecDuplicate(v,&ctx->Q);
845: ctx->allocated = PETSC_TRUE;
846: return(0);
847: }