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