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: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
41: if (lsb->phi == 0.0) {
42: MatSolve_LMVMBFGS(B, F, dX);
43: return 0;
44: }
45: if (lsb->phi == 1.0) {
46: MatSolve_LMVMDFP(B, F, dX);
47: return 0;
48: }
50: VecCheckSameSize(F, 2, dX, 3);
51: VecCheckMatCompatible(B, dX, 3, F, 2);
53: if (lsb->needP) {
54: /* Start the loop for (P[k] = (B_k) * S[k]) */
55: for (i = 0; i <= lmvm->k; ++i) {
56: MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);
57: for (j = 0; j <= i - 1; ++j) {
58: /* Compute the necessary dot products */
59: VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);
60: VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
61: VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);
62: VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
63: /* Compute the pure BFGS component of the forward product */
64: 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: VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);
68: VecDot(lsb->work, lmvm->S[i], &wtsi);
69: VecAXPY(lsb->P[i], lsb->phi * lsb->stp[j] * PetscRealPart(wtsi), lsb->work);
70: }
71: }
72: 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: MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]);
81: for (j = 0; j <= i - 1; ++j) {
82: /* Compute the necessary dot products */
83: VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi);
84: VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
85: VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi);
86: VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
87: /* Compute the pure DFP component of the inverse application*/
88: 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: VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]);
92: VecDot(lsb->work, lmvm->Y[i], &wtyi);
93: VecAXPY(lsb->Q[i], lsb->psi[j] * lsb->ytq[j] * PetscRealPart(wtyi), lsb->work);
94: }
95: }
96: 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: MatSymBrdnApplyJ0Inv(B, F, dX);
112: for (i = 0; i <= lmvm->k; ++i) {
113: /* Compute the necessary dot products -- store yTs and yTp for inner iterations later */
114: VecDotBegin(lmvm->Y[i], dX, &ytx);
115: VecDotBegin(lmvm->S[i], F, &stf);
116: VecDotEnd(lmvm->Y[i], dX, &ytx);
117: VecDotEnd(lmvm->S[i], F, &stf);
118: /* Compute the pure DFP component */
119: VecAXPBYPCZ(dX, -PetscRealPart(ytx) / lsb->ytq[i], PetscRealPart(stf) / lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]);
120: /* Tack on the convexly scaled extras */
121: VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]);
122: VecDot(lsb->work, F, &wtf);
123: VecAXPY(dX, lsb->psi[i] * lsb->ytq[i] * PetscRealPart(wtf), lsb->work);
124: }
126: return 0;
127: }
129: /*------------------------------------------------------------*/
131: /*
132: The forward-product below is the matrix-free implementation of
133: Equation 16 in Dennis and Wolkowicz "Sizing and Least Change Secant
134: Methods" (http://www.caam.rice.edu/caam/trs/90/TR90-05.pdf).
136: P[i] = (B_i)*S[i] terms are computed ahead of time whenever
137: the matrix is updated with a new (S[i], Y[i]) pair. This allows
138: repeated calls of MatMult inside KSP solvers without unnecessarily
139: recomputing P[i] terms in expensive nested-loops.
141: Z <- J0 * X
143: for i=0,1,2,...,k
144: # P[i] = (B_k) * S[i]
146: rho = 1.0 / (Y[i]^T S[i])
147: alpha = rho * (Y[i]^T F)
148: zeta = 1.0 / (S[i]^T P[i])
149: gamma = zeta * (S[i]^T dX)
151: dX <- dX - (gamma * P[i]) + (alpha * S[i])
152: W <- (rho * Y[i]) - (zeta * P[i])
153: dX <- dX + (phi * (S[i]^T P[i]) * (W^T F) * W)
154: end
155: */
156: static PetscErrorCode MatMult_LMVMSymBrdn(Mat B, Vec X, Vec Z)
157: {
158: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
159: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
160: PetscInt i, j;
161: PetscScalar sjtpi, yjtsi, wtsi, stz, ytx, wtx, stp;
163: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
164: if (lsb->phi == 0.0) {
165: MatMult_LMVMBFGS(B, X, Z);
166: return 0;
167: }
168: if (lsb->phi == 1.0) {
169: MatMult_LMVMDFP(B, X, Z);
170: return 0;
171: }
173: VecCheckSameSize(X, 2, Z, 3);
174: VecCheckMatCompatible(B, X, 2, Z, 3);
176: if (lsb->needP) {
177: /* Start the loop for (P[k] = (B_k) * S[k]) */
178: for (i = 0; i <= lmvm->k; ++i) {
179: MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);
180: for (j = 0; j <= i - 1; ++j) {
181: /* Compute the necessary dot products */
182: VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);
183: VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
184: VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);
185: VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
186: /* Compute the pure BFGS component of the forward product */
187: VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi) / lsb->stp[j], PetscRealPart(yjtsi) / lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);
188: /* Tack on the convexly scaled extras to the forward product */
189: if (lsb->phi > 0.0) {
190: VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);
191: VecDot(lsb->work, lmvm->S[i], &wtsi);
192: VecAXPY(lsb->P[i], lsb->phi * lsb->stp[j] * PetscRealPart(wtsi), lsb->work);
193: }
194: }
195: VecDot(lmvm->S[i], lsb->P[i], &stp);
196: lsb->stp[i] = PetscRealPart(stp);
197: }
198: lsb->needP = PETSC_FALSE;
199: }
201: /* Start the outer iterations for (B * X) */
202: MatSymBrdnApplyJ0Fwd(B, X, Z);
203: for (i = 0; i <= lmvm->k; ++i) {
204: /* Compute the necessary dot products */
205: VecDotBegin(lmvm->S[i], Z, &stz);
206: VecDotBegin(lmvm->Y[i], X, &ytx);
207: VecDotEnd(lmvm->S[i], Z, &stz);
208: VecDotEnd(lmvm->Y[i], X, &ytx);
209: /* Compute the pure BFGS component */
210: VecAXPBYPCZ(Z, -PetscRealPart(stz) / lsb->stp[i], PetscRealPart(ytx) / lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]);
211: /* Tack on the convexly scaled extras */
212: VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]);
213: VecDot(lsb->work, X, &wtx);
214: VecAXPY(Z, lsb->phi * lsb->stp[i] * PetscRealPart(wtx), lsb->work);
215: }
216: return 0;
217: }
219: /*------------------------------------------------------------*/
221: static PetscErrorCode MatUpdate_LMVMSymBrdn(Mat B, Vec X, Vec F)
222: {
223: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
224: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
225: Mat_LMVM *dbase;
226: Mat_DiagBrdn *dctx;
227: PetscInt old_k, i;
228: PetscReal curvtol;
229: PetscScalar curvature, ytytmp, ststmp;
231: if (!lmvm->m) return 0;
232: if (lmvm->prev_set) {
233: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
234: VecAYPX(lmvm->Xprev, -1.0, X);
235: VecAYPX(lmvm->Fprev, -1.0, F);
236: /* Test if the updates can be accepted */
237: VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
238: VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
239: VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
240: VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
241: if (PetscRealPart(ststmp) < lmvm->eps) {
242: curvtol = 0.0;
243: } else {
244: curvtol = lmvm->eps * PetscRealPart(ststmp);
245: }
246: if (PetscRealPart(curvature) > curvtol) {
247: /* Update is good, accept it */
248: lsb->watchdog = 0;
249: lsb->needP = lsb->needQ = PETSC_TRUE;
250: old_k = lmvm->k;
251: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
252: /* If we hit the memory limit, shift the yts, yty and sts arrays */
253: if (old_k == lmvm->k) {
254: for (i = 0; i <= lmvm->k - 1; ++i) {
255: lsb->yts[i] = lsb->yts[i + 1];
256: lsb->yty[i] = lsb->yty[i + 1];
257: lsb->sts[i] = lsb->sts[i + 1];
258: }
259: }
260: /* Update history of useful scalars */
261: VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);
262: lsb->yts[lmvm->k] = PetscRealPart(curvature);
263: lsb->yty[lmvm->k] = PetscRealPart(ytytmp);
264: lsb->sts[lmvm->k] = PetscRealPart(ststmp);
265: /* Compute the scalar scale if necessary */
266: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) MatSymBrdnComputeJ0Scalar(B);
267: } else {
268: /* Update is bad, skip it */
269: ++lmvm->nrejects;
270: ++lsb->watchdog;
271: }
272: } else {
273: switch (lsb->scale_type) {
274: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
275: dbase = (Mat_LMVM *)lsb->D->data;
276: dctx = (Mat_DiagBrdn *)dbase->ctx;
277: VecSet(dctx->invD, lsb->delta);
278: break;
279: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
280: lsb->sigma = lsb->delta;
281: break;
282: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
283: lsb->sigma = 1.0;
284: break;
285: default:
286: break;
287: }
288: }
290: /* Update the scaling */
291: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) MatLMVMUpdate(lsb->D, X, F);
293: if (lsb->watchdog > lsb->max_seq_rejects) {
294: MatLMVMReset(B, PETSC_FALSE);
295: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) MatLMVMReset(lsb->D, PETSC_FALSE);
296: }
298: /* Save the solution and function to be used in the next update */
299: VecCopy(X, lmvm->Xprev);
300: VecCopy(F, lmvm->Fprev);
301: lmvm->prev_set = PETSC_TRUE;
302: return 0;
303: }
305: /*------------------------------------------------------------*/
307: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
308: {
309: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
310: Mat_SymBrdn *blsb = (Mat_SymBrdn *)bdata->ctx;
311: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
312: Mat_SymBrdn *mlsb = (Mat_SymBrdn *)mdata->ctx;
313: PetscInt i;
315: mlsb->phi = blsb->phi;
316: mlsb->needP = blsb->needP;
317: mlsb->needQ = blsb->needQ;
318: for (i = 0; i <= bdata->k; ++i) {
319: mlsb->stp[i] = blsb->stp[i];
320: mlsb->ytq[i] = blsb->ytq[i];
321: mlsb->yts[i] = blsb->yts[i];
322: mlsb->psi[i] = blsb->psi[i];
323: VecCopy(blsb->P[i], mlsb->P[i]);
324: VecCopy(blsb->Q[i], mlsb->Q[i]);
325: }
326: mlsb->scale_type = blsb->scale_type;
327: mlsb->alpha = blsb->alpha;
328: mlsb->beta = blsb->beta;
329: mlsb->rho = blsb->rho;
330: mlsb->delta = blsb->delta;
331: mlsb->sigma_hist = blsb->sigma_hist;
332: mlsb->watchdog = blsb->watchdog;
333: mlsb->max_seq_rejects = blsb->max_seq_rejects;
334: switch (blsb->scale_type) {
335: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
336: mlsb->sigma = blsb->sigma;
337: break;
338: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
339: MatCopy(blsb->D, mlsb->D, SAME_NONZERO_PATTERN);
340: break;
341: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
342: mlsb->sigma = 1.0;
343: break;
344: default:
345: break;
346: }
347: return 0;
348: }
350: /*------------------------------------------------------------*/
352: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
353: {
354: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
355: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
356: Mat_LMVM *dbase;
357: Mat_DiagBrdn *dctx;
359: lsb->watchdog = 0;
360: lsb->needP = lsb->needQ = PETSC_TRUE;
361: if (lsb->allocated) {
362: if (destructive) {
363: VecDestroy(&lsb->work);
364: PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);
365: PetscFree(lsb->psi);
366: VecDestroyVecs(lmvm->m, &lsb->P);
367: VecDestroyVecs(lmvm->m, &lsb->Q);
368: switch (lsb->scale_type) {
369: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
370: MatLMVMReset(lsb->D, PETSC_TRUE);
371: break;
372: default:
373: break;
374: }
375: lsb->allocated = PETSC_FALSE;
376: } else {
377: PetscMemzero(lsb->psi, lmvm->m);
378: switch (lsb->scale_type) {
379: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
380: lsb->sigma = lsb->delta;
381: break;
382: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
383: MatLMVMReset(lsb->D, PETSC_FALSE);
384: dbase = (Mat_LMVM *)lsb->D->data;
385: dctx = (Mat_DiagBrdn *)dbase->ctx;
386: VecSet(dctx->invD, lsb->delta);
387: break;
388: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
389: lsb->sigma = 1.0;
390: break;
391: default:
392: break;
393: }
394: }
395: }
396: MatReset_LMVM(B, destructive);
397: return 0;
398: }
400: /*------------------------------------------------------------*/
402: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
403: {
404: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
405: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
407: MatAllocate_LMVM(B, X, F);
408: if (!lsb->allocated) {
409: VecDuplicate(X, &lsb->work);
410: if (lmvm->m > 0) {
411: PetscMalloc5(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts);
412: PetscCalloc1(lmvm->m, &lsb->psi);
413: VecDuplicateVecs(X, lmvm->m, &lsb->P);
414: VecDuplicateVecs(X, lmvm->m, &lsb->Q);
415: }
416: switch (lsb->scale_type) {
417: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
418: MatLMVMAllocate(lsb->D, X, F);
419: break;
420: default:
421: break;
422: }
423: lsb->allocated = PETSC_TRUE;
424: }
425: return 0;
426: }
428: /*------------------------------------------------------------*/
430: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
431: {
432: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
433: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
435: if (lsb->allocated) {
436: VecDestroy(&lsb->work);
437: PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);
438: PetscFree(lsb->psi);
439: VecDestroyVecs(lmvm->m, &lsb->P);
440: VecDestroyVecs(lmvm->m, &lsb->Q);
441: lsb->allocated = PETSC_FALSE;
442: }
443: MatDestroy(&lsb->D);
444: PetscFree(lmvm->ctx);
445: MatDestroy_LMVM(B);
446: return 0;
447: }
449: /*------------------------------------------------------------*/
451: static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
452: {
453: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
454: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
455: PetscInt n, N;
457: MatSetUp_LMVM(B);
458: if (!lsb->allocated) {
459: VecDuplicate(lmvm->Xprev, &lsb->work);
460: if (lmvm->m > 0) {
461: PetscMalloc5(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts);
462: PetscCalloc1(lmvm->m, &lsb->psi);
463: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->P);
464: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->Q);
465: }
466: switch (lsb->scale_type) {
467: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
468: MatGetLocalSize(B, &n, &n);
469: MatGetSize(B, &N, &N);
470: MatSetSizes(lsb->D, n, n, N, N);
471: MatSetUp(lsb->D);
472: break;
473: default:
474: break;
475: }
476: lsb->allocated = PETSC_TRUE;
477: }
478: return 0;
479: }
481: /*------------------------------------------------------------*/
483: PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
484: {
485: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
486: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
487: PetscBool isascii;
489: PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii);
490: if (isascii) {
491: PetscViewerASCIIPrintf(pv, "Scale type: %s\n", MatLMVMSymBroydenScaleTypes[lsb->scale_type]);
492: PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", lsb->sigma_hist);
493: PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)lsb->alpha, (double)lsb->beta, (double)lsb->rho);
494: PetscViewerASCIIPrintf(pv, "Convex factors: phi=%g, theta=%g\n", (double)lsb->phi, (double)lsb->theta);
495: }
496: MatView_LMVM(B, pv);
497: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) MatView(lsb->D, pv);
498: return 0;
499: }
501: /*------------------------------------------------------------*/
503: PetscErrorCode MatSetFromOptions_LMVMSymBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
504: {
505: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
506: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
508: MatSetFromOptions_LMVM(B, PetscOptionsObject);
509: PetscOptionsHeadBegin(PetscOptionsObject, "Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
510: PetscOptionsReal("-mat_lmvm_phi", "(developer) convex ratio between BFGS and DFP components of the update", "", lsb->phi, &lsb->phi, NULL);
512: MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject);
513: PetscOptionsHeadEnd();
514: return 0;
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: PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", lsb->beta, &lsb->beta, NULL);
527: PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", lsb->theta, &lsb->theta, NULL);
529: PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", lsb->rho, &lsb->rho, NULL);
531: PetscOptionsReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", lsb->alpha, &lsb->alpha, NULL);
533: 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);
534: PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBrdnScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg);
535: if (flg) MatLMVMSymBroydenSetScaleType(B, stype);
536: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
537: const char *prefix;
539: MatGetOptionsPrefix(B, &prefix);
540: MatSetOptionsPrefix(lsb->D, prefix);
541: MatAppendOptionsPrefix(lsb->D, "J0_");
542: MatSetFromOptions(lsb->D);
543: dbase = (Mat_LMVM *)lsb->D->data;
544: dctx = (Mat_DiagBrdn *)dbase->ctx;
545: dctx->delta_min = lsb->delta_min;
546: dctx->delta_max = lsb->delta_max;
547: dctx->theta = lsb->theta;
548: dctx->rho = lsb->rho;
549: dctx->alpha = lsb->alpha;
550: dctx->beta = lsb->beta;
551: dctx->sigma_hist = lsb->sigma_hist;
552: }
553: return 0;
554: }
556: /*------------------------------------------------------------*/
558: PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
559: {
560: Mat_LMVM *lmvm;
561: Mat_SymBrdn *lsb;
563: MatCreate_LMVM(B);
564: PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN);
565: MatSetOption(B, MAT_SPD, PETSC_TRUE);
566: MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE);
567: B->ops->view = MatView_LMVMSymBrdn;
568: B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
569: B->ops->setup = MatSetUp_LMVMSymBrdn;
570: B->ops->destroy = MatDestroy_LMVMSymBrdn;
571: B->ops->solve = MatSolve_LMVMSymBrdn;
573: lmvm = (Mat_LMVM *)B->data;
574: lmvm->square = PETSC_TRUE;
575: lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
576: lmvm->ops->reset = MatReset_LMVMSymBrdn;
577: lmvm->ops->update = MatUpdate_LMVMSymBrdn;
578: lmvm->ops->mult = MatMult_LMVMSymBrdn;
579: lmvm->ops->copy = MatCopy_LMVMSymBrdn;
581: PetscNew(&lsb);
582: lmvm->ctx = (void *)lsb;
583: lsb->allocated = PETSC_FALSE;
584: lsb->needP = lsb->needQ = PETSC_TRUE;
585: lsb->phi = 0.125;
586: lsb->theta = 0.125;
587: lsb->alpha = 1.0;
588: lsb->rho = 1.0;
589: lsb->beta = 0.5;
590: lsb->sigma = 1.0;
591: lsb->delta = 1.0;
592: lsb->delta_min = 1e-7;
593: lsb->delta_max = 100.0;
594: lsb->sigma_hist = 1;
595: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
596: lsb->watchdog = 0;
597: lsb->max_seq_rejects = lmvm->m / 2;
599: MatCreate(PetscObjectComm((PetscObject)B), &lsb->D);
600: MatSetType(lsb->D, MATLMVMDIAGBROYDEN);
601: return 0;
602: }
604: /*------------------------------------------------------------*/
606: /*@
607: MatLMVMSymBroydenSetDelta - Sets the starting value for the diagonal scaling vector computed
608: in the SymBrdn approximations (also works for BFGS and DFP).
610: Input Parameters:
611: + B - LMVM matrix
612: - delta - initial value for diagonal scaling
614: Level: intermediate
616: @*/
618: PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
619: {
620: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
621: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
622: PetscBool is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn;
624: PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);
625: PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);
626: PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn);
627: PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn);
629: lsb->delta = PetscAbsReal(PetscRealPart(delta));
630: lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
631: lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
632: return 0;
633: }
635: /*------------------------------------------------------------*/
637: /*@
638: MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.
640: Input Parameters:
641: + snes - the iterative context
642: - rtype - restart type
644: Options Database Key:
645: . -mat_lmvm_scale_type <none,scalar,diagonal> - set the scaling type
647: Level: intermediate
649: MatLMVMSymBrdnScaleTypes:
650: + MAT_LMVM_SYMBROYDEN_SCALE_NONE - initial Hessian is the identity matrix
651: . MAT_LMVM_SYMBROYDEN_SCALE_SCALAR - use the Shanno scalar as the initial Hessian
652: - MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL - use a diagonalized BFGS update as the initial Hessian
654: .seealso: [](chapter_ksp), `MATLMVMSYMBROYDEN`, `MatCreateLMVMSymBroyden()`
655: @*/
656: PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
657: {
658: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
659: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
662: lsb->scale_type = stype;
663: return 0;
664: }
666: /*------------------------------------------------------------*/
668: /*@
669: MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
670: for approximating Jacobians. L-SymBrdn is a convex combination of L-DFP and
671: L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor
672: phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed
673: to be symmetric positive-definite.
675: The provided local and global sizes must match the solution and function vectors
676: used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have
677: storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
678: parallel. 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: Collective
685: Input Parameters:
686: + comm - MPI communicator, set to PETSC_COMM_SELF
687: . n - number of local rows for storage vectors
688: - N - global size of the storage vectors
690: Output Parameter:
691: . B - the matrix
693: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
694: paradigm instead of this routine directly.
696: Options Database Keys:
697: + -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
698: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
699: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
700: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
701: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
702: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
703: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
705: Level: intermediate
707: .seealso: [](chapter_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
708: `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`
709: @*/
710: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
711: {
712: MatCreate(comm, B);
713: MatSetSizes(*B, n, n, N, N);
714: MatSetType(*B, MATLMVMSYMBROYDEN);
715: MatSetUp(*B);
716: return 0;
717: }
719: /*------------------------------------------------------------*/
721: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
722: {
723: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
724: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
726: if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
727: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
728: MatLMVMApplyJ0Fwd(B, X, Z);
729: } else {
730: switch (lsb->scale_type) {
731: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
732: VecCopy(X, Z);
733: VecScale(Z, 1.0 / lsb->sigma);
734: break;
735: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
736: MatMult(lsb->D, X, Z);
737: break;
738: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
739: default:
740: VecCopy(X, Z);
741: break;
742: }
743: }
744: return 0;
745: }
747: /*------------------------------------------------------------*/
749: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
750: {
751: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
752: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
754: if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
755: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
756: MatLMVMApplyJ0Inv(B, F, dX);
757: } else {
758: switch (lsb->scale_type) {
759: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
760: VecCopy(F, dX);
761: VecScale(dX, lsb->sigma);
762: break;
763: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
764: MatSolve(lsb->D, F, dX);
765: break;
766: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
767: default:
768: VecCopy(F, dX);
769: break;
770: }
771: }
772: return 0;
773: }
775: /*------------------------------------------------------------*/
777: PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
778: {
779: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
780: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
781: PetscInt i, start;
782: PetscReal a, b, c, sig1, sig2, signew;
784: if (lsb->sigma_hist == 0) {
785: signew = 1.0;
786: } else {
787: start = PetscMax(0, lmvm->k - lsb->sigma_hist + 1);
788: signew = 0.0;
789: if (lsb->alpha == 1.0) {
790: for (i = start; i <= lmvm->k; ++i) signew += lsb->yts[i] / lsb->yty[i];
791: } else if (lsb->alpha == 0.5) {
792: for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yty[i];
793: signew = PetscSqrtReal(signew);
794: } else if (lsb->alpha == 0.0) {
795: for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yts[i];
796: } else {
797: /* compute coefficients of the quadratic */
798: a = b = c = 0.0;
799: for (i = start; i <= lmvm->k; ++i) {
800: a += lsb->yty[i];
801: b += lsb->yts[i];
802: c += lsb->sts[i];
803: }
804: a *= lsb->alpha;
805: b *= -(2.0 * lsb->alpha - 1.0);
806: c *= lsb->alpha - 1.0;
807: /* use quadratic formula to find roots */
808: sig1 = (-b + PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
809: sig2 = (-b - PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
810: /* accept the positive root as the scalar */
811: if (sig1 > 0.0) {
812: signew = sig1;
813: } else if (sig2 > 0.0) {
814: signew = sig2;
815: } else {
816: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
817: }
818: }
819: }
820: lsb->sigma = lsb->rho * signew + (1.0 - lsb->rho) * lsb->sigma;
821: return 0;
822: }