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: }