Actual source code: symbrdn.c
1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
2: #include <../src/ksp/ksp/utils/lmvm/dense/denseqn.h>
3: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
4: #include <petsc/private/kspimpl.h>
5: #include <petscdevice.h>
7: const char *const MatLMVMSymBroydenScaleTypes[] = {"NONE", "SCALAR", "DIAGONAL", "USER", "MatLMVMSymBrdnScaleType", "MAT_LMVM_SYMBROYDEN_SCALING_", NULL};
9: /*
10: The solution method below is the matrix-free implementation of
11: Equation 8.6a in Dennis and More "Quasi-Newton Methods, Motivation
12: and Theory" (https://epubs.siam.org/doi/abs/10.1137/1019005).
14: Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
15: the matrix is updated with a new (S[i], Y[i]) pair. This allows
16: repeated calls of MatSolve without incurring redundant computation.
18: dX <- J0^{-1} * F
20: for i=0,1,2,...,k
21: # Q[i] = (B_i)^T{-1} Y[i]
23: rho = 1.0 / (Y[i]^T S[i])
24: alpha = rho * (S[i]^T F)
25: zeta = 1.0 / (Y[i]^T Q[i])
26: gamma = zeta * (Y[i]^T dX)
28: dX <- dX - (gamma * Q[i]) + (alpha * Y[i])
29: W <- (rho * S[i]) - (zeta * Q[i])
30: dX <- dX + (psi[i] * (Y[i]^T Q[i]) * (W^T F) * W)
31: end
32: */
33: static PetscErrorCode MatSolve_LMVMSymBrdn(Mat B, Vec F, Vec dX)
34: {
35: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
36: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
37: PetscInt i, j;
38: PetscReal numer;
39: PetscScalar sjtpi, yjtsi, wtsi, yjtqi, sjtyi, wtyi, ytx, stf, wtf, stp, ytq;
41: PetscFunctionBegin;
42: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
43: if (lsb->phi == 0.0) {
44: PetscCall(MatSolve_LMVMBFGS(B, F, dX));
45: PetscFunctionReturn(PETSC_SUCCESS);
46: }
47: if (lsb->phi == 1.0) {
48: PetscCall(MatSolve_LMVMDFP(B, F, dX));
49: PetscFunctionReturn(PETSC_SUCCESS);
50: }
52: VecCheckSameSize(F, 2, dX, 3);
53: VecCheckMatCompatible(B, dX, 3, F, 2);
55: if (lsb->needP) {
56: /* Start the loop for (P[k] = (B_k) * S[k]) */
57: for (i = 0; i <= lmvm->k; ++i) {
58: PetscCall(MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]));
59: /* Compute the necessary dot products */
60: PetscCall(VecMDot(lmvm->S[i], i, lmvm->Y, lsb->workscalar));
61: for (j = 0; j < i; ++j) {
62: yjtsi = lsb->workscalar[j];
63: PetscCall(VecDot(lmvm->S[j], lsb->P[i], &sjtpi));
64: /* Compute the pure BFGS component of the forward product */
65: PetscCall(VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi) / lsb->stp[j], PetscRealPart(yjtsi) / lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]));
66: /* Tack on the convexly scaled extras to the forward product */
67: if (lsb->phi > 0.0) {
68: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]));
69: PetscCall(VecDot(lsb->work, lmvm->S[i], &wtsi));
70: PetscCall(VecAXPY(lsb->P[i], lsb->phi * lsb->stp[j] * PetscRealPart(wtsi), lsb->work));
71: }
72: }
73: PetscCall(VecDot(lmvm->S[i], lsb->P[i], &stp));
74: lsb->stp[i] = PetscRealPart(stp);
75: }
76: lsb->needP = PETSC_FALSE;
77: }
78: if (lsb->needQ) {
79: /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
80: for (i = 0; i <= lmvm->k; ++i) {
81: PetscCall(MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]));
82: /* Compute the necessary dot products */
83: PetscCall(VecMDot(lmvm->Y[i], i, lmvm->S, lsb->workscalar));
84: for (j = 0; j < i; ++j) {
85: sjtyi = lsb->workscalar[j];
86: PetscCall(VecDot(lmvm->Y[j], lsb->Q[i], &yjtqi));
87: /* Compute the pure DFP component of the inverse application*/
88: PetscCall(VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi) / lsb->ytq[j], PetscRealPart(sjtyi) / lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]));
89: /* Tack on the convexly scaled extras to the inverse application*/
90: if (lsb->psi[j] > 0.0) {
91: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]));
92: PetscCall(VecDot(lsb->work, lmvm->Y[i], &wtyi));
93: PetscCall(VecAXPY(lsb->Q[i], lsb->psi[j] * lsb->ytq[j] * PetscRealPart(wtyi), lsb->work));
94: }
95: }
96: PetscCall(VecDot(lmvm->Y[i], lsb->Q[i], &ytq));
97: lsb->ytq[i] = PetscRealPart(ytq);
98: if (lsb->phi == 1.0) {
99: lsb->psi[i] = 0.0;
100: } else if (lsb->phi == 0.0) {
101: lsb->psi[i] = 1.0;
102: } else {
103: numer = (1.0 - lsb->phi) * lsb->yts[i] * lsb->yts[i];
104: lsb->psi[i] = numer / (numer + (lsb->phi * lsb->ytq[i] * lsb->stp[i]));
105: }
106: }
107: lsb->needQ = PETSC_FALSE;
108: }
110: /* Start the outer iterations for ((B^{-1}) * dX) */
111: PetscCall(MatSymBrdnApplyJ0Inv(B, F, dX));
112: /* Get all the dot products we need */
113: PetscCall(VecMDot(F, lmvm->k + 1, lmvm->S, lsb->workscalar));
114: for (i = 0; i <= lmvm->k; ++i) {
115: stf = lsb->workscalar[i];
116: PetscCall(VecDot(lmvm->Y[i], dX, &ytx));
117: /* Compute the pure DFP component */
118: PetscCall(VecAXPBYPCZ(dX, -PetscRealPart(ytx) / lsb->ytq[i], PetscRealPart(stf) / lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]));
119: /* Tack on the convexly scaled extras */
120: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]));
121: PetscCall(VecDot(lsb->work, F, &wtf));
122: PetscCall(VecAXPY(dX, lsb->psi[i] * lsb->ytq[i] * PetscRealPart(wtf), lsb->work));
123: }
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: /*
128: The forward-product below is the matrix-free implementation of
129: Equation 16 in Dennis and Wolkowicz "Sizing and Least Change Secant
130: Methods" (http://www.caam.rice.edu/caam/trs/90/TR90-05.pdf).
132: P[i] = (B_i)*S[i] terms are computed ahead of time whenever
133: the matrix is updated with a new (S[i], Y[i]) pair. This allows
134: repeated calls of MatMult inside KSP solvers without unnecessarily
135: recomputing P[i] terms in expensive nested-loops.
137: Z <- J0 * X
139: for i=0,1,2,...,k
140: # P[i] = (B_k) * S[i]
142: rho = 1.0 / (Y[i]^T S[i])
143: alpha = rho * (Y[i]^T F)
144: zeta = 1.0 / (S[i]^T P[i])
145: gamma = zeta * (S[i]^T dX)
147: dX <- dX - (gamma * P[i]) + (alpha * S[i])
148: W <- (rho * Y[i]) - (zeta * P[i])
149: dX <- dX + (phi * (S[i]^T P[i]) * (W^T F) * W)
150: end
151: */
152: static PetscErrorCode MatMult_LMVMSymBrdn(Mat B, Vec X, Vec Z)
153: {
154: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
155: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
156: PetscInt i, j;
157: PetscScalar sjtpi, yjtsi, wtsi, stz, ytx, wtx, stp;
159: PetscFunctionBegin;
160: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
161: if (lsb->phi == 0.0) {
162: PetscCall(MatMult_LMVMBFGS(B, X, Z));
163: PetscFunctionReturn(PETSC_SUCCESS);
164: }
165: if (lsb->phi == 1.0) {
166: PetscCall(MatMult_LMVMDFP(B, X, Z));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: VecCheckSameSize(X, 2, Z, 3);
171: VecCheckMatCompatible(B, X, 2, Z, 3);
173: if (lsb->needP) {
174: /* Start the loop for (P[k] = (B_k) * S[k]) */
175: for (i = 0; i <= lmvm->k; ++i) {
176: PetscCall(MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]));
177: /* Compute the necessary dot products */
178: PetscCall(VecMDot(lmvm->S[i], i, lmvm->Y, lsb->workscalar));
179: for (j = 0; j < i; ++j) {
180: yjtsi = lsb->workscalar[j];
181: PetscCall(VecDot(lmvm->S[j], lsb->P[i], &sjtpi));
182: /* Compute the pure BFGS component of the forward product */
183: PetscCall(VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi) / lsb->stp[j], PetscRealPart(yjtsi) / lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]));
184: /* Tack on the convexly scaled extras to the forward product */
185: if (lsb->phi > 0.0) {
186: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]));
187: PetscCall(VecDot(lsb->work, lmvm->S[i], &wtsi));
188: PetscCall(VecAXPY(lsb->P[i], lsb->phi * lsb->stp[j] * PetscRealPart(wtsi), lsb->work));
189: }
190: }
191: PetscCall(VecDot(lmvm->S[i], lsb->P[i], &stp));
192: lsb->stp[i] = PetscRealPart(stp);
193: }
194: lsb->needP = PETSC_FALSE;
195: }
197: /* Start the outer iterations for (B * X) */
198: PetscCall(MatSymBrdnApplyJ0Fwd(B, X, Z));
199: /* Get all the dot products we need */
200: PetscCall(VecMDot(X, lmvm->k + 1, lmvm->Y, lsb->workscalar));
201: for (i = 0; i <= lmvm->k; ++i) {
202: ytx = lsb->workscalar[i];
203: PetscCall(VecDot(lmvm->S[i], Z, &stz));
204: /* Compute the pure BFGS component */
205: PetscCall(VecAXPBYPCZ(Z, -PetscRealPart(stz) / lsb->stp[i], PetscRealPart(ytx) / lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]));
206: /* Tack on the convexly scaled extras */
207: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]));
208: PetscCall(VecDot(lsb->work, X, &wtx));
209: PetscCall(VecAXPY(Z, lsb->phi * lsb->stp[i] * PetscRealPart(wtx), lsb->work));
210: }
211: PetscFunctionReturn(PETSC_SUCCESS);
212: }
214: static PetscErrorCode MatUpdate_LMVMSymBrdn(Mat B, Vec X, Vec F)
215: {
216: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
217: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
218: Mat_LMVM *dbase;
219: Mat_DiagBrdn *dctx;
220: PetscInt old_k, i;
221: PetscReal curvtol, ytytmp;
222: PetscScalar curvature, ststmp;
224: PetscFunctionBegin;
225: if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
226: if (lmvm->prev_set) {
227: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
228: PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
229: PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
231: /* Test if the updates can be accepted */
232: PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ytytmp));
233: if (ytytmp < lmvm->eps) curvtol = 0.0;
234: else curvtol = lmvm->eps * ytytmp;
236: if (PetscRealPart(curvature) > curvtol) {
237: /* Update is good, accept it */
238: lsb->watchdog = 0;
239: lsb->needP = lsb->needQ = PETSC_TRUE;
240: old_k = lmvm->k;
241: PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
242: /* If we hit the memory limit, shift the yts, yty and sts arrays */
243: if (old_k == lmvm->k) {
244: for (i = 0; i <= lmvm->k - 1; ++i) {
245: lsb->yts[i] = lsb->yts[i + 1];
246: lsb->yty[i] = lsb->yty[i + 1];
247: lsb->sts[i] = lsb->sts[i + 1];
248: }
249: }
250: /* Update history of useful scalars */
251: PetscCall(VecDot(lmvm->S[lmvm->k], lmvm->S[lmvm->k], &ststmp));
252: lsb->yts[lmvm->k] = PetscRealPart(curvature);
253: lsb->yty[lmvm->k] = ytytmp;
254: lsb->sts[lmvm->k] = PetscRealPart(ststmp);
255: /* Compute the scalar scale if necessary */
256: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(MatSymBrdnComputeJ0Scalar(B));
257: } else {
258: /* Update is bad, skip it */
259: ++lmvm->nrejects;
260: ++lsb->watchdog;
261: }
262: } else {
263: switch (lsb->scale_type) {
264: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
265: dbase = (Mat_LMVM *)lsb->D->data;
266: dctx = (Mat_DiagBrdn *)dbase->ctx;
267: PetscCall(VecSet(dctx->invD, lsb->delta));
268: break;
269: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
270: lsb->sigma = lsb->delta;
271: break;
272: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
273: lsb->sigma = 1.0;
274: break;
275: default:
276: break;
277: }
278: }
280: /* Update the scaling */
281: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMUpdate(lsb->D, X, F));
283: if (lsb->watchdog > lsb->max_seq_rejects) {
284: PetscCall(MatLMVMReset(B, PETSC_FALSE));
285: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMReset(lsb->D, PETSC_FALSE));
286: }
288: /* Save the solution and function to be used in the next update */
289: PetscCall(VecCopy(X, lmvm->Xprev));
290: PetscCall(VecCopy(F, lmvm->Fprev));
291: lmvm->prev_set = PETSC_TRUE;
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
296: {
297: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
298: Mat_SymBrdn *blsb = (Mat_SymBrdn *)bdata->ctx;
299: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
300: Mat_SymBrdn *mlsb = (Mat_SymBrdn *)mdata->ctx;
301: PetscInt i;
303: PetscFunctionBegin;
304: mlsb->phi = blsb->phi;
305: mlsb->needP = blsb->needP;
306: mlsb->needQ = blsb->needQ;
307: for (i = 0; i <= bdata->k; ++i) {
308: mlsb->stp[i] = blsb->stp[i];
309: mlsb->ytq[i] = blsb->ytq[i];
310: mlsb->yts[i] = blsb->yts[i];
311: mlsb->psi[i] = blsb->psi[i];
312: PetscCall(VecCopy(blsb->P[i], mlsb->P[i]));
313: PetscCall(VecCopy(blsb->Q[i], mlsb->Q[i]));
314: }
315: mlsb->scale_type = blsb->scale_type;
316: mlsb->alpha = blsb->alpha;
317: mlsb->beta = blsb->beta;
318: mlsb->rho = blsb->rho;
319: mlsb->delta = blsb->delta;
320: mlsb->sigma_hist = blsb->sigma_hist;
321: mlsb->watchdog = blsb->watchdog;
322: mlsb->max_seq_rejects = blsb->max_seq_rejects;
323: switch (blsb->scale_type) {
324: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
325: mlsb->sigma = blsb->sigma;
326: break;
327: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
328: PetscCall(MatCopy(blsb->D, mlsb->D, SAME_NONZERO_PATTERN));
329: break;
330: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
331: mlsb->sigma = 1.0;
332: break;
333: default:
334: break;
335: }
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
340: {
341: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
342: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
343: Mat_LMVM *dbase;
344: Mat_DiagBrdn *dctx;
346: PetscFunctionBegin;
347: lsb->watchdog = 0;
348: lsb->needP = lsb->needQ = PETSC_TRUE;
349: if (lsb->allocated) {
350: if (destructive) {
351: PetscCall(VecDestroy(&lsb->work));
352: PetscCall(PetscFree6(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts, lsb->workscalar));
353: PetscCall(PetscFree(lsb->psi));
354: PetscCall(VecDestroyVecs(lmvm->m, &lsb->P));
355: PetscCall(VecDestroyVecs(lmvm->m, &lsb->Q));
356: switch (lsb->scale_type) {
357: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
358: PetscCall(MatLMVMReset(lsb->D, PETSC_TRUE));
359: break;
360: default:
361: break;
362: }
363: lsb->allocated = PETSC_FALSE;
364: } else {
365: PetscCall(PetscMemzero(lsb->psi, lmvm->m));
366: switch (lsb->scale_type) {
367: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
368: lsb->sigma = lsb->delta;
369: break;
370: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
371: PetscCall(MatLMVMReset(lsb->D, PETSC_FALSE));
372: dbase = (Mat_LMVM *)lsb->D->data;
373: dctx = (Mat_DiagBrdn *)dbase->ctx;
374: PetscCall(VecSet(dctx->invD, lsb->delta));
375: break;
376: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
377: lsb->sigma = 1.0;
378: break;
379: default:
380: break;
381: }
382: }
383: }
384: PetscCall(MatReset_LMVM(B, destructive));
385: PetscFunctionReturn(PETSC_SUCCESS);
386: }
388: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
389: {
390: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
391: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
393: PetscFunctionBegin;
394: PetscCall(MatAllocate_LMVM(B, X, F));
395: if (!lsb->allocated) {
396: PetscCall(VecDuplicate(X, &lsb->work));
397: if (lmvm->m > 0) {
398: PetscCall(PetscMalloc6(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts, lmvm->m, &lsb->workscalar));
399: PetscCall(PetscCalloc1(lmvm->m, &lsb->psi));
400: PetscCall(VecDuplicateVecs(X, lmvm->m, &lsb->P));
401: PetscCall(VecDuplicateVecs(X, lmvm->m, &lsb->Q));
402: }
403: switch (lsb->scale_type) {
404: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
405: PetscCall(MatLMVMAllocate(lsb->D, X, F));
406: break;
407: default:
408: break;
409: }
410: lsb->allocated = PETSC_TRUE;
411: }
412: PetscFunctionReturn(PETSC_SUCCESS);
413: }
415: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
416: {
417: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
418: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
420: PetscFunctionBegin;
421: if (lsb->allocated) {
422: PetscCall(VecDestroy(&lsb->work));
423: PetscCall(PetscFree6(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts, lsb->workscalar));
424: PetscCall(PetscFree(lsb->psi));
425: PetscCall(VecDestroyVecs(lmvm->m, &lsb->P));
426: PetscCall(VecDestroyVecs(lmvm->m, &lsb->Q));
427: lsb->allocated = PETSC_FALSE;
428: }
429: PetscCall(MatDestroy(&lsb->D));
430: PetscCall(PetscFree(lmvm->ctx));
431: PetscCall(MatDestroy_LMVM(B));
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
436: {
437: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
438: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
439: PetscInt n, N;
441: PetscFunctionBegin;
442: PetscCall(MatSetUp_LMVM(B));
443: if (!lsb->allocated) {
444: PetscCall(VecDuplicate(lmvm->Xprev, &lsb->work));
445: if (lmvm->m > 0) {
446: PetscCall(PetscMalloc6(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts, lmvm->m, &lsb->workscalar));
447: PetscCall(PetscCalloc1(lmvm->m, &lsb->psi));
448: PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->P));
449: PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->Q));
450: }
451: switch (lsb->scale_type) {
452: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
453: PetscCall(MatGetLocalSize(B, &n, &n));
454: PetscCall(MatGetSize(B, &N, &N));
455: PetscCall(MatSetSizes(lsb->D, n, n, N, N));
456: PetscCall(MatSetUp(lsb->D));
457: break;
458: default:
459: break;
460: }
461: lsb->allocated = PETSC_TRUE;
462: }
463: PetscFunctionReturn(PETSC_SUCCESS);
464: }
466: PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
467: {
468: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
469: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
470: PetscBool isascii;
472: PetscFunctionBegin;
473: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
474: if (isascii) {
475: PetscCall(PetscViewerASCIIPrintf(pv, "Scale type: %s\n", MatLMVMSymBroydenScaleTypes[lsb->scale_type]));
476: PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", lsb->sigma_hist));
477: PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)lsb->alpha, (double)lsb->beta, (double)lsb->rho));
478: PetscCall(PetscViewerASCIIPrintf(pv, "Convex factors: phi=%g, theta=%g\n", (double)lsb->phi, (double)lsb->theta));
479: }
480: PetscCall(MatView_LMVM(B, pv));
481: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatView(lsb->D, pv));
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: PetscErrorCode MatSetFromOptions_LMVMSymBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
486: {
487: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
488: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
490: PetscFunctionBegin;
491: PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
492: PetscOptionsHeadBegin(PetscOptionsObject, "Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
493: PetscCall(PetscOptionsRangeReal("-mat_lmvm_phi", "(developer) convex ratio between BFGS and DFP components of the update", "", lsb->phi, &lsb->phi, NULL, 0.0, 1.0));
494: PetscCall(MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject));
495: PetscOptionsHeadEnd();
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: PetscErrorCode MatSetFromOptions_LMVMSymBrdn_Private(Mat B, PetscOptionItems *PetscOptionsObject)
500: {
501: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
502: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
503: Mat_LMVM *dbase;
504: Mat_DiagBrdn *dctx;
505: MatLMVMSymBroydenScaleType stype = lsb->scale_type;
506: PetscBool flg;
508: PetscFunctionBegin;
509: PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", lsb->beta, &lsb->beta, NULL));
510: PetscCall(PetscOptionsRangeReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", lsb->theta, &lsb->theta, NULL, 0.0, 1.0));
511: PetscCall(PetscOptionsRangeReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", lsb->rho, &lsb->rho, NULL, 0.0, 1.0));
512: PetscCall(PetscOptionsRangeReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", lsb->alpha, &lsb->alpha, NULL, 0.0, 1.0));
513: PetscCall(PetscOptionsBoundedInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", lsb->sigma_hist, &lsb->sigma_hist, NULL, 1));
514: PetscCall(PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBrdnScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg));
515: if (flg) PetscCall(MatLMVMSymBroydenSetScaleType(B, stype));
516: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
517: const char *prefix;
519: PetscCall(MatGetOptionsPrefix(B, &prefix));
520: PetscCall(MatSetOptionsPrefix(lsb->D, prefix));
521: PetscCall(MatAppendOptionsPrefix(lsb->D, "J0_"));
522: PetscCall(MatSetFromOptions(lsb->D));
523: dbase = (Mat_LMVM *)lsb->D->data;
524: dctx = (Mat_DiagBrdn *)dbase->ctx;
525: dctx->delta_min = lsb->delta_min;
526: dctx->delta_max = lsb->delta_max;
527: dctx->theta = lsb->theta;
528: dctx->rho = lsb->rho;
529: dctx->alpha = lsb->alpha;
530: dctx->beta = lsb->beta;
531: dctx->sigma_hist = lsb->sigma_hist;
532: }
533: PetscFunctionReturn(PETSC_SUCCESS);
534: }
536: PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
537: {
538: Mat_LMVM *lmvm;
539: Mat_SymBrdn *lsb;
541: PetscFunctionBegin;
542: PetscCall(MatCreate_LMVM(B));
543: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN));
544: PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
545: PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
546: B->ops->view = MatView_LMVMSymBrdn;
547: B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
548: B->ops->setup = MatSetUp_LMVMSymBrdn;
549: B->ops->destroy = MatDestroy_LMVMSymBrdn;
551: lmvm = (Mat_LMVM *)B->data;
552: lmvm->square = PETSC_TRUE;
553: lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
554: lmvm->ops->reset = MatReset_LMVMSymBrdn;
555: lmvm->ops->update = MatUpdate_LMVMSymBrdn;
556: lmvm->ops->mult = MatMult_LMVMSymBrdn;
557: lmvm->ops->solve = MatSolve_LMVMSymBrdn;
558: lmvm->ops->copy = MatCopy_LMVMSymBrdn;
560: PetscCall(PetscNew(&lsb));
561: lmvm->ctx = (void *)lsb;
562: lsb->allocated = PETSC_FALSE;
563: lsb->needP = lsb->needQ = PETSC_TRUE;
564: lsb->phi = 0.125;
565: lsb->theta = 0.125;
566: lsb->alpha = 1.0;
567: lsb->rho = 1.0;
568: lsb->beta = 0.5;
569: lsb->sigma = 1.0;
570: lsb->delta = 1.0;
571: lsb->delta_min = 1e-7;
572: lsb->delta_max = 100.0;
573: lsb->sigma_hist = 1;
574: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
575: lsb->watchdog = 0;
576: lsb->max_seq_rejects = lmvm->m / 2;
578: PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &lsb->D));
579: PetscCall(MatSetType(lsb->D, MATLMVMDIAGBROYDEN));
580: PetscFunctionReturn(PETSC_SUCCESS);
581: }
583: /*@
584: MatLMVMSymBroydenSetDelta - Sets the starting value for the diagonal scaling vector computed
585: in the SymBrdn approximations (also works for BFGS and DFP).
587: Input Parameters:
588: + B - `MATLMVM` matrix
589: - delta - initial value for diagonal scaling
591: Level: intermediate
593: .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`
594: @*/
595: PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
596: {
597: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
598: PetscBool is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn, is_dbfgs, is_ddfp, is_dqn;
599: PetscReal del_min, del_max, del_buf;
601: PetscFunctionBegin;
602: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs));
603: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDBFGS, &is_dbfgs));
604: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDDFP, &is_ddfp));
605: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDQN, &is_dqn));
606: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp));
607: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn));
608: PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn));
610: if (is_bfgs || is_dfp || is_symbrdn || is_symbadbrdn) {
611: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
613: lsb = (Mat_SymBrdn *)lmvm->ctx;
614: del_min = lsb->delta_min;
615: del_max = lsb->delta_max;
616: } else if (is_dbfgs || is_ddfp || is_dqn) {
617: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
618: Mat_LMVM *dbase = (Mat_LMVM *)lqn->diag_qn->data;
619: Mat_DiagBrdn *diagctx = (Mat_DiagBrdn *)dbase->ctx;
621: del_min = diagctx->delta_min;
622: del_max = diagctx->delta_max;
623: } else {
624: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "diagonal scaling only available for SymBrdn-derived types (DBFGS, BFGS, DDFP, DFP, SymBrdn, SymBadBrdn");
625: }
627: del_buf = PetscAbsReal(PetscRealPart(delta));
628: del_buf = PetscMin(del_buf, del_max);
629: del_buf = PetscMax(del_buf, del_min);
630: if (is_dbfgs || is_ddfp || is_dqn) {
631: Mat_DQN *lqn = (Mat_DQN *)lmvm->ctx;
632: Mat_LMVM *dbase = (Mat_LMVM *)lqn->diag_qn->data;
633: Mat_DiagBrdn *diagctx = (Mat_DiagBrdn *)dbase->ctx;
635: diagctx->delta = del_buf;
636: } else {
637: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
639: lsb = (Mat_SymBrdn *)lmvm->ctx;
640: lsb->delta = del_buf;
641: }
642: PetscFunctionReturn(PETSC_SUCCESS);
643: }
645: /*@
646: MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.
648: Input Parameters:
649: + B - the `MATLMVM` matrix
650: - stype - scale type, see `MatLMVMSymBroydenScaleType`
652: Options Database Key:
653: . -mat_lmvm_scale_type <none,scalar,diagonal> - set the scaling type
655: Level: intermediate
657: MatLMVMSymBrdnScaleTypes\:
658: + `MAT_LMVM_SYMBROYDEN_SCALE_NONE` - initial Hessian is the identity matrix
659: . `MAT_LMVM_SYMBROYDEN_SCALE_SCALAR` - use the Shanno scalar as the initial Hessian
660: - `MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL` - use a diagonalized BFGS update as the initial Hessian
662: .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`, `MatCreateLMVMSymBroyden()`, `MatLMVMSymBroydenScaleType`
663: @*/
664: PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
665: {
666: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
667: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
669: PetscFunctionBegin;
671: lsb->scale_type = stype;
672: PetscFunctionReturn(PETSC_SUCCESS);
673: }
675: /*@
676: MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
677: for approximating Jacobians.
679: Collective
681: Input Parameters:
682: + comm - MPI communicator, set to `PETSC_COMM_SELF`
683: . n - number of local rows for storage vectors
684: - N - global size of the storage vectors
686: Output Parameter:
687: . B - the matrix
689: Options Database Keys:
690: + -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
691: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
692: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
693: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
694: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
695: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
696: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
698: Level: intermediate
700: Notes:
701: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
702: paradigm instead of this routine directly.
704: L-SymBrdn is a convex combination of L-DFP and
705: L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor
706: phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed
707: to be symmetric positive-definite.
709: To use the L-SymBrdn matrix with other vector types, the matrix must be
710: created using MatCreate() and MatSetType(), followed by `MatLMVMAllocate()`.
711: This ensures that the internal storage and work vectors are duplicated from the
712: correct type of vector.
714: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
715: `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`
716: @*/
717: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
718: {
719: PetscFunctionBegin;
720: PetscCall(KSPInitializePackage());
721: PetscCall(MatCreate(comm, B));
722: PetscCall(MatSetSizes(*B, n, n, N, N));
723: PetscCall(MatSetType(*B, MATLMVMSYMBROYDEN));
724: PetscCall(MatSetUp(*B));
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
729: {
730: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
731: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
733: PetscFunctionBegin;
734: if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
735: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
736: PetscCall(MatLMVMApplyJ0Fwd(B, X, Z));
737: } else {
738: switch (lsb->scale_type) {
739: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
740: PetscCall(VecAXPBY(Z, 1.0 / lsb->sigma, 0.0, X));
741: break;
742: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
743: PetscCall(MatMult(lsb->D, X, Z));
744: break;
745: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
746: default:
747: PetscCall(VecCopy(X, Z));
748: break;
749: }
750: }
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
755: {
756: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
757: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
759: PetscFunctionBegin;
760: if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
761: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
762: PetscCall(MatLMVMApplyJ0Inv(B, F, dX));
763: } else {
764: switch (lsb->scale_type) {
765: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
766: PetscCall(VecAXPBY(dX, lsb->sigma, 0.0, F));
767: break;
768: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
769: PetscCall(MatSolve(lsb->D, F, dX));
770: break;
771: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
772: default:
773: PetscCall(VecCopy(F, dX));
774: break;
775: }
776: }
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
781: {
782: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
783: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
784: PetscInt i, start;
785: PetscReal a, b, c, sig1, sig2, signew;
787: PetscFunctionBegin;
788: if (lsb->sigma_hist == 0) {
789: signew = 1.0;
790: } else {
791: start = PetscMax(0, lmvm->k - lsb->sigma_hist + 1);
792: signew = 0.0;
793: if (lsb->alpha == 1.0) {
794: for (i = start; i <= lmvm->k; ++i) signew += lsb->yts[i] / lsb->yty[i];
795: } else if (lsb->alpha == 0.5) {
796: for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yty[i];
797: signew = PetscSqrtReal(signew);
798: } else if (lsb->alpha == 0.0) {
799: for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yts[i];
800: } else {
801: /* compute coefficients of the quadratic */
802: a = b = c = 0.0;
803: for (i = start; i <= lmvm->k; ++i) {
804: a += lsb->yty[i];
805: b += lsb->yts[i];
806: c += lsb->sts[i];
807: }
808: a *= lsb->alpha;
809: b *= -(2.0 * lsb->alpha - 1.0);
810: c *= lsb->alpha - 1.0;
811: /* use quadratic formula to find roots */
812: sig1 = (-b + PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
813: sig2 = (-b - PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
814: /* accept the positive root as the scalar */
815: if (sig1 > 0.0) {
816: signew = sig1;
817: } else if (sig2 > 0.0) {
818: signew = sig2;
819: } else {
820: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
821: }
822: }
823: }
824: lsb->sigma = lsb->rho * signew + (1.0 - lsb->rho) * lsb->sigma;
825: PetscFunctionReturn(PETSC_SUCCESS);
826: }