Actual source code: diagbrdn.c
1: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
3: static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
4: {
5: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
6: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
8: PetscFunctionBegin;
9: PetscCall(VecPointwiseMult(dX, ldb->invD, F));
10: PetscFunctionReturn(PETSC_SUCCESS);
11: }
13: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
14: {
15: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
16: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
18: PetscFunctionBegin;
19: PetscCall(VecPointwiseDivide(Z, X, ldb->invD));
20: PetscFunctionReturn(PETSC_SUCCESS);
21: }
23: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
24: {
25: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
26: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
27: PetscInt old_k, i, start;
28: PetscScalar yty, curvature, ytDy, stDs, ytDs;
29: PetscReal curvtol, sigma, yy_sum, ss_sum, ys_sum, denom, ststmp;
30: PetscReal stDsr, ytDyr;
32: PetscFunctionBegin;
33: if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
34: if (lmvm->prev_set) {
35: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
36: PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
37: PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
39: /* Test if the updates can be accepted */
40: PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ststmp));
41: if (ststmp < lmvm->eps) curvtol = 0.0;
42: else curvtol = lmvm->eps * ststmp;
44: /* Test the curvature for the update */
45: if (PetscRealPart(curvature) > curvtol) {
46: /* Update is good so we accept it */
47: old_k = lmvm->k;
48: PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
49: /* If we hit the memory limit, shift the yty and yts arrays */
50: if (old_k == lmvm->k) {
51: for (i = 0; i <= lmvm->k - 1; ++i) {
52: ldb->yty[i] = ldb->yty[i + 1];
53: ldb->yts[i] = ldb->yts[i + 1];
54: ldb->sts[i] = ldb->sts[i + 1];
55: }
56: }
57: /* Accept dot products into the history */
58: PetscCall(VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty));
59: ldb->yty[lmvm->k] = PetscRealPart(yty);
60: ldb->yts[lmvm->k] = PetscRealPart(curvature);
61: ldb->sts[lmvm->k] = ststmp;
62: if (ldb->forward) {
63: /* We are doing diagonal scaling of the forward Hessian B */
64: /* BFGS = DFP = inv(D); */
65: PetscCall(VecCopy(ldb->invD, ldb->invDnew));
66: PetscCall(VecReciprocal(ldb->invDnew));
68: /* V = y*y */
69: PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]));
71: /* W = inv(D)*s */
72: PetscCall(VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]));
73: PetscCall(VecDot(ldb->W, lmvm->S[lmvm->k], &stDs));
75: /* Safeguard stDs */
76: stDs = PetscMax(PetscRealPart(stDs), ldb->tol);
78: if (1.0 != ldb->theta) {
79: /* BFGS portion of the update */
80: /* U = (inv(D)*s)*(inv(D)*s) */
81: PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));
83: /* Assemble */
84: PetscCall(VecAXPBY(ldb->BFGS, -1.0 / stDs, 0.0, ldb->U));
85: }
86: if (0.0 != ldb->theta) {
87: /* DFP portion of the update */
88: /* U = inv(D)*s*y */
89: PetscCall(VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]));
91: /* Assemble */
92: PetscCall(VecAXPBY(ldb->DFP, stDs / ldb->yts[lmvm->k], 0.0, ldb->V));
93: PetscCall(VecAXPY(ldb->DFP, -2.0, ldb->U));
94: }
96: if (0.0 == ldb->theta) {
97: PetscCall(VecAXPY(ldb->invDnew, 1.0, ldb->BFGS));
98: } else if (1.0 == ldb->theta) {
99: PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->DFP));
100: } else {
101: /* Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
102: PetscCall(VecAXPBYPCZ(ldb->invDnew, 1.0 - ldb->theta, (ldb->theta) / ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP));
103: }
105: PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
106: /* Obtain inverse and ensure positive definite */
107: PetscCall(VecReciprocal(ldb->invDnew));
108: PetscCall(VecAbs(ldb->invDnew));
110: } else {
111: /* Inverse Hessian update instead. */
112: PetscCall(VecCopy(ldb->invD, ldb->invDnew));
114: /* V = s*s */
115: PetscCall(VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]));
117: /* W = D*y */
118: PetscCall(VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]));
119: PetscCall(VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy));
121: /* Safeguard ytDy */
122: ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);
124: if (1.0 != ldb->theta) {
125: /* BFGS portion of the update */
126: /* U = s*Dy */
127: PetscCall(VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]));
129: /* Assemble */
130: PetscCall(VecAXPBY(ldb->BFGS, ytDy / ldb->yts[lmvm->k], 0.0, ldb->V));
131: PetscCall(VecAXPY(ldb->BFGS, -2.0, ldb->U));
132: }
133: if (0.0 != ldb->theta) {
134: /* DFP portion of the update */
136: /* U = (inv(D)*y)*(inv(D)*y) */
137: PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));
139: /* Assemble */
140: PetscCall(VecAXPBY(ldb->DFP, -1.0 / ytDy, 0.0, ldb->U));
141: }
143: if (0.0 == ldb->theta) {
144: PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->BFGS));
145: } else if (1.0 == ldb->theta) {
146: PetscCall(VecAXPY(ldb->invDnew, 1.0, ldb->DFP));
147: } else {
148: /* Broyden update U=(1-theta)*P + theta*Q */
149: PetscCall(VecAXPBYPCZ(ldb->invDnew, (1.0 - ldb->theta) / ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP));
150: }
151: PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
152: /* Ensure positive definite */
153: PetscCall(VecAbs(ldb->invDnew));
154: }
155: if (ldb->sigma_hist > 0) {
156: /* Start with re-scaling on the newly computed diagonal */
157: if (0.5 == ldb->beta) {
158: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
159: PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew));
160: PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));
162: PetscCall(VecDot(ldb->V, lmvm->Y[0], &ytDy));
163: PetscCall(VecDot(ldb->W, lmvm->S[0], &stDs));
165: ss_sum = PetscRealPart(stDs);
166: yy_sum = PetscRealPart(ytDy);
167: ys_sum = ldb->yts[0];
168: } else {
169: PetscCall(VecCopy(ldb->invDnew, ldb->U));
170: PetscCall(VecReciprocal(ldb->U));
172: /* Compute summations for scalar scaling */
173: yy_sum = 0; /* No safeguard required */
174: ys_sum = 0; /* No safeguard required */
175: ss_sum = 0; /* No safeguard required */
176: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
177: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
178: PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U));
179: PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));
181: PetscCall(VecDot(ldb->W, lmvm->S[i], &stDs));
182: PetscCall(VecDot(ldb->V, lmvm->Y[i], &ytDy));
184: ss_sum += PetscRealPart(stDs);
185: ys_sum += ldb->yts[i];
186: yy_sum += PetscRealPart(ytDy);
187: }
188: }
189: } else if (0.0 == ldb->beta) {
190: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
191: /* Compute summations for scalar scaling */
192: PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));
194: PetscCall(VecDotNorm2(lmvm->Y[0], ldb->W, &ytDs, &stDsr));
196: ys_sum = PetscRealPart(ytDs);
197: ss_sum = stDsr;
198: yy_sum = ldb->yty[0];
199: } else {
200: PetscCall(VecCopy(ldb->invDnew, ldb->U));
201: PetscCall(VecReciprocal(ldb->U));
203: /* Compute summations for scalar scaling */
204: yy_sum = 0; /* No safeguard required */
205: ys_sum = 0; /* No safeguard required */
206: ss_sum = 0; /* No safeguard required */
207: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
208: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
209: PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));
211: PetscCall(VecDotNorm2(lmvm->Y[i], ldb->W, &ytDs, &stDsr));
213: ss_sum += stDsr;
214: ys_sum += PetscRealPart(ytDs);
215: yy_sum += ldb->yty[i];
216: }
217: }
218: } else if (1.0 == ldb->beta) {
219: /* Compute summations for scalar scaling */
220: yy_sum = 0; /* No safeguard required */
221: ys_sum = 0; /* No safeguard required */
222: ss_sum = 0; /* No safeguard required */
223: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
224: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
225: PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew));
227: PetscCall(VecDotNorm2(lmvm->S[i], ldb->V, &ytDs, &ytDyr));
229: yy_sum += ytDyr;
230: ys_sum += PetscRealPart(ytDs);
231: ss_sum += ldb->sts[i];
232: }
233: } else {
234: PetscCall(VecCopy(ldb->invDnew, ldb->U));
235: PetscCall(VecPow(ldb->U, ldb->beta - 1));
237: /* Compute summations for scalar scaling */
238: yy_sum = 0; /* No safeguard required */
239: ys_sum = 0; /* No safeguard required */
240: ss_sum = 0; /* No safeguard required */
241: start = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
242: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
243: PetscCall(VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]));
244: PetscCall(VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]));
246: PetscCall(VecDotNorm2(ldb->W, ldb->V, &ytDs, &ytDyr));
247: PetscCall(VecDot(ldb->W, ldb->W, &stDs));
249: yy_sum += ytDyr;
250: ys_sum += PetscRealPart(ytDs);
251: ss_sum += PetscRealPart(stDs);
252: }
253: }
255: if (0.0 == ldb->alpha) {
256: /* Safeguard ys_sum */
257: ys_sum = PetscMax(ldb->tol, ys_sum);
259: sigma = ss_sum / ys_sum;
260: } else if (1.0 == ldb->alpha) {
261: /* yy_sum is never 0; if it were, we'd be at the minimum */
262: sigma = ys_sum / yy_sum;
263: } else {
264: denom = 2.0 * ldb->alpha * yy_sum;
266: /* Safeguard denom */
267: denom = PetscMax(ldb->tol, denom);
269: 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;
270: }
271: } else {
272: sigma = 1.0;
273: }
274: /* If Q has small values, then Q^(r_beta - 1)
275: can have very large values. Hence, ys_sum
276: and ss_sum can be infinity. In this case,
277: sigma can either be not-a-number or infinity. */
279: if (PetscIsInfOrNanScalar(sigma)) {
280: /* sigma is not-a-number; skip rescaling */
281: } else if (0.0 == sigma) {
282: /* sigma is zero; this is a bad case; skip rescaling */
283: } else {
284: /* sigma is positive */
285: PetscCall(VecScale(ldb->invDnew, sigma));
286: }
288: /* Combine the old diagonal and the new diagonal using a convex limiter */
289: if (1.0 == ldb->rho) {
290: PetscCall(VecCopy(ldb->invDnew, ldb->invD));
291: } else if (ldb->rho) PetscCall(VecAXPBY(ldb->invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew));
292: } else {
293: PetscCall(MatLMVMReset(B, PETSC_FALSE));
294: }
295: /* End DiagBrdn update */
296: }
297: /* Save the solution and function to be used in the next update */
298: PetscCall(VecCopy(X, lmvm->Xprev));
299: PetscCall(VecCopy(F, lmvm->Fprev));
300: lmvm->prev_set = PETSC_TRUE;
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
304: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
305: {
306: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
307: Mat_DiagBrdn *bctx = (Mat_DiagBrdn *)bdata->ctx;
308: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
309: Mat_DiagBrdn *mctx = (Mat_DiagBrdn *)mdata->ctx;
310: PetscInt i;
312: PetscFunctionBegin;
313: mctx->theta = bctx->theta;
314: mctx->alpha = bctx->alpha;
315: mctx->beta = bctx->beta;
316: mctx->rho = bctx->rho;
317: mctx->delta = bctx->delta;
318: mctx->delta_min = bctx->delta_min;
319: mctx->delta_max = bctx->delta_max;
320: mctx->tol = bctx->tol;
321: mctx->sigma = bctx->sigma;
322: mctx->sigma_hist = bctx->sigma_hist;
323: mctx->forward = bctx->forward;
324: PetscCall(VecCopy(bctx->invD, mctx->invD));
325: for (i = 0; i <= bdata->k; ++i) {
326: mctx->yty[i] = bctx->yty[i];
327: mctx->yts[i] = bctx->yts[i];
328: mctx->sts[i] = bctx->sts[i];
329: }
330: PetscFunctionReturn(PETSC_SUCCESS);
331: }
333: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
334: {
335: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
336: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
337: PetscBool isascii;
339: PetscFunctionBegin;
340: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
341: if (isascii) {
342: PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", ldb->sigma_hist));
343: PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho));
344: PetscCall(PetscViewerASCIIPrintf(pv, "Convex factor: theta=%g\n", (double)ldb->theta));
345: }
346: PetscCall(MatView_LMVM(B, pv));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: static PetscErrorCode MatSetFromOptions_DiagBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
351: {
352: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
353: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
355: PetscFunctionBegin;
356: PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
357: PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
358: PetscCall(PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL));
359: PetscCall(PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL));
360: PetscCall(PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL));
361: PetscCall(PetscOptionsReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL));
362: PetscCall(PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL));
363: PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL));
364: PetscCall(PetscOptionsInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", ldb->sigma_hist, &ldb->sigma_hist, NULL));
365: PetscOptionsHeadEnd();
366: PetscCheck(!(ldb->theta < 0.0) && !(ldb->theta > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
367: PetscCheck(!(ldb->alpha < 0.0) && !(ldb->alpha > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
368: PetscCheck(!(ldb->rho < 0.0) && !(ldb->rho > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex update limiter in the J0 scaling cannot be outside the range of [0, 1]");
369: PetscCheck(ldb->sigma_hist >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
374: {
375: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
376: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
378: PetscFunctionBegin;
379: PetscCall(VecSet(ldb->invD, ldb->delta));
380: if (destructive && ldb->allocated) {
381: PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
382: PetscCall(VecDestroy(&ldb->invDnew));
383: PetscCall(VecDestroy(&ldb->invD));
384: PetscCall(VecDestroy(&ldb->BFGS));
385: PetscCall(VecDestroy(&ldb->DFP));
386: PetscCall(VecDestroy(&ldb->U));
387: PetscCall(VecDestroy(&ldb->V));
388: PetscCall(VecDestroy(&ldb->W));
389: ldb->allocated = PETSC_FALSE;
390: }
391: PetscCall(MatReset_LMVM(B, destructive));
392: PetscFunctionReturn(PETSC_SUCCESS);
393: }
395: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
396: {
397: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
398: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
400: PetscFunctionBegin;
401: PetscCall(MatAllocate_LMVM(B, X, F));
402: if (!ldb->allocated) {
403: PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
404: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
405: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
406: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
407: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
408: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
409: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
410: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
411: ldb->allocated = PETSC_TRUE;
412: }
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
417: {
418: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
419: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
421: PetscFunctionBegin;
422: if (ldb->allocated) {
423: PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
424: PetscCall(VecDestroy(&ldb->invDnew));
425: PetscCall(VecDestroy(&ldb->invD));
426: PetscCall(VecDestroy(&ldb->BFGS));
427: PetscCall(VecDestroy(&ldb->DFP));
428: PetscCall(VecDestroy(&ldb->U));
429: PetscCall(VecDestroy(&ldb->V));
430: PetscCall(VecDestroy(&ldb->W));
431: ldb->allocated = PETSC_FALSE;
432: }
433: PetscCall(PetscFree(lmvm->ctx));
434: PetscCall(MatDestroy_LMVM(B));
435: PetscFunctionReturn(PETSC_SUCCESS);
436: }
438: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
439: {
440: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
441: Mat_DiagBrdn *ldb = (Mat_DiagBrdn *)lmvm->ctx;
443: PetscFunctionBegin;
444: PetscCall(MatSetUp_LMVM(B));
445: if (!ldb->allocated) {
446: PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
447: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
448: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
449: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
450: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
451: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
452: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
453: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
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->solve = MatSolve_DiagBrdn;
471: B->ops->view = MatView_DiagBrdn;
473: lmvm = (Mat_LMVM *)B->data;
474: lmvm->square = PETSC_TRUE;
475: lmvm->m = 1;
476: lmvm->ops->allocate = MatAllocate_DiagBrdn;
477: lmvm->ops->reset = MatReset_DiagBrdn;
478: lmvm->ops->mult = MatMult_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(MatCreate(comm, B));
551: PetscCall(MatSetSizes(*B, n, n, N, N));
552: PetscCall(MatSetType(*B, MATLMVMDIAGBROYDEN));
553: PetscCall(MatSetUp(*B));
554: PetscFunctionReturn(PETSC_SUCCESS);
555: }