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