Actual source code: diagbrdn.c
1: #include <petscdevice.h>
2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
4: static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
5: {
6: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
7: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
9: PetscFunctionBegin;
10: PetscCall(VecPointwiseMult(dX, ldb->invD, F));
11: PetscFunctionReturn(PETSC_SUCCESS);
12: }
14: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
15: {
16: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
17: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
19: PetscFunctionBegin;
20: PetscCall(VecPointwiseDivide(Z, X, ldb->invD));
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
25: {
26: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
27: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
28: PetscInt old_k, i, start;
29: PetscScalar curvature, ytDy, sts, stDs, ytDs;
30: PetscReal curvtol, sigma, yy_sum, ss_sum, ys_sum, denom, ytytmp;
31: PetscReal stDsr, ytDyr;
33: PetscFunctionBegin;
34: if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
35: if (lmvm->prev_set) {
36: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
37: PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
38: PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
40: /* Test if the updates can be accepted */
41: PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ytytmp));
42: if (ytytmp < lmvm->eps) curvtol = 0.0;
43: else curvtol = lmvm->eps * ytytmp;
45: /* Test the curvature for the update */
46: if (PetscRealPart(curvature) > curvtol) {
47: /* Update is good so we accept it */
48: old_k = lmvm->k;
49: PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
50: /* If we hit the memory limit, shift the yty and yts arrays */
51: if (old_k == lmvm->k) {
52: for (i = 0; i <= lmvm->k - 1; ++i) {
53: ldb->yty[i] = ldb->yty[i + 1];
54: ldb->yts[i] = ldb->yts[i + 1];
55: ldb->sts[i] = ldb->sts[i + 1];
56: }
57: }
58: /* Accept dot products into the history */
59: PetscCall(VecDot(lmvm->S[lmvm->k], lmvm->S[lmvm->k], &sts));
60: ldb->yty[lmvm->k] = ytytmp;
61: ldb->yts[lmvm->k] = PetscRealPart(curvature);
62: ldb->sts[lmvm->k] = PetscRealPart(sts);
63: if (ldb->forward) {
64: /* We are doing diagonal scaling of the forward Hessian B */
65: /* BFGS = DFP = inv(D); */
66: PetscCall(VecCopy(ldb->invD, ldb->invDnew));
67: PetscCall(VecReciprocal(ldb->invDnew));
69: /* V = y*y */
70: PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]));
72: /* W = inv(D)*s */
73: PetscCall(VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]));
74: PetscCall(VecDot(ldb->W, lmvm->S[lmvm->k], &stDs));
76: /* Safeguard stDs */
77: stDs = PetscMax(PetscRealPart(stDs), ldb->tol);
79: if (1.0 != ldb->theta) {
80: /* BFGS portion of the update */
81: /* U = (inv(D)*s)*(inv(D)*s) */
82: PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));
84: /* Assemble */
85: PetscCall(VecAXPBY(ldb->BFGS, -1.0 / stDs, 0.0, ldb->U));
86: }
87: if (0.0 != ldb->theta) {
88: /* DFP portion of the update */
89: /* U = inv(D)*s*y */
90: PetscCall(VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]));
92: /* Assemble */
93: PetscCall(VecAXPBY(ldb->DFP, stDs / ldb->yts[lmvm->k], 0.0, ldb->V));
94: PetscCall(VecAXPY(ldb->DFP, -2.0, ldb->U));
95: }
97: if (0.0 == ldb->theta) {
98: PetscCall(VecAXPY(ldb->invDnew, 1.0, ldb->BFGS));
99: } else if (1.0 == ldb->theta) {
100: PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->DFP));
101: } else {
102: /* Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
103: PetscCall(VecAXPBYPCZ(ldb->invDnew, 1.0 - ldb->theta, (ldb->theta) / ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP));
104: }
106: PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
107: /* Obtain inverse and ensure positive definite */
108: PetscCall(VecReciprocal(ldb->invDnew));
109: PetscCall(VecAbs(ldb->invDnew));
111: } else {
112: /* Inverse Hessian update instead. */
113: PetscCall(VecCopy(ldb->invD, ldb->invDnew));
115: /* V = s*s */
116: PetscCall(VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]));
118: /* W = D*y */
119: PetscCall(VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]));
120: PetscCall(VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy));
122: /* Safeguard ytDy */
123: ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);
125: if (1.0 != ldb->theta) {
126: /* BFGS portion of the update */
127: /* U = s*Dy */
128: PetscCall(VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]));
130: /* Assemble */
131: PetscCall(VecAXPBY(ldb->BFGS, ytDy / ldb->yts[lmvm->k], 0.0, ldb->V));
132: PetscCall(VecAXPY(ldb->BFGS, -2.0, ldb->U));
133: }
134: if (0.0 != ldb->theta) {
135: /* DFP portion of the update */
137: /* U = (inv(D)*y)*(inv(D)*y) */
138: PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));
140: /* Assemble */
141: PetscCall(VecAXPBY(ldb->DFP, -1.0 / ytDy, 0.0, ldb->U));
142: }
144: if (0.0 == ldb->theta) {
145: PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->BFGS));
146: } else if (1.0 == ldb->theta) {
147: PetscCall(VecAXPY(ldb->invDnew, 1.0, ldb->DFP));
148: } else {
149: /* Broyden update U=(1-theta)*P + theta*Q */
150: PetscCall(VecAXPBYPCZ(ldb->invDnew, (1.0 - ldb->theta) / ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP));
151: }
152: PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
153: /* Ensure positive definite */
154: PetscCall(VecAbs(ldb->invDnew));
155: }
156: if (ldb->sigma_hist > 0) {
157: /* Start with re-scaling on the newly computed diagonal */
158: if (0.5 == ldb->beta) {
159: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
160: PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew));
161: PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));
163: PetscCall(VecDot(ldb->V, lmvm->Y[0], &ytDy));
164: PetscCall(VecDot(ldb->W, lmvm->S[0], &stDs));
166: ss_sum = PetscRealPart(stDs);
167: yy_sum = PetscRealPart(ytDy);
168: ys_sum = ldb->yts[0];
169: } else {
170: PetscCall(VecCopy(ldb->invDnew, ldb->U));
171: PetscCall(VecReciprocal(ldb->U));
173: /* Compute summations for scalar scaling */
174: yy_sum = 0; /* No safeguard required */
175: ys_sum = 0; /* No safeguard required */
176: ss_sum = 0; /* No safeguard required */
177: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
178: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
179: PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U));
180: PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));
182: PetscCall(VecDot(ldb->W, lmvm->S[i], &stDs));
183: PetscCall(VecDot(ldb->V, lmvm->Y[i], &ytDy));
185: ss_sum += PetscRealPart(stDs);
186: ys_sum += ldb->yts[i];
187: yy_sum += PetscRealPart(ytDy);
188: }
189: }
190: } else if (0.0 == ldb->beta) {
191: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
192: /* Compute summations for scalar scaling */
193: PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));
195: PetscCall(VecDotNorm2(lmvm->Y[0], ldb->W, &ytDs, &stDsr));
197: ys_sum = PetscRealPart(ytDs);
198: ss_sum = stDsr;
199: yy_sum = ldb->yty[0];
200: } else {
201: PetscCall(VecCopy(ldb->invDnew, ldb->U));
202: PetscCall(VecReciprocal(ldb->U));
204: /* Compute summations for scalar scaling */
205: yy_sum = 0; /* No safeguard required */
206: ys_sum = 0; /* No safeguard required */
207: ss_sum = 0; /* No safeguard required */
208: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
209: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
210: PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));
212: PetscCall(VecDotNorm2(lmvm->Y[i], ldb->W, &ytDs, &stDsr));
214: ss_sum += stDsr;
215: ys_sum += PetscRealPart(ytDs);
216: yy_sum += ldb->yty[i];
217: }
218: }
219: } else if (1.0 == ldb->beta) {
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: PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew));
228: PetscCall(VecDotNorm2(lmvm->S[i], ldb->V, &ytDs, &ytDyr));
230: yy_sum += ytDyr;
231: ys_sum += PetscRealPart(ytDs);
232: ss_sum += ldb->sts[i];
233: }
234: } else {
235: PetscCall(VecCopy(ldb->invDnew, ldb->U));
236: PetscCall(VecPow(ldb->U, ldb->beta - 1));
238: /* Compute summations for scalar scaling */
239: yy_sum = 0; /* No safeguard required */
240: ys_sum = 0; /* No safeguard required */
241: ss_sum = 0; /* No safeguard required */
242: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
243: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
244: PetscCall(VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]));
245: PetscCall(VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]));
247: PetscCall(VecDotNorm2(ldb->W, ldb->V, &ytDs, &ytDyr));
248: PetscCall(VecDot(ldb->W, ldb->W, &stDs));
250: yy_sum += ytDyr;
251: ys_sum += PetscRealPart(ytDs);
252: ss_sum += PetscRealPart(stDs);
253: }
254: }
256: if (0.0 == ldb->alpha) {
257: /* Safeguard ys_sum */
258: ys_sum = PetscMax(ldb->tol, ys_sum);
260: sigma = ss_sum / ys_sum;
261: } else if (1.0 == ldb->alpha) {
262: /* yy_sum is never 0; if it were, we'd be at the minimum */
263: sigma = ys_sum / yy_sum;
264: } else {
265: denom = 2.0 * ldb->alpha * yy_sum;
267: /* Safeguard denom */
268: denom = PetscMax(ldb->tol, denom);
270: 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;
271: }
272: } else {
273: sigma = 1.0;
274: }
275: /* If Q has small values, then Q^(r_beta - 1)
276: can have very large values. Hence, ys_sum
277: and ss_sum can be infinity. In this case,
278: sigma can either be not-a-number or infinity. */
280: if (PetscIsInfOrNanScalar(sigma)) {
281: /* sigma is not-a-number; skip rescaling */
282: } else if (0.0 == sigma) {
283: /* sigma is zero; this is a bad case; skip rescaling */
284: } else {
285: /* sigma is positive */
286: PetscCall(VecScale(ldb->invDnew, sigma));
287: }
289: /* Combine the old diagonal and the new diagonal using a convex limiter */
290: if (1.0 == ldb->rho) {
291: PetscCall(VecCopy(ldb->invDnew, ldb->invD));
292: } else if (ldb->rho) PetscCall(VecAXPBY(ldb->invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew));
293: } else {
294: PetscCall(MatLMVMReset(B, PETSC_FALSE));
295: }
296: /* End DiagBrdn update */
297: }
298: /* Save the solution and function to be used in the next update */
299: PetscCall(VecCopy(X, lmvm->Xprev));
300: PetscCall(VecCopy(F, lmvm->Fprev));
301: lmvm->prev_set = PETSC_TRUE;
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
306: {
307: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
308: Mat_DiagBrdn *bctx = (Mat_DiagBrdn *)bdata->ctx;
309: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
310: Mat_DiagBrdn *mctx = (Mat_DiagBrdn *)mdata->ctx;
311: PetscInt i;
313: PetscFunctionBegin;
314: mctx->theta = bctx->theta;
315: mctx->alpha = bctx->alpha;
316: mctx->beta = bctx->beta;
317: mctx->rho = bctx->rho;
318: mctx->delta = bctx->delta;
319: mctx->delta_min = bctx->delta_min;
320: mctx->delta_max = bctx->delta_max;
321: mctx->tol = bctx->tol;
322: mctx->sigma = bctx->sigma;
323: mctx->sigma_hist = bctx->sigma_hist;
324: mctx->forward = bctx->forward;
325: PetscCall(VecCopy(bctx->invD, mctx->invD));
326: for (i = 0; i <= bdata->k; ++i) {
327: mctx->yty[i] = bctx->yty[i];
328: mctx->yts[i] = bctx->yts[i];
329: mctx->sts[i] = bctx->sts[i];
330: }
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
335: {
336: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
337: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
338: PetscBool isascii;
340: PetscFunctionBegin;
341: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
342: if (isascii) {
343: PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", ldb->sigma_hist));
344: PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho));
345: PetscCall(PetscViewerASCIIPrintf(pv, "Convex factor: theta=%g\n", (double)ldb->theta));
346: }
347: PetscCall(MatView_LMVM(B, pv));
348: PetscFunctionReturn(PETSC_SUCCESS);
349: }
351: static PetscErrorCode MatSetFromOptions_DiagBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
352: {
353: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
354: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
356: PetscFunctionBegin;
357: PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
358: PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
359: PetscCall(PetscOptionsRangeReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL, 0.0, 1.0));
360: PetscCall(PetscOptionsRangeReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL, 0.0, 1.0));
361: PetscCall(PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL));
362: PetscCall(PetscOptionsRangeReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL, 0.0, 1.0));
363: PetscCall(PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL));
364: PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL));
365: PetscCall(PetscOptionsBoundedInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", ldb->sigma_hist, &ldb->sigma_hist, NULL, 0));
366: PetscOptionsHeadEnd();
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: static PetscErrorCode MatSetUp_DiagBrdn(Mat);
371: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
372: {
373: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
374: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
376: PetscFunctionBegin;
377: if (!ldb->allocated) PetscCall(MatSetUp_DiagBrdn(B));
378: PetscCall(VecSet(ldb->invD, ldb->delta));
379: if (destructive && ldb->allocated) {
380: PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
381: PetscCall(VecDestroy(&ldb->invDnew));
382: PetscCall(VecDestroy(&ldb->invD));
383: PetscCall(VecDestroy(&ldb->BFGS));
384: PetscCall(VecDestroy(&ldb->DFP));
385: PetscCall(VecDestroy(&ldb->U));
386: PetscCall(VecDestroy(&ldb->V));
387: PetscCall(VecDestroy(&ldb->W));
388: ldb->allocated = PETSC_FALSE;
389: }
390: PetscCall(MatReset_LMVM(B, destructive));
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
395: {
396: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
397: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
399: PetscFunctionBegin;
400: PetscCall(MatAllocate_LMVM(B, X, F));
401: if (!ldb->allocated) {
402: PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
403: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
404: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
405: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
406: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
407: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
408: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
409: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
410: ldb->allocated = PETSC_TRUE;
411: }
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
416: {
417: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
418: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
420: PetscFunctionBegin;
421: if (ldb->allocated) {
422: PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
423: PetscCall(VecDestroy(&ldb->invDnew));
424: PetscCall(VecDestroy(&ldb->invD));
425: PetscCall(VecDestroy(&ldb->BFGS));
426: PetscCall(VecDestroy(&ldb->DFP));
427: PetscCall(VecDestroy(&ldb->U));
428: PetscCall(VecDestroy(&ldb->V));
429: PetscCall(VecDestroy(&ldb->W));
430: ldb->allocated = PETSC_FALSE;
431: }
432: PetscCall(PetscFree(lmvm->ctx));
433: PetscCall(MatDestroy_LMVM(B));
434: PetscFunctionReturn(PETSC_SUCCESS);
435: }
437: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
438: {
439: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
440: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
442: PetscFunctionBegin;
443: PetscCall(MatSetUp_LMVM(B));
444: if (!ldb->allocated) {
445: PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
446: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
447: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
448: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
449: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
450: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
451: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
452: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
453: PetscCall(VecSet(ldb->invD, ldb->delta));
454: ldb->allocated = PETSC_TRUE;
455: }
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
460: {
461: Mat_LMVM *lmvm;
462: Mat_DiagBrdn *ldb;
464: PetscFunctionBegin;
465: PetscCall(MatCreate_LMVM(B));
466: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN));
467: B->ops->setup = MatSetUp_DiagBrdn;
468: B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
469: B->ops->destroy = MatDestroy_DiagBrdn;
470: B->ops->view = MatView_DiagBrdn;
472: lmvm = (Mat_LMVM *)B->data;
473: lmvm->square = PETSC_TRUE;
474: lmvm->m = 1;
475: lmvm->ops->allocate = MatAllocate_DiagBrdn;
476: lmvm->ops->reset = MatReset_DiagBrdn;
477: lmvm->ops->mult = MatMult_DiagBrdn;
478: lmvm->ops->solve = MatSolve_DiagBrdn;
479: lmvm->ops->update = MatUpdate_DiagBrdn;
480: lmvm->ops->copy = MatCopy_DiagBrdn;
482: PetscCall(PetscNew(&ldb));
483: lmvm->ctx = (void *)ldb;
484: ldb->theta = 0.0;
485: ldb->alpha = 1.0;
486: ldb->rho = 1.0;
487: ldb->forward = PETSC_TRUE;
488: ldb->beta = 0.5;
489: ldb->sigma = 1.0;
490: ldb->delta = 1.0;
491: ldb->delta_min = 1e-7;
492: ldb->delta_max = 100.0;
493: ldb->tol = 1e-8;
494: ldb->sigma_hist = 1;
495: ldb->allocated = PETSC_FALSE;
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: /*@
500: MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
501: for approximating Hessians.
503: Collective
505: Input Parameters:
506: + comm - MPI communicator
507: . n - number of local rows for storage vectors
508: - N - global size of the storage vectors
510: Output Parameter:
511: . B - the matrix
513: Options Database Keys:
514: + -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
515: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
516: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
517: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
518: . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
519: . -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
520: - -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal
522: Level: intermediate
524: Notes:
525: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
526: paradigm instead of this routine directly.
528: It consists of a convex combination of DFP and BFGS
529: diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
530: To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
531: We also ensure positive definiteness by taking the `VecAbs()` of the final vector.
533: There are two ways of approximating the diagonal: using the forward (B) update
534: schemes for BFGS and DFP and then taking the inverse, or directly working with
535: the inverse (H) update schemes for the BFGS and DFP updates, derived using the
536: Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a
537: parameter below.
539: In order to use the DiagBrdn matrix with other vector types, i.e. doing matrix-vector products
540: and matrix solves, the matrix must first be created using `MatCreate()` and `MatSetType()`,
541: followed by `MatLMVMAllocate()`. Then it will be available for updating
542: (via `MatLMVMUpdate()`) in one's favored solver implementation.
544: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDIAGBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
545: `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMSymBrdn()`
546: @*/
547: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
548: {
549: PetscFunctionBegin;
550: PetscCall(KSPInitializePackage());
551: PetscCall(MatCreate(comm, B));
552: PetscCall(MatSetSizes(*B, n, n, N, N));
553: PetscCall(MatSetType(*B, MATLMVMDIAGBROYDEN));
554: PetscCall(MatSetUp(*B));
555: PetscFunctionReturn(PETSC_SUCCESS);
556: }