Actual source code: diagbrdn.c
petsc-3.13.6 2020-09-29
1: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
3: /*------------------------------------------------------------*/
5: static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
6: {
7: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
8: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
9: PetscErrorCode ierr;
12: VecCheckSameSize(F, 2, dX, 3);
13: VecCheckMatCompatible(B, dX, 3, F, 2);
14: VecPointwiseMult(dX, ldb->invD, F);
15: return(0);
16: }
18: /*------------------------------------------------------------*/
20: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
21: {
22: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
23: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
24: PetscErrorCode ierr;
27: VecCheckSameSize(X, 2, Z, 3);
28: VecCheckMatCompatible(B, X, 2, Z, 3);
29: VecPointwiseDivide(Z, X, ldb->invD);
30: return(0);
31: }
33: /*------------------------------------------------------------*/
35: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
36: {
37: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
38: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
39: PetscErrorCode ierr;
40: PetscInt old_k, i, start;
41: PetscScalar yty, ststmp, curvature, ytDy, stDs, ytDs;
42: PetscReal curvtol, sigma, yy_sum, ss_sum, ys_sum, denom;
45: if (!lmvm->m) return(0);
46: if (lmvm->prev_set) {
47: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
48: VecAYPX(lmvm->Xprev, -1.0, X);
49: VecAYPX(lmvm->Fprev, -1.0, F);
50: /* Compute tolerance for accepting the update */
51: VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
52: VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
53: VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
54: VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
55: if (PetscRealPart(ststmp) < lmvm->eps){
56: curvtol = 0.0;
57: } else {
58: curvtol = lmvm->eps * PetscRealPart(ststmp);
59: }
60: /* Test the curvature for the update */
61: if (PetscRealPart(curvature) > curvtol) {
62: /* Update is good so we accept it */
63: old_k = lmvm->k;
64: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
65: /* If we hit the memory limit, shift the yty and yts arrays */
66: if (old_k == lmvm->k) {
67: for (i = 0; i <= lmvm->k-1; ++i) {
68: ldb->yty[i] = ldb->yty[i+1];
69: ldb->yts[i] = ldb->yts[i+1];
70: ldb->sts[i] = ldb->sts[i+1];
71: }
72: }
73: /* Accept dot products into the history */
74: VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
75: ldb->yty[lmvm->k] = PetscRealPart(yty);
76: ldb->yts[lmvm->k] = PetscRealPart(curvature);
77: ldb->sts[lmvm->k] = PetscRealPart(ststmp);
78: if (ldb->forward) {
79: /* We are doing diagonal scaling of the forward Hessian B */
80: /* BFGS = DFP = inv(D); */
81: VecCopy(ldb->invD, ldb->invDnew);
82: VecReciprocal(ldb->invDnew);
84: /* V = y*y */
85: VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]);
87: /* W = inv(D)*s */
88: VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]);
89: VecDot(ldb->W, lmvm->S[lmvm->k], &stDs);
91: /* Safeguard stDs */
92: stDs = PetscMax(PetscRealPart(stDs), ldb->tol);
94: if (1.0 != ldb->theta) {
95: /* BFGS portion of the update */
96: /* U = (inv(D)*s)*(inv(D)*s) */
97: VecPointwiseMult(ldb->U, ldb->W, ldb->W);
99: /* Assemble */
100: VecAXPBY(ldb->BFGS, -1.0/stDs, 0.0, ldb->U);
101: }
102: if (0.0 != ldb->theta) {
103: /* DFP portion of the update */
104: /* U = inv(D)*s*y */
105: VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]);
107: /* Assemble */
108: VecAXPBY(ldb->DFP, stDs/ldb->yts[lmvm->k], 0.0, ldb->V);
109: VecAXPY(ldb->DFP, -2.0, ldb->U);
110: }
112: if (0.0 == ldb->theta) {
113: VecAXPY(ldb->invDnew, 1.0, ldb->BFGS);
114: } else if (1.0 == ldb->theta) {
115: VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->DFP);
116: } else {
117: /* Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
118: VecAXPBYPCZ(ldb->invDnew, 1.0-ldb->theta, (ldb->theta)/ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP);
119: }
121: VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
122: /* Obtain inverse and ensure positive definite */
123: VecReciprocal(ldb->invDnew);
124: VecAbs(ldb->invDnew);
126: } else {
127: /* Inverse Hessian update instead. */
128: VecCopy(ldb->invD, ldb->invDnew);
130: /* V = s*s */
131: VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]);
133: /* W = D*y */
134: VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]);
135: VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy);
137: /* Safeguard ytDy */
138: ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);
140: if (1.0 != ldb->theta) {
141: /* BFGS portion of the update */
142: /* U = s*Dy */
143: VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]);
145: /* Assemble */
146: VecAXPBY(ldb->BFGS, ytDy/ldb->yts[lmvm->k], 0.0, ldb->V);
147: VecAXPY(ldb->BFGS, -2.0, ldb->U);
148: }
149: if (0.0 != ldb->theta) {
150: /* DFP portion of the update */
152: /* U = (inv(D)*y)*(inv(D)*y) */
153: VecPointwiseMult(ldb->U, ldb->W, ldb->W);
155: /* Assemble */
156: VecAXPBY(ldb->DFP, -1.0/ytDy, 0.0, ldb->U);
157: }
159: if (0.0 == ldb->theta) {
160: VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->BFGS);
161: } else if (1.0 == ldb->theta) {
162: VecAXPY(ldb->invDnew, 1.0, ldb->DFP);
163: } else {
164: /* Broyden update U=(1-theta)*P + theta*Q */
165: VecAXPBYPCZ(ldb->invDnew, (1.0-ldb->theta)/ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP);
166: }
167: VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
168: /* Ensure positive definite */
169: VecAbs(ldb->invDnew);
170: }
171: if (ldb->sigma_hist > 0){
172: /* Start with re-scaling on the newly computed diagonal */
173: if (0.5 == ldb->beta) {
174: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
175: VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew);
176: VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);
178: VecDotBegin(ldb->V, lmvm->Y[0], &ytDy);
179: VecDotBegin(ldb->W, lmvm->S[0], &stDs);
180: VecDotEnd(ldb->V, lmvm->Y[0], &ytDy);
181: VecDotEnd(ldb->W, lmvm->S[0], &stDs);
183: ss_sum = PetscRealPart(stDs);
184: yy_sum = PetscRealPart(ytDy);
185: ys_sum = ldb->yts[0];
186: } else {
187: VecCopy(ldb->invDnew, ldb->U);
188: VecReciprocal(ldb->U);
190: /* Compute summations for scalar scaling */
191: yy_sum = 0; /* No safeguard required */
192: ys_sum = 0; /* No safeguard required */
193: ss_sum = 0; /* No safeguard required */
194: start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
195: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
196: VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U);
197: VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);
199: VecDotBegin(ldb->W, lmvm->S[i], &stDs);
200: VecDotBegin(ldb->V, lmvm->Y[i], &ytDy);
201: VecDotEnd(ldb->W, lmvm->S[i], &stDs);
202: VecDotEnd(ldb->V, lmvm->Y[i], &ytDy);
204: ss_sum += PetscRealPart(stDs);
205: ys_sum += ldb->yts[i];
206: yy_sum += PetscRealPart(ytDy);
207: }
208: }
209: } else if (0.0 == ldb->beta) {
210: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
211: /* Compute summations for scalar scaling */
212: VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);
214: VecDotBegin(ldb->W, lmvm->Y[0], &ytDs);
215: VecDotBegin(ldb->W, ldb->W, &stDs);
216: VecDotEnd(ldb->W, lmvm->Y[0], &ytDs);
217: VecDotEnd(ldb->W, ldb->W, &stDs);
219: ys_sum = PetscRealPart(ytDs);
220: ss_sum = PetscRealPart(stDs);
221: yy_sum = ldb->yty[0];
222: } else {
223: VecCopy(ldb->invDnew, ldb->U);
224: VecReciprocal(ldb->U);
226: /* Compute summations for scalar scaling */
227: yy_sum = 0; /* No safeguard required */
228: ys_sum = 0; /* No safeguard required */
229: ss_sum = 0; /* No safeguard required */
230: start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
231: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
232: VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);
234: VecDotBegin(ldb->W, lmvm->Y[i], &ytDs);
235: VecDotBegin(ldb->W, ldb->W, &stDs);
236: VecDotEnd(ldb->W, lmvm->Y[i], &ytDs);
237: VecDotEnd(ldb->W, ldb->W, &stDs);
239: ss_sum += PetscRealPart(stDs);
240: ys_sum += PetscRealPart(ytDs);
241: yy_sum += ldb->yty[i];
242: }
243: }
244: } else if (1.0 == ldb->beta) {
245: /* Compute summations for scalar scaling */
246: yy_sum = 0; /* No safeguard required */
247: ys_sum = 0; /* No safeguard required */
248: ss_sum = 0; /* No safeguard required */
249: start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
250: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
251: VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew);
253: VecDotBegin(ldb->V, lmvm->S[i], &ytDs);
254: VecDotBegin(ldb->V, ldb->V, &ytDy);
255: VecDotEnd(ldb->V, lmvm->S[i], &ytDs);
256: VecDotEnd(ldb->V, ldb->V, &ytDy);
258: yy_sum += PetscRealPart(ytDy);
259: ys_sum += PetscRealPart(ytDs);
260: ss_sum += ldb->sts[i];
261: }
262: } else {
263: VecCopy(ldb->invDnew, ldb->U);
264: VecPow(ldb->U, ldb->beta-1);
266: /* Compute summations for scalar scaling */
267: yy_sum = 0; /* No safeguard required */
268: ys_sum = 0; /* No safeguard required */
269: ss_sum = 0; /* No safeguard required */
270: start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
271: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
272: VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]);
273: VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]);
275: VecDotBegin(ldb->V, ldb->W, &ytDs);
276: VecDotBegin(ldb->V, ldb->V, &ytDy);
277: VecDotBegin(ldb->W, ldb->W, &stDs);
278: VecDotEnd(ldb->V, ldb->W, &ytDs);
279: VecDotEnd(ldb->V, ldb->V, &ytDy);
280: VecDotEnd(ldb->W, ldb->W, &stDs);
282: yy_sum += PetscRealPart(ytDy);
283: ys_sum += PetscRealPart(ytDs);
284: ss_sum += PetscRealPart(stDs);
285: }
286: }
288: if (0.0 == ldb->alpha) {
289: /* Safeguard ys_sum */
290: ys_sum = PetscMax(ldb->tol, ys_sum);
292: sigma = ss_sum / ys_sum;
293: } else if (1.0 == ldb->alpha) {
294: /* yy_sum is never 0; if it were, we'd be at the minimum */
295: sigma = ys_sum / yy_sum;
296: } else {
297: denom = 2.0*ldb->alpha*yy_sum;
299: /* Safeguard denom */
300: denom = PetscMax(ldb->tol, denom);
302: sigma = ((2.0*ldb->alpha-1)*ys_sum + PetscSqrtReal((2.0*ldb->alpha-1)*(2.0*ldb->alpha-1)*ys_sum*ys_sum - 4.0*ldb->alpha*(ldb->alpha-1)*yy_sum*ss_sum)) / denom;
303: }
304: } else {
305: sigma = 1.0;
306: }
307: /* If Q has small values, then Q^(r_beta - 1)
308: can have very large values. Hence, ys_sum
309: and ss_sum can be infinity. In this case,
310: sigma can either be not-a-number or infinity. */
312: if (PetscIsInfOrNanScalar(sigma)) {
313: /* sigma is not-a-number; skip rescaling */
314: } else if (0.0 == sigma) {
315: /* sigma is zero; this is a bad case; skip rescaling */
316: } else {
317: /* sigma is positive */
318: VecScale(ldb->invDnew, sigma);
319: }
321: /* Combine the old diagonal and the new diagonal using a convex limiter */
322: if (1.0 == ldb->rho) {
323: VecCopy(ldb->invDnew, ldb->invD);
324: } else if (ldb->rho) {
325: VecAXPBY(ldb->invD, 1.0-ldb->rho, ldb->rho, ldb->invDnew);
326: }
327: } else {
328: MatLMVMReset(B, PETSC_FALSE);
329: }
330: /* End DiagBrdn update */
332: }
333: /* Save the solution and function to be used in the next update */
334: VecCopy(X, lmvm->Xprev);
335: VecCopy(F, lmvm->Fprev);
336: lmvm->prev_set = PETSC_TRUE;
337: return(0);
338: }
340: /*------------------------------------------------------------*/
342: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
343: {
344: Mat_LMVM *bdata = (Mat_LMVM*)B->data;
345: Mat_DiagBrdn *bctx = (Mat_DiagBrdn*)bdata->ctx;
346: Mat_LMVM *mdata = (Mat_LMVM*)M->data;
347: Mat_DiagBrdn *mctx = (Mat_DiagBrdn*)mdata->ctx;
348: PetscErrorCode ierr;
349: PetscInt i;
352: mctx->theta = bctx->theta;
353: mctx->alpha = bctx->alpha;
354: mctx->beta = bctx->beta;
355: mctx->rho = bctx->rho;
356: mctx->delta = bctx->delta;
357: mctx->delta_min = bctx->delta_min;
358: mctx->delta_max = bctx->delta_max;
359: mctx->tol = bctx->tol;
360: mctx->sigma = bctx->sigma;
361: mctx->sigma_hist = bctx->sigma_hist;
362: mctx->forward = bctx->forward;
363: VecCopy(bctx->invD, mctx->invD);
364: for (i=0; i<=bdata->k; ++i) {
365: mctx->yty[i] = bctx->yty[i];
366: mctx->yts[i] = bctx->yts[i];
367: mctx->sts[i] = bctx->sts[i];
368: }
369: return(0);
370: }
372: /*------------------------------------------------------------*/
374: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
375: {
376: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
377: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
378: PetscErrorCode ierr;
379: PetscBool isascii;
382: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
383: if (isascii) {
384: PetscViewerASCIIPrintf(pv,"Scale history: %d\n",ldb->sigma_hist);
385: PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)ldb->alpha, (double)ldb->beta, (double)ldb->rho);
386: PetscViewerASCIIPrintf(pv,"Convex factor: theta=%g\n", (double)ldb->theta);
387: }
388: MatView_LMVM(B, pv);
389: return(0);
390: }
392: /*------------------------------------------------------------*/
394: static PetscErrorCode MatSetFromOptions_DiagBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
395: {
396: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
397: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
398: PetscErrorCode ierr;
401: MatSetFromOptions_LMVM(PetscOptionsObject, B);
402: PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
403: PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",ldb->theta,&ldb->theta,NULL);
404: PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",ldb->rho,&ldb->rho,NULL);
405: PetscOptionsReal("-mat_lmvm_tol","(developer) tolerance for bounding rescaling denominator","",ldb->tol,&ldb->tol,NULL);
406: PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",ldb->alpha,&ldb->alpha,NULL);
407: PetscOptionsBool("-mat_lmvm_forward","Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.","",ldb->forward,&ldb->forward,NULL);
408: PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",ldb->beta,&ldb->beta,NULL);
409: PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",ldb->sigma_hist,&ldb->sigma_hist,NULL);
410: PetscOptionsTail();
411: if ((ldb->theta < 0.0) || (ldb->theta > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
412: if ((ldb->alpha < 0.0) || (ldb->alpha > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
413: if ((ldb->rho < 0.0) || (ldb->rho > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex update limiter in the J0 scaling cannot be outside the range of [0, 1]");
414: if (ldb->sigma_hist < 0) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
415: return(0);
416: }
418: /*------------------------------------------------------------*/
420: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
421: {
422: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
423: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
424: PetscErrorCode ierr;
427: VecSet(ldb->invD, ldb->delta);
428: if (destructive && ldb->allocated) {
429: PetscFree3(ldb->yty, ldb->yts, ldb->sts);
430: VecDestroy(&ldb->invDnew);
431: VecDestroy(&ldb->invD);
432: VecDestroy(&ldb->BFGS);
433: VecDestroy(&ldb->DFP);
434: VecDestroy(&ldb->U);
435: VecDestroy(&ldb->V);
436: VecDestroy(&ldb->W);
437: ldb->allocated = PETSC_FALSE;
438: }
439: MatReset_LMVM(B, destructive);
440: return(0);
441: }
443: /*------------------------------------------------------------*/
445: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
446: {
447: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
448: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
449: PetscErrorCode ierr;
450:
452: MatAllocate_LMVM(B, X, F);
453: if (!ldb->allocated) {
454: PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
455: VecDuplicate(lmvm->Xprev, &ldb->invDnew);
456: VecDuplicate(lmvm->Xprev, &ldb->invD);
457: VecDuplicate(lmvm->Xprev, &ldb->BFGS);
458: VecDuplicate(lmvm->Xprev, &ldb->DFP);
459: VecDuplicate(lmvm->Xprev, &ldb->U);
460: VecDuplicate(lmvm->Xprev, &ldb->V);
461: VecDuplicate(lmvm->Xprev, &ldb->W);
462: ldb->allocated = PETSC_TRUE;
463: }
464: return(0);
465: }
467: /*------------------------------------------------------------*/
469: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
470: {
471: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
472: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
473: PetscErrorCode ierr;
476: if (ldb->allocated) {
477: PetscFree3(ldb->yty, ldb->yts, ldb->sts);
478: VecDestroy(&ldb->invDnew);
479: VecDestroy(&ldb->invD);
480: VecDestroy(&ldb->BFGS);
481: VecDestroy(&ldb->DFP);
482: VecDestroy(&ldb->U);
483: VecDestroy(&ldb->V);
484: VecDestroy(&ldb->W);
485: ldb->allocated = PETSC_FALSE;
486: }
487: PetscFree(lmvm->ctx);
488: MatDestroy_LMVM(B);
489: return(0);
490: }
492: /*------------------------------------------------------------*/
494: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
495: {
496: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
497: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
498: PetscErrorCode ierr;
501: MatSetUp_LMVM(B);
502: if (!ldb->allocated) {
503: PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
504: VecDuplicate(lmvm->Xprev, &ldb->invDnew);
505: VecDuplicate(lmvm->Xprev, &ldb->invD);
506: VecDuplicate(lmvm->Xprev, &ldb->BFGS);
507: VecDuplicate(lmvm->Xprev, &ldb->DFP);
508: VecDuplicate(lmvm->Xprev, &ldb->U);
509: VecDuplicate(lmvm->Xprev, &ldb->V);
510: VecDuplicate(lmvm->Xprev, &ldb->W);
511: ldb->allocated = PETSC_TRUE;
512: }
513: return(0);
514: }
516: /*------------------------------------------------------------*/
518: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
519: {
520: Mat_LMVM *lmvm;
521: Mat_DiagBrdn *ldb;
522: PetscErrorCode ierr;
525: MatCreate_LMVM(B);
526: PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN);
527: B->ops->setup = MatSetUp_DiagBrdn;
528: B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
529: B->ops->destroy = MatDestroy_DiagBrdn;
530: B->ops->solve = MatSolve_DiagBrdn;
531: B->ops->view = MatView_DiagBrdn;
533: lmvm = (Mat_LMVM*)B->data;
534: lmvm->square = PETSC_TRUE;
535: lmvm->m = 1;
536: lmvm->ops->allocate = MatAllocate_DiagBrdn;
537: lmvm->ops->reset = MatReset_DiagBrdn;
538: lmvm->ops->mult = MatMult_DiagBrdn;
539: lmvm->ops->update = MatUpdate_DiagBrdn;
540: lmvm->ops->copy = MatCopy_DiagBrdn;
542: PetscNewLog(B, &ldb);
543: lmvm->ctx = (void*)ldb;
544: ldb->theta = 0.0;
545: ldb->alpha = 1.0;
546: ldb->rho = 1.0;
547: ldb->forward = PETSC_TRUE;
548: ldb->beta = 0.5;
549: ldb->sigma = 1.0;
550: ldb->delta = 1.0;
551: ldb->delta_min = 1e-7;
552: ldb->delta_max = 100.0;
553: ldb->tol = 1e-8;
554: ldb->sigma_hist = 1;
555: ldb->allocated = PETSC_FALSE;
556: return(0);
557: }
559: /*------------------------------------------------------------*/
561: /*@
562: MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
563: for approximating Hessians. It consists of a convex combination of DFP and BFGS
564: diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
565: To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
566: We also ensure positive definiteness by taking the VecAbs() of the final vector.
568: There are two ways of approximating the diagonal: using the forward (B) update
569: schemes for BFGS and DFP and then taking the inverse, or directly working with
570: the inverse (H) update schemes for the BFGS and DFP updates, derived using the
571: Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a
572: parameter below.
574: In order to use the DiagBrdn matrix with other vector types, i.e. doing MatMults
575: and MatSolves, the matrix must first be created using MatCreate() and MatSetType(),
576: followed by MatLMVMAllocate(). Then it will be available for updating
577: (via MatLMVMUpdate) in one's favored solver implementation.
578: This allows for MPI compatibility.
580: Collective
582: Input Parameters:
583: + comm - MPI communicator, set to PETSC_COMM_SELF
584: . n - number of local rows for storage vectors
585: - N - global size of the storage vectors
587: Output Parameter:
588: . B - the matrix
590: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
591: paradigm instead of this routine directly.
593: Options Database Keys:
594: + -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
595: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
596: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
597: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
598: . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
599: . -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
600: - -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal
602: Level: intermediate
604: .seealso: MatCreate(), MATLMVM, MATLMVMDIAGBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
605: MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMSymBrdn()
606: @*/
607: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
608: {
609: PetscErrorCode ierr;
612: MatCreate(comm, B);
613: MatSetSizes(*B, n, n, N, N);
614: MatSetType(*B, MATLMVMDIAGBROYDEN);
615: MatSetUp(*B);
616: return(0);
617: }