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