Actual source code: lmvmmat.c
petsc-3.5.4 2015-05-23
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: PetscOptionsInt("-tao_lmm_vectors", "vectors to use for approximation", "", ctx->lm, &ctx->lm, 0);
71: PetscOptionsReal("-tao_lmm_limit_mu", "mu limiting factor", "", ctx->mu, &ctx->mu, 0);
72: PetscOptionsReal("-tao_lmm_limit_nu", "nu limiting factor", "", ctx->nu, &ctx->nu, 0);
73: PetscOptionsReal("-tao_lmm_broyden_phi", "phi factor for Broyden scaling", "", ctx->phi, &ctx->phi, 0);
74: PetscOptionsReal("-tao_lmm_scalar_alpha", "alpha factor for scalar scaling", "",ctx->s_alpha, &ctx->s_alpha, 0);
75: PetscOptionsReal("-tao_lmm_rescale_alpha", "alpha factor for rescaling diagonal", "", ctx->r_alpha, &ctx->r_alpha, 0);
76: PetscOptionsReal("-tao_lmm_rescale_beta", "beta factor for rescaling diagonal", "", ctx->r_beta, &ctx->r_beta, 0);
77: PetscOptionsInt("-tao_lmm_scalar_history", "amount of history for scalar scaling", "", ctx->scalar_history, &ctx->scalar_history, 0);
78: PetscOptionsInt("-tao_lmm_rescale_history", "amount of history for rescaling diagonal", "", ctx->rescale_history, &ctx->rescale_history, 0);
79: PetscOptionsReal("-tao_lmm_eps", "rejection tolerance", "", ctx->eps, &ctx->eps, 0);
80: PetscOptionsEList("-tao_lmm_scale_type", "scale type", "", Scale_Table, MatLMVM_Scale_Types, Scale_Table[ctx->scaleType], &ctx->scaleType, 0);
81: PetscOptionsEList("-tao_lmm_rescale_type", "rescale type", "", Rescale_Table, MatLMVM_Rescale_Types, Rescale_Table[ctx->rScaleType], &ctx->rScaleType, 0);
82: PetscOptionsEList("-tao_lmm_limit_type", "limit type", "", Limit_Table, MatLMVM_Limit_Types, Limit_Table[ctx->limitType], &ctx->limitType, 0);
83: PetscOptionsReal("-tao_lmm_delta_min", "minimum delta value", "", ctx->delta_min, &ctx->delta_min, 0);
84: PetscOptionsReal("-tao_lmm_delta_max", "maximum delta value", "", ctx->delta_max, &ctx->delta_max, 0);
85:
86: /* Complete configuration */
87: ctx->rescale_history = PetscMin(ctx->rescale_history, ctx->lm);
88: PetscMalloc1(ctx->lm+1,&ctx->rho);
89: PetscMalloc1(ctx->lm+1,&ctx->beta);
91: nhistory = PetscMax(ctx->scalar_history,1);
92: PetscMalloc1(nhistory,&ctx->yy_history);
93: PetscMalloc1(nhistory,&ctx->ys_history);
94: PetscMalloc1(nhistory,&ctx->ss_history);
96: nhistory = PetscMax(ctx->rescale_history,1);
97: PetscMalloc1(nhistory,&ctx->yy_rhistory);
98: PetscMalloc1(nhistory,&ctx->ys_rhistory);
99: PetscMalloc1(nhistory,&ctx->ss_rhistory);
101: /* Finish initializations */
102: ctx->lmnow = 0;
103: ctx->iter = 0;
104: ctx->nupdates = 0;
105: ctx->nrejects = 0;
106: ctx->delta = 1.0;
108: ctx->Gprev = 0;
109: ctx->Xprev = 0;
111: ctx->scale = 0;
112: ctx->useScale = PETSC_FALSE;
114: ctx->H0 = 0;
115: ctx->useDefaultH0=PETSC_TRUE;
117: MatCreateShell(comm, n, n, N, N, ctx, A);
118: MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
119: MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
120: return(0);
121: }
125: extern PetscErrorCode MatLMVMSolve(Mat A, Vec b, Vec x)
126: {
127: PetscReal sq, yq, dd;
128: PetscInt ll;
129: PetscBool scaled;
130: MatLMVMCtx *shell;
137: MatShellGetContext(A,(void**)&shell);
138: if (shell->lmnow < 1) {
139: shell->rho[0] = 1.0;
140: }
142: VecCopy(b,x);
143: for (ll = 0; ll < shell->lmnow; ++ll) {
144: VecDot(x,shell->S[ll],&sq);
145: shell->beta[ll] = sq * shell->rho[ll];
146: VecAXPY(x,-shell->beta[ll],shell->Y[ll]);
147: }
149: scaled = PETSC_FALSE;
150: if (!scaled && !shell->useDefaultH0 && shell->H0) {
151: MatSolve(shell->H0,x,shell->U);
152: VecDot(x,shell->U,&dd);
153: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
154: /* Accept Hessian solve */
155: VecCopy(shell->U,x);
156: scaled = PETSC_TRUE;
157: }
158: }
160: if (!scaled && shell->useScale) {
161: VecPointwiseMult(shell->U,x,shell->scale);
162: VecDot(x,shell->U,&dd);
163: if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
164: /* Accept scaling */
165: VecCopy(shell->U,x);
166: scaled = PETSC_TRUE;
167: }
168: }
170: if (!scaled) {
171: switch(shell->scaleType) {
172: case MatLMVM_Scale_None:
173: break;
175: case MatLMVM_Scale_Scalar:
176: VecScale(x,shell->sigma);
177: break;
179: case MatLMVM_Scale_Broyden:
180: VecPointwiseMult(x,x,shell->D);
181: break;
182: }
183: }
184: for (ll = shell->lmnow-1; ll >= 0; --ll) {
185: VecDot(x,shell->Y[ll],&yq);
186: VecAXPY(x,shell->beta[ll]-yq*shell->rho[ll],shell->S[ll]);
187: }
188: return(0);
189: }
193: extern PetscErrorCode MatView_LMVM(Mat A, PetscViewer pv)
194: {
195: PetscBool isascii;
197: MatLMVMCtx *lmP;
200: MatShellGetContext(A,(void**)&lmP);
201: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
202: if (isascii) {
203: PetscViewerASCIIPrintf(pv,"LMVM Matrix\n");
204: PetscViewerASCIIPrintf(pv," Number of vectors: %D\n",lmP->lm);
205: PetscViewerASCIIPrintf(pv," scale type: %s\n",Scale_Table[lmP->scaleType]);
206: PetscViewerASCIIPrintf(pv," rescale type: %s\n",Rescale_Table[lmP->rScaleType]);
207: PetscViewerASCIIPrintf(pv," limit type: %s\n",Limit_Table[lmP->limitType]);
208: PetscViewerASCIIPrintf(pv," updates: %D\n",lmP->nupdates);
209: PetscViewerASCIIPrintf(pv," rejects: %D\n",lmP->nrejects);
210: }
211: return(0);
212: }
216: extern PetscErrorCode MatDestroy_LMVM(Mat M)
217: {
218: MatLMVMCtx *ctx;
222: MatShellGetContext(M,(void**)&ctx);
223: if (ctx->allocated) {
224: if (ctx->Xprev) {
225: PetscObjectDereference((PetscObject)ctx->Xprev);
226: }
227: if (ctx->Gprev) {
228: PetscObjectDereference((PetscObject)ctx->Gprev);
229: }
231: VecDestroyVecs(ctx->lm+1,&ctx->S);
232: VecDestroyVecs(ctx->lm+1,&ctx->Y);
233: VecDestroy(&ctx->D);
234: VecDestroy(&ctx->U);
235: VecDestroy(&ctx->V);
236: VecDestroy(&ctx->W);
237: VecDestroy(&ctx->P);
238: VecDestroy(&ctx->Q);
239: if (ctx->scale) {
240: VecDestroy(&ctx->scale);
241: }
242: }
243: PetscFree(ctx->rho);
244: PetscFree(ctx->beta);
245: PetscFree(ctx->yy_history);
246: PetscFree(ctx->ys_history);
247: PetscFree(ctx->ss_history);
248: PetscFree(ctx->yy_rhistory);
249: PetscFree(ctx->ys_rhistory);
250: PetscFree(ctx->ss_rhistory);
251: PetscFree(ctx);
252: return(0);
253: }
257: extern PetscErrorCode MatLMVMReset(Mat M)
258: {
260: MatLMVMCtx *ctx;
261: PetscInt i;
264: MatShellGetContext(M,(void**)&ctx);
265: if (ctx->Gprev) {
266: PetscObjectDereference((PetscObject)ctx->Gprev);
267: }
268: if (ctx->Xprev) {
269: PetscObjectDereference((PetscObject)ctx->Xprev);
270: }
271: ctx->Gprev = ctx->Y[ctx->lm];
272: ctx->Xprev = ctx->S[ctx->lm];
273: PetscObjectReference((PetscObject)ctx->Gprev);
274: PetscObjectReference((PetscObject)ctx->Xprev);
275: for (i=0; i<ctx->lm; ++i) {
276: ctx->rho[i] = 0.0;
277: }
278: ctx->rho[0] = 1.0;
280: /* Set the scaling and diagonal scaling matrix */
281: switch(ctx->scaleType) {
282: case MatLMVM_Scale_None:
283: ctx->sigma = 1.0;
284: break;
285: case MatLMVM_Scale_Scalar:
286: ctx->sigma = ctx->delta;
287: break;
288: case MatLMVM_Scale_Broyden:
289: VecSet(ctx->D,ctx->delta);
290: break;
291: }
293: ctx->iter=0;
294: ctx->nupdates=0;
295: ctx->lmnow=0;
296: return(0);
297: }
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: }
705: extern PetscErrorCode MatLMVMSetDelta(Mat m, PetscReal d)
706: {
707: MatLMVMCtx *ctx;
709: PetscBool same;
713: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
714: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
715: MatShellGetContext(m,(void**)&ctx);
716: ctx->delta = PetscAbsReal(d);
717: ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
718: ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
719: return(0);
720: }
724: extern PetscErrorCode MatLMVMSetScale(Mat m, Vec s)
725: {
726: MatLMVMCtx *ctx;
728: PetscBool same;
732: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
733: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
734: MatShellGetContext(m,(void**)&ctx);
736: VecDestroy(&ctx->scale);
737: if (s) {
738: VecDuplicate(s,&ctx->scale);
739: }
740: return(0);
741: }
745: extern PetscErrorCode MatLMVMGetRejects(Mat m, PetscInt *nrejects)
746: {
747: MatLMVMCtx *ctx;
749: PetscBool same;
753: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
754: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
755: MatShellGetContext(m,(void**)&ctx);
756: *nrejects = ctx->nrejects;
757: return(0);
758: }
762: extern PetscErrorCode MatLMVMSetH0(Mat m, Mat A)
763: {
765: return(0);
766: }
770: extern PetscErrorCode MatLMVMGetX0(Mat m, Vec x)
771: {
773: return(0);
774: }
778: extern PetscErrorCode MatLMVMSetPrev(Mat M, Vec x, Vec g)
779: {
780: MatLMVMCtx *ctx;
782: PetscBool same;
787: PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
788: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
789: MatShellGetContext(M,(void**)&ctx);
790: if (ctx->nupdates == 0) {
791: MatLMVMUpdate(M,x,g);
792: } else {
793: VecCopy(x,ctx->Xprev);
794: VecCopy(g,ctx->Gprev);
795: /* TODO scaling specific terms */
796: }
797: return(0);
798: }
802: extern PetscErrorCode MatLMVMRefine(Mat coarse, Mat op, Vec fineX, Vec fineG)
803: {
805: PetscBool same;
812: PetscObjectTypeCompare((PetscObject)coarse,MATSHELL,&same);
813: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
814: PetscObjectTypeCompare((PetscObject)op,MATSHELL,&same);
815: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
816: return(0);
817: }
821: extern PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v)
822: {
824: MatLMVMCtx *ctx;
825: PetscBool same;
830: PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
831: if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
832: MatShellGetContext(m,(void**)&ctx);
834: /* Perform allocations */
835: VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
836: VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
837: VecDuplicate(v,&ctx->D);
838: VecDuplicate(v,&ctx->U);
839: VecDuplicate(v,&ctx->V);
840: VecDuplicate(v,&ctx->W);
841: VecDuplicate(v,&ctx->P);
842: VecDuplicate(v,&ctx->Q);
843: ctx->allocated = PETSC_TRUE;
844: return(0);
845: }