Actual source code: symbadbrdn.c
1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
4: /*------------------------------------------------------------*/
6: static PetscErrorCode MatSolve_LMVMSymBadBrdn(Mat B, Vec F, Vec dX)
7: {
8: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
9: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
10: PetscInt i, j;
11: PetscScalar yjtqi, sjtyi, wtyi, ytx, stf, wtf, ytq;
13: PetscFunctionBegin;
14: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
15: if (lsb->phi == 0.0) {
16: PetscCall(MatSolve_LMVMBFGS(B, F, dX));
17: PetscFunctionReturn(PETSC_SUCCESS);
18: }
19: if (lsb->phi == 1.0) {
20: PetscCall(MatSolve_LMVMDFP(B, F, dX));
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: VecCheckSameSize(F, 2, dX, 3);
25: VecCheckMatCompatible(B, dX, 3, F, 2);
27: if (lsb->needQ) {
28: /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
29: for (i = 0; i <= lmvm->k; ++i) {
30: PetscCall(MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]));
31: for (j = 0; j <= i - 1; ++j) {
32: /* Compute the necessary dot products */
33: PetscCall(VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi));
34: PetscCall(VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi));
35: PetscCall(VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi));
36: PetscCall(VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi));
37: /* Compute the pure DFP component of the inverse application*/
38: PetscCall(VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi) / lsb->ytq[j], PetscRealPart(sjtyi) / lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]));
39: /* Tack on the convexly scaled extras to the inverse application*/
40: if (lsb->psi[j] > 0.0) {
41: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]));
42: PetscCall(VecDot(lsb->work, lmvm->Y[i], &wtyi));
43: PetscCall(VecAXPY(lsb->Q[i], lsb->phi * lsb->ytq[j] * PetscRealPart(wtyi), lsb->work));
44: }
45: }
46: PetscCall(VecDot(lmvm->Y[i], lsb->Q[i], &ytq));
47: lsb->ytq[i] = PetscRealPart(ytq);
48: }
49: lsb->needQ = PETSC_FALSE;
50: }
52: /* Start the outer iterations for ((B^{-1}) * dX) */
53: PetscCall(MatSymBrdnApplyJ0Inv(B, F, dX));
54: for (i = 0; i <= lmvm->k; ++i) {
55: /* Compute the necessary dot products -- store yTs and yTp for inner iterations later */
56: PetscCall(VecDotBegin(lmvm->Y[i], dX, &ytx));
57: PetscCall(VecDotBegin(lmvm->S[i], F, &stf));
58: PetscCall(VecDotEnd(lmvm->Y[i], dX, &ytx));
59: PetscCall(VecDotEnd(lmvm->S[i], F, &stf));
60: /* Compute the pure DFP component */
61: PetscCall(VecAXPBYPCZ(dX, -PetscRealPart(ytx) / lsb->ytq[i], PetscRealPart(stf) / lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]));
62: /* Tack on the convexly scaled extras */
63: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]));
64: PetscCall(VecDot(lsb->work, F, &wtf));
65: PetscCall(VecAXPY(dX, lsb->phi * lsb->ytq[i] * PetscRealPart(wtf), lsb->work));
66: }
68: PetscFunctionReturn(PETSC_SUCCESS);
69: }
71: /*------------------------------------------------------------*/
73: static PetscErrorCode MatMult_LMVMSymBadBrdn(Mat B, Vec X, Vec Z)
74: {
75: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
76: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
77: PetscInt i, j;
78: PetscReal numer;
79: PetscScalar sjtpi, sjtyi, yjtsi, yjtqi, wtsi, wtyi, stz, ytx, ytq, wtx, stp;
81: PetscFunctionBegin;
82: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
83: if (lsb->phi == 0.0) {
84: PetscCall(MatMult_LMVMBFGS(B, X, Z));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
87: if (lsb->phi == 1.0) {
88: PetscCall(MatMult_LMVMDFP(B, X, Z));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: VecCheckSameSize(X, 2, Z, 3);
93: VecCheckMatCompatible(B, X, 2, Z, 3);
95: if (lsb->needQ) {
96: /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
97: for (i = 0; i <= lmvm->k; ++i) {
98: PetscCall(MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]));
99: for (j = 0; j <= i - 1; ++j) {
100: /* Compute the necessary dot products */
101: PetscCall(VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi));
102: PetscCall(VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi));
103: PetscCall(VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi));
104: PetscCall(VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi));
105: /* Compute the pure DFP component of the inverse application*/
106: PetscCall(VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi) / lsb->ytq[j], PetscRealPart(sjtyi) / lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]));
107: /* Tack on the convexly scaled extras to the inverse application*/
108: if (lsb->psi[j] > 0.0) {
109: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]));
110: PetscCall(VecDot(lsb->work, lmvm->Y[i], &wtyi));
111: PetscCall(VecAXPY(lsb->Q[i], lsb->phi * lsb->ytq[j] * PetscRealPart(wtyi), lsb->work));
112: }
113: }
114: PetscCall(VecDot(lmvm->Y[i], lsb->Q[i], &ytq));
115: lsb->ytq[i] = PetscRealPart(ytq);
116: }
117: lsb->needQ = PETSC_FALSE;
118: }
119: if (lsb->needP) {
120: /* Start the loop for (P[k] = (B_k) * S[k]) */
121: for (i = 0; i <= lmvm->k; ++i) {
122: PetscCall(MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]));
123: for (j = 0; j <= i - 1; ++j) {
124: /* Compute the necessary dot products */
125: PetscCall(VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi));
126: PetscCall(VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi));
127: PetscCall(VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi));
128: PetscCall(VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi));
129: /* Compute the pure BFGS component of the forward product */
130: PetscCall(VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi) / lsb->stp[j], PetscRealPart(yjtsi) / lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]));
131: /* Tack on the convexly scaled extras to the forward product */
132: if (lsb->phi > 0.0) {
133: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]));
134: PetscCall(VecDot(lsb->work, lmvm->S[i], &wtsi));
135: PetscCall(VecAXPY(lsb->P[i], lsb->psi[j] * lsb->stp[j] * PetscRealPart(wtsi), lsb->work));
136: }
137: }
138: PetscCall(VecDot(lmvm->S[i], lsb->P[i], &stp));
139: lsb->stp[i] = PetscRealPart(stp);
140: if (lsb->phi == 1.0) {
141: lsb->psi[i] = 0.0;
142: } else if (lsb->phi == 0.0) {
143: lsb->psi[i] = 1.0;
144: } else {
145: numer = (1.0 - lsb->phi) * lsb->yts[i] * lsb->yts[i];
146: lsb->psi[i] = numer / (numer + (lsb->phi * lsb->ytq[i] * lsb->stp[i]));
147: }
148: }
149: lsb->needP = PETSC_FALSE;
150: }
152: /* Start the outer iterations for (B * X) */
153: PetscCall(MatSymBrdnApplyJ0Fwd(B, X, Z));
154: for (i = 0; i <= lmvm->k; ++i) {
155: /* Compute the necessary dot products */
156: PetscCall(VecDotBegin(lmvm->S[i], Z, &stz));
157: PetscCall(VecDotBegin(lmvm->Y[i], X, &ytx));
158: PetscCall(VecDotEnd(lmvm->S[i], Z, &stz));
159: PetscCall(VecDotEnd(lmvm->Y[i], X, &ytx));
160: /* Compute the pure BFGS component */
161: PetscCall(VecAXPBYPCZ(Z, -PetscRealPart(stz) / lsb->stp[i], PetscRealPart(ytx) / lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]));
162: /* Tack on the convexly scaled extras */
163: PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]));
164: PetscCall(VecDot(lsb->work, X, &wtx));
165: PetscCall(VecAXPY(Z, lsb->psi[i] * lsb->stp[i] * PetscRealPart(wtx), lsb->work));
166: }
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: /*------------------------------------------------------------*/
172: static PetscErrorCode MatSetFromOptions_LMVMSymBadBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
173: {
174: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
175: Mat_SymBrdn *lsb = (Mat_SymBrdn *)lmvm->ctx;
176: Mat_LMVM *dbase;
177: Mat_DiagBrdn *dctx;
179: PetscFunctionBegin;
180: PetscCall(MatSetFromOptions_LMVMSymBrdn(B, PetscOptionsObject));
181: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
182: dbase = (Mat_LMVM *)lsb->D->data;
183: dctx = (Mat_DiagBrdn *)dbase->ctx;
184: dctx->forward = PETSC_FALSE;
185: }
186: PetscFunctionReturn(PETSC_SUCCESS);
187: }
189: /*------------------------------------------------------------*/
191: PetscErrorCode MatCreate_LMVMSymBadBrdn(Mat B)
192: {
193: Mat_LMVM *lmvm;
195: PetscFunctionBegin;
196: PetscCall(MatCreate_LMVMSymBrdn(B));
197: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBADBROYDEN));
198: B->ops->setfromoptions = MatSetFromOptions_LMVMSymBadBrdn;
199: B->ops->solve = MatSolve_LMVMSymBadBrdn;
201: lmvm = (Mat_LMVM *)B->data;
202: lmvm->ops->mult = MatMult_LMVMSymBadBrdn;
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*------------------------------------------------------------*/
208: /*@
209: MatCreateLMVMSymBadBroyden - Creates a limited-memory Symmetric "Bad" Broyden-type matrix used
210: for approximating Jacobians. L-SymBadBrdn is a convex combination of L-DFP and
211: L-BFGS such that SymBadBrdn^{-1} = (1 - phi)*BFGS^{-1} + phi*DFP^{-1}. The combination factor
212: phi is restricted to the range [0, 1], where the L-SymBadBrdn matrix is guaranteed
213: to be symmetric positive-definite. Note that this combination is on the inverses and not
214: on the forwards. For forward convex combinations, use the L-SymBrdn matrix.
216: To use the L-SymBrdn matrix with other vector types, the matrix must be
217: created using MatCreate() and MatSetType(), followed by `MatLMVMAllocate()`.
218: This ensures that the internal storage and work vectors are duplicated from the
219: correct type of vector.
221: Collective
223: Input Parameters:
224: + comm - MPI communicator
225: . n - number of local rows for storage vectors
226: - N - global size of the storage vectors
228: Output Parameter:
229: . B - the matrix
231: Options Database Keys:
232: + -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
233: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
234: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
235: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
236: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
237: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
238: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
240: Level: intermediate
242: Note:
243: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
244: paradigm instead of this routine directly.
246: .seealso: [](ch_ksp), [LMVM Matrices](sec_matlmvm), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
247: `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`
248: @*/
249: PetscErrorCode MatCreateLMVMSymBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
250: {
251: PetscFunctionBegin;
252: PetscCall(MatCreate(comm, B));
253: PetscCall(MatSetSizes(*B, n, n, N, N));
254: PetscCall(MatSetType(*B, MATLMVMSYMBADBROYDEN));
255: PetscCall(MatSetUp(*B));
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }