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