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