Actual source code: diagbrdn.c
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;
10: VecCheckSameSize(F, 2, dX, 3);
11: VecCheckMatCompatible(B, dX, 3, F, 2);
12: VecPointwiseMult(dX, ldb->invD, F);
13: return 0;
14: }
16: /*------------------------------------------------------------*/
18: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
19: {
20: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
21: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
23: VecCheckSameSize(X, 2, Z, 3);
24: VecCheckMatCompatible(B, X, 2, Z, 3);
25: VecPointwiseDivide(Z, X, ldb->invD);
26: return 0;
27: }
29: /*------------------------------------------------------------*/
31: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
32: {
33: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
34: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
35: PetscInt old_k, i, start;
36: PetscScalar yty, ststmp, curvature, ytDy, stDs, ytDs;
37: PetscReal curvtol, sigma, yy_sum, ss_sum, ys_sum, denom;
39: if (!lmvm->m) return 0;
40: if (lmvm->prev_set) {
41: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
42: VecAYPX(lmvm->Xprev, -1.0, X);
43: VecAYPX(lmvm->Fprev, -1.0, F);
44: /* Compute tolerance for accepting the update */
45: VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
46: VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
47: VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
48: VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
49: if (PetscRealPart(ststmp) < lmvm->eps) {
50: curvtol = 0.0;
51: } else {
52: curvtol = lmvm->eps * PetscRealPart(ststmp);
53: }
54: /* Test the curvature for the update */
55: if (PetscRealPart(curvature) > curvtol) {
56: /* Update is good so we accept it */
57: old_k = lmvm->k;
58: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
59: /* If we hit the memory limit, shift the yty and yts arrays */
60: if (old_k == lmvm->k) {
61: for (i = 0; i <= lmvm->k - 1; ++i) {
62: ldb->yty[i] = ldb->yty[i + 1];
63: ldb->yts[i] = ldb->yts[i + 1];
64: ldb->sts[i] = ldb->sts[i + 1];
65: }
66: }
67: /* Accept dot products into the history */
68: VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
69: ldb->yty[lmvm->k] = PetscRealPart(yty);
70: ldb->yts[lmvm->k] = PetscRealPart(curvature);
71: ldb->sts[lmvm->k] = PetscRealPart(ststmp);
72: if (ldb->forward) {
73: /* We are doing diagonal scaling of the forward Hessian B */
74: /* BFGS = DFP = inv(D); */
75: VecCopy(ldb->invD, ldb->invDnew);
76: VecReciprocal(ldb->invDnew);
78: /* V = y*y */
79: VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]);
81: /* W = inv(D)*s */
82: VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]);
83: VecDot(ldb->W, lmvm->S[lmvm->k], &stDs);
85: /* Safeguard stDs */
86: stDs = PetscMax(PetscRealPart(stDs), ldb->tol);
88: if (1.0 != ldb->theta) {
89: /* BFGS portion of the update */
90: /* U = (inv(D)*s)*(inv(D)*s) */
91: VecPointwiseMult(ldb->U, ldb->W, ldb->W);
93: /* Assemble */
94: VecAXPBY(ldb->BFGS, -1.0 / stDs, 0.0, ldb->U);
95: }
96: if (0.0 != ldb->theta) {
97: /* DFP portion of the update */
98: /* U = inv(D)*s*y */
99: VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]);
101: /* Assemble */
102: VecAXPBY(ldb->DFP, stDs / ldb->yts[lmvm->k], 0.0, ldb->V);
103: VecAXPY(ldb->DFP, -2.0, ldb->U);
104: }
106: if (0.0 == ldb->theta) {
107: VecAXPY(ldb->invDnew, 1.0, ldb->BFGS);
108: } else if (1.0 == ldb->theta) {
109: VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->DFP);
110: } else {
111: /* Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
112: VecAXPBYPCZ(ldb->invDnew, 1.0 - ldb->theta, (ldb->theta) / ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP);
113: }
115: VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V);
116: /* Obtain inverse and ensure positive definite */
117: VecReciprocal(ldb->invDnew);
118: VecAbs(ldb->invDnew);
120: } else {
121: /* Inverse Hessian update instead. */
122: VecCopy(ldb->invD, ldb->invDnew);
124: /* V = s*s */
125: VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]);
127: /* W = D*y */
128: VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]);
129: VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy);
131: /* Safeguard ytDy */
132: ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);
134: if (1.0 != ldb->theta) {
135: /* BFGS portion of the update */
136: /* U = s*Dy */
137: VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]);
139: /* Assemble */
140: VecAXPBY(ldb->BFGS, ytDy / ldb->yts[lmvm->k], 0.0, ldb->V);
141: VecAXPY(ldb->BFGS, -2.0, ldb->U);
142: }
143: if (0.0 != ldb->theta) {
144: /* DFP portion of the update */
146: /* U = (inv(D)*y)*(inv(D)*y) */
147: VecPointwiseMult(ldb->U, ldb->W, ldb->W);
149: /* Assemble */
150: VecAXPBY(ldb->DFP, -1.0 / ytDy, 0.0, ldb->U);
151: }
153: if (0.0 == ldb->theta) {
154: VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->BFGS);
155: } else if (1.0 == ldb->theta) {
156: VecAXPY(ldb->invDnew, 1.0, ldb->DFP);
157: } else {
158: /* Broyden update U=(1-theta)*P + theta*Q */
159: VecAXPBYPCZ(ldb->invDnew, (1.0 - ldb->theta) / ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP);
160: }
161: VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V);
162: /* Ensure positive definite */
163: VecAbs(ldb->invDnew);
164: }
165: if (ldb->sigma_hist > 0) {
166: /* Start with re-scaling on the newly computed diagonal */
167: if (0.5 == ldb->beta) {
168: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
169: VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew);
170: VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);
172: VecDotBegin(ldb->V, lmvm->Y[0], &ytDy);
173: VecDotBegin(ldb->W, lmvm->S[0], &stDs);
174: VecDotEnd(ldb->V, lmvm->Y[0], &ytDy);
175: VecDotEnd(ldb->W, lmvm->S[0], &stDs);
177: ss_sum = PetscRealPart(stDs);
178: yy_sum = PetscRealPart(ytDy);
179: ys_sum = ldb->yts[0];
180: } else {
181: VecCopy(ldb->invDnew, ldb->U);
182: VecReciprocal(ldb->U);
184: /* Compute summations for scalar scaling */
185: yy_sum = 0; /* No safeguard required */
186: ys_sum = 0; /* No safeguard required */
187: ss_sum = 0; /* No safeguard required */
188: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
189: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
190: VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U);
191: VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);
193: VecDotBegin(ldb->W, lmvm->S[i], &stDs);
194: VecDotBegin(ldb->V, lmvm->Y[i], &ytDy);
195: VecDotEnd(ldb->W, lmvm->S[i], &stDs);
196: VecDotEnd(ldb->V, lmvm->Y[i], &ytDy);
198: ss_sum += PetscRealPart(stDs);
199: ys_sum += ldb->yts[i];
200: yy_sum += PetscRealPart(ytDy);
201: }
202: }
203: } else if (0.0 == ldb->beta) {
204: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
205: /* Compute summations for scalar scaling */
206: VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);
208: VecDotBegin(ldb->W, lmvm->Y[0], &ytDs);
209: VecDotBegin(ldb->W, ldb->W, &stDs);
210: VecDotEnd(ldb->W, lmvm->Y[0], &ytDs);
211: VecDotEnd(ldb->W, ldb->W, &stDs);
213: ys_sum = PetscRealPart(ytDs);
214: ss_sum = PetscRealPart(stDs);
215: yy_sum = ldb->yty[0];
216: } else {
217: VecCopy(ldb->invDnew, ldb->U);
218: VecReciprocal(ldb->U);
220: /* Compute summations for scalar scaling */
221: yy_sum = 0; /* No safeguard required */
222: ys_sum = 0; /* No safeguard required */
223: ss_sum = 0; /* No safeguard required */
224: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
225: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
226: VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);
228: VecDotBegin(ldb->W, lmvm->Y[i], &ytDs);
229: VecDotBegin(ldb->W, ldb->W, &stDs);
230: VecDotEnd(ldb->W, lmvm->Y[i], &ytDs);
231: VecDotEnd(ldb->W, ldb->W, &stDs);
233: ss_sum += PetscRealPart(stDs);
234: ys_sum += PetscRealPart(ytDs);
235: yy_sum += ldb->yty[i];
236: }
237: }
238: } else if (1.0 == ldb->beta) {
239: /* Compute summations for scalar scaling */
240: yy_sum = 0; /* No safeguard required */
241: ys_sum = 0; /* No safeguard required */
242: ss_sum = 0; /* No safeguard required */
243: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
244: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
245: VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew);
247: VecDotBegin(ldb->V, lmvm->S[i], &ytDs);
248: VecDotBegin(ldb->V, ldb->V, &ytDy);
249: VecDotEnd(ldb->V, lmvm->S[i], &ytDs);
250: VecDotEnd(ldb->V, ldb->V, &ytDy);
252: yy_sum += PetscRealPart(ytDy);
253: ys_sum += PetscRealPart(ytDs);
254: ss_sum += ldb->sts[i];
255: }
256: } else {
257: VecCopy(ldb->invDnew, ldb->U);
258: VecPow(ldb->U, ldb->beta - 1);
260: /* Compute summations for scalar scaling */
261: yy_sum = 0; /* No safeguard required */
262: ys_sum = 0; /* No safeguard required */
263: ss_sum = 0; /* No safeguard required */
264: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
265: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
266: VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]);
267: VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]);
269: VecDotBegin(ldb->V, ldb->W, &ytDs);
270: VecDotBegin(ldb->V, ldb->V, &ytDy);
271: VecDotBegin(ldb->W, ldb->W, &stDs);
272: VecDotEnd(ldb->V, ldb->W, &ytDs);
273: VecDotEnd(ldb->V, ldb->V, &ytDy);
274: VecDotEnd(ldb->W, ldb->W, &stDs);
276: yy_sum += PetscRealPart(ytDy);
277: ys_sum += PetscRealPart(ytDs);
278: ss_sum += PetscRealPart(stDs);
279: }
280: }
282: if (0.0 == ldb->alpha) {
283: /* Safeguard ys_sum */
284: ys_sum = PetscMax(ldb->tol, ys_sum);
286: sigma = ss_sum / ys_sum;
287: } else if (1.0 == ldb->alpha) {
288: /* yy_sum is never 0; if it were, we'd be at the minimum */
289: sigma = ys_sum / yy_sum;
290: } else {
291: denom = 2.0 * ldb->alpha * yy_sum;
293: /* Safeguard denom */
294: denom = PetscMax(ldb->tol, denom);
296: 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;
297: }
298: } else {
299: sigma = 1.0;
300: }
301: /* If Q has small values, then Q^(r_beta - 1)
302: can have very large values. Hence, ys_sum
303: and ss_sum can be infinity. In this case,
304: sigma can either be not-a-number or infinity. */
306: if (PetscIsInfOrNanScalar(sigma)) {
307: /* sigma is not-a-number; skip rescaling */
308: } else if (0.0 == sigma) {
309: /* sigma is zero; this is a bad case; skip rescaling */
310: } else {
311: /* sigma is positive */
312: VecScale(ldb->invDnew, sigma);
313: }
315: /* Combine the old diagonal and the new diagonal using a convex limiter */
316: if (1.0 == ldb->rho) {
317: VecCopy(ldb->invDnew, ldb->invD);
318: } else if (ldb->rho) VecAXPBY(ldb->invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew);
319: } else {
320: MatLMVMReset(B, PETSC_FALSE);
321: }
322: /* End DiagBrdn update */
323: }
324: /* Save the solution and function to be used in the next update */
325: VecCopy(X, lmvm->Xprev);
326: VecCopy(F, lmvm->Fprev);
327: lmvm->prev_set = PETSC_TRUE;
328: return 0;
329: }
331: /*------------------------------------------------------------*/
333: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
334: {
335: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
336: Mat_DiagBrdn *bctx = (Mat_DiagBrdn *)bdata->ctx;
337: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
338: Mat_DiagBrdn *mctx = (Mat_DiagBrdn *)mdata->ctx;
339: PetscInt i;
341: mctx->theta = bctx->theta;
342: mctx->alpha = bctx->alpha;
343: mctx->beta = bctx->beta;
344: mctx->rho = bctx->rho;
345: mctx->delta = bctx->delta;
346: mctx->delta_min = bctx->delta_min;
347: mctx->delta_max = bctx->delta_max;
348: mctx->tol = bctx->tol;
349: mctx->sigma = bctx->sigma;
350: mctx->sigma_hist = bctx->sigma_hist;
351: mctx->forward = bctx->forward;
352: VecCopy(bctx->invD, mctx->invD);
353: for (i = 0; i <= bdata->k; ++i) {
354: mctx->yty[i] = bctx->yty[i];
355: mctx->yts[i] = bctx->yts[i];
356: mctx->sts[i] = bctx->sts[i];
357: }
358: return 0;
359: }
361: /*------------------------------------------------------------*/
363: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
364: {
365: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
366: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
367: PetscBool isascii;
369: PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii);
370: if (isascii) {
371: PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", ldb->sigma_hist);
372: PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho);
373: PetscViewerASCIIPrintf(pv, "Convex factor: theta=%g\n", (double)ldb->theta);
374: }
375: MatView_LMVM(B, pv);
376: return 0;
377: }
379: /*------------------------------------------------------------*/
381: static PetscErrorCode MatSetFromOptions_DiagBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
382: {
383: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
384: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
386: MatSetFromOptions_LMVM(B, PetscOptionsObject);
387: PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
388: PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL);
389: PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL);
390: PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL);
391: PetscOptionsReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL);
392: PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL);
393: PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL);
394: PetscOptionsInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", ldb->sigma_hist, &ldb->sigma_hist, NULL);
395: PetscOptionsHeadEnd();
400: return 0;
401: }
403: /*------------------------------------------------------------*/
405: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
406: {
407: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
408: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
410: VecSet(ldb->invD, ldb->delta);
411: if (destructive && ldb->allocated) {
412: PetscFree3(ldb->yty, ldb->yts, ldb->sts);
413: VecDestroy(&ldb->invDnew);
414: VecDestroy(&ldb->invD);
415: VecDestroy(&ldb->BFGS);
416: VecDestroy(&ldb->DFP);
417: VecDestroy(&ldb->U);
418: VecDestroy(&ldb->V);
419: VecDestroy(&ldb->W);
420: ldb->allocated = PETSC_FALSE;
421: }
422: MatReset_LMVM(B, destructive);
423: return 0;
424: }
426: /*------------------------------------------------------------*/
428: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
429: {
430: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
431: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
433: MatAllocate_LMVM(B, X, F);
434: if (!ldb->allocated) {
435: PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
436: VecDuplicate(lmvm->Xprev, &ldb->invDnew);
437: VecDuplicate(lmvm->Xprev, &ldb->invD);
438: VecDuplicate(lmvm->Xprev, &ldb->BFGS);
439: VecDuplicate(lmvm->Xprev, &ldb->DFP);
440: VecDuplicate(lmvm->Xprev, &ldb->U);
441: VecDuplicate(lmvm->Xprev, &ldb->V);
442: VecDuplicate(lmvm->Xprev, &ldb->W);
443: ldb->allocated = PETSC_TRUE;
444: }
445: return 0;
446: }
448: /*------------------------------------------------------------*/
450: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
451: {
452: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
453: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
455: if (ldb->allocated) {
456: PetscFree3(ldb->yty, ldb->yts, ldb->sts);
457: VecDestroy(&ldb->invDnew);
458: VecDestroy(&ldb->invD);
459: VecDestroy(&ldb->BFGS);
460: VecDestroy(&ldb->DFP);
461: VecDestroy(&ldb->U);
462: VecDestroy(&ldb->V);
463: VecDestroy(&ldb->W);
464: ldb->allocated = PETSC_FALSE;
465: }
466: PetscFree(lmvm->ctx);
467: MatDestroy_LMVM(B);
468: return 0;
469: }
471: /*------------------------------------------------------------*/
473: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
474: {
475: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
476: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
478: MatSetUp_LMVM(B);
479: if (!ldb->allocated) {
480: PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
481: VecDuplicate(lmvm->Xprev, &ldb->invDnew);
482: VecDuplicate(lmvm->Xprev, &ldb->invD);
483: VecDuplicate(lmvm->Xprev, &ldb->BFGS);
484: VecDuplicate(lmvm->Xprev, &ldb->DFP);
485: VecDuplicate(lmvm->Xprev, &ldb->U);
486: VecDuplicate(lmvm->Xprev, &ldb->V);
487: VecDuplicate(lmvm->Xprev, &ldb->W);
488: ldb->allocated = PETSC_TRUE;
489: }
490: return 0;
491: }
493: /*------------------------------------------------------------*/
495: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
496: {
497: Mat_LMVM *lmvm;
498: Mat_DiagBrdn *ldb;
500: MatCreate_LMVM(B);
501: PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN);
502: B->ops->setup = MatSetUp_DiagBrdn;
503: B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
504: B->ops->destroy = MatDestroy_DiagBrdn;
505: B->ops->solve = MatSolve_DiagBrdn;
506: B->ops->view = MatView_DiagBrdn;
508: lmvm = (Mat_LMVM *)B->data;
509: lmvm->square = PETSC_TRUE;
510: lmvm->m = 1;
511: lmvm->ops->allocate = MatAllocate_DiagBrdn;
512: lmvm->ops->reset = MatReset_DiagBrdn;
513: lmvm->ops->mult = MatMult_DiagBrdn;
514: lmvm->ops->update = MatUpdate_DiagBrdn;
515: lmvm->ops->copy = MatCopy_DiagBrdn;
517: PetscNew(&ldb);
518: lmvm->ctx = (void *)ldb;
519: ldb->theta = 0.0;
520: ldb->alpha = 1.0;
521: ldb->rho = 1.0;
522: ldb->forward = PETSC_TRUE;
523: ldb->beta = 0.5;
524: ldb->sigma = 1.0;
525: ldb->delta = 1.0;
526: ldb->delta_min = 1e-7;
527: ldb->delta_max = 100.0;
528: ldb->tol = 1e-8;
529: ldb->sigma_hist = 1;
530: ldb->allocated = PETSC_FALSE;
531: return 0;
532: }
534: /*------------------------------------------------------------*/
536: /*@
537: MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
538: for approximating Hessians. It consists of a convex combination of DFP and BFGS
539: diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
540: To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
541: We also ensure positive definiteness by taking the `VecAbs()` of the final vector.
543: There are two ways of approximating the diagonal: using the forward (B) update
544: schemes for BFGS and DFP and then taking the inverse, or directly working with
545: the inverse (H) update schemes for the BFGS and DFP updates, derived using the
546: Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a
547: parameter below.
549: In order to use the DiagBrdn matrix with other vector types, i.e. doing matrix-vector products
550: and matrix solves, the matrix must first be created using `MatCreate()` and `MatSetType()`,
551: followed by `MatLMVMAllocate()`. Then it will be available for updating
552: (via `MatLMVMUpdate()`) in one's favored solver implementation.
554: Collective
556: Input Parameters:
557: + comm - MPI communicator
558: . n - number of local rows for storage vectors
559: - N - global size of the storage vectors
561: Output Parameter:
562: . B - the matrix
564: Options Database Keys:
565: + -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
566: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
567: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
568: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
569: . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
570: . -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
571: - -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal
573: Level: intermediate
575: Note:
576: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
577: paradigm instead of this routine directly.
579: .seealso: [](chapter_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDIAGBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
580: `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMSymBrdn()`
581: @*/
582: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
583: {
584: MatCreate(comm, B);
585: MatSetSizes(*B, n, n, N, N);
586: MatSetType(*B, MATLMVMDIAGBROYDEN);
587: MatSetUp(*B);
588: return 0;
589: }