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