Actual source code: badbrdn.c

  1: #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h>

  3: /*------------------------------------------------------------*/

  5: /*
  6:   The solution method is the matrix-free implementation of the inverse Hessian in
  7:   Equation 6 on page 312 of Griewank "Broyden Updating, The Good and The Bad!"
  8:   (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).

 10:   Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
 11:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
 12:   repeated calls of MatSolve without incurring redundant computation.

 14:   dX <- J0^{-1} * F

 16:   for i=0,1,2,...,k
 17:     # Q[i] = (B_i)^{-1} * Y[i]
 18:     tau = (Y[i]^T F) / (Y[i]^T Y[i])
 19:     dX <- dX + (tau * (S[i] - Q[i]))
 20:   end
 21:  */

 23: static PetscErrorCode MatSolve_LMVMBadBrdn(Mat B, Vec F, Vec dX)
 24: {
 25:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
 26:   Mat_Brdn   *lbb  = (Mat_Brdn *)lmvm->ctx;
 27:   PetscInt    i, j;
 28:   PetscScalar yjtyi, ytf;

 30:   VecCheckSameSize(F, 2, dX, 3);
 31:   VecCheckMatCompatible(B, dX, 3, F, 2);

 33:   if (lbb->needQ) {
 34:     /* Pre-compute (Q[i] = (B_i)^{-1} * Y[i]) */
 35:     for (i = 0; i <= lmvm->k; ++i) {
 36:       MatLMVMApplyJ0Inv(B, lmvm->Y[i], lbb->Q[i]);
 37:       for (j = 0; j <= i - 1; ++j) {
 38:         VecDot(lmvm->Y[j], lmvm->Y[i], &yjtyi);
 39:         VecAXPBYPCZ(lbb->Q[i], PetscRealPart(yjtyi) / lbb->yty[j], -PetscRealPart(yjtyi) / lbb->yty[j], 1.0, lmvm->S[j], lbb->Q[j]);
 40:       }
 41:     }
 42:     lbb->needQ = PETSC_FALSE;
 43:   }

 45:   MatLMVMApplyJ0Inv(B, F, dX);
 46:   for (i = 0; i <= lmvm->k; ++i) {
 47:     VecDot(lmvm->Y[i], F, &ytf);
 48:     VecAXPBYPCZ(dX, PetscRealPart(ytf) / lbb->yty[i], -PetscRealPart(ytf) / lbb->yty[i], 1.0, lmvm->S[i], lbb->Q[i]);
 49:   }
 50:   return 0;
 51: }

 53: /*------------------------------------------------------------*/

 55: /*
 56:   The forward product is the matrix-free implementation of the direct update in
 57:   Equation 6 on page 302 of Griewank "Broyden Updating, The Good and The Bad!"
 58:   (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).

 60:   P[i] = (B_i)*S[i] terms are computed ahead of time whenever
 61:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
 62:   repeated calls of MatMult inside KSP solvers without unnecessarily
 63:   recomputing P[i] terms in expensive nested-loops.

 65:   Z <- J0 * X

 67:   for i=0,1,2,...,k
 68:     # P[i] = B_i * S[i]
 69:     tau = (Y[i]^T X) / (Y[i]^T S[i])
 70:     dX <- dX + (tau * (Y[i] - P[i]))
 71:   end
 72:  */

 74: static PetscErrorCode MatMult_LMVMBadBrdn(Mat B, Vec X, Vec Z)
 75: {
 76:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
 77:   Mat_Brdn   *lbb  = (Mat_Brdn *)lmvm->ctx;
 78:   PetscInt    i, j;
 79:   PetscScalar yjtsi, ytx;

 81:   VecCheckSameSize(X, 2, Z, 3);
 82:   VecCheckMatCompatible(B, X, 2, Z, 3);

 84:   if (lbb->needP) {
 85:     /* Pre-compute (P[i] = (B_i) * S[i]) */
 86:     for (i = 0; i <= lmvm->k; ++i) {
 87:       MatLMVMApplyJ0Fwd(B, lmvm->S[i], lbb->P[i]);
 88:       for (j = 0; j <= i - 1; ++j) {
 89:         VecDot(lmvm->Y[j], lmvm->S[i], &yjtsi);
 90:         VecAXPBYPCZ(lbb->P[i], PetscRealPart(yjtsi) / lbb->yts[j], -PetscRealPart(yjtsi) / lbb->yts[j], 1.0, lmvm->Y[j], lbb->P[j]);
 91:       }
 92:     }
 93:     lbb->needP = PETSC_FALSE;
 94:   }

 96:   MatLMVMApplyJ0Fwd(B, X, Z);
 97:   for (i = 0; i <= lmvm->k; ++i) {
 98:     VecDot(lmvm->Y[i], X, &ytx);
 99:     VecAXPBYPCZ(Z, PetscRealPart(ytx) / lbb->yts[i], -PetscRealPart(ytx) / lbb->yts[i], 1.0, lmvm->Y[i], lbb->P[i]);
100:   }
101:   return 0;
102: }

104: /*------------------------------------------------------------*/

106: static PetscErrorCode MatUpdate_LMVMBadBrdn(Mat B, Vec X, Vec F)
107: {
108:   Mat_LMVM   *lmvm = (Mat_LMVM *)B->data;
109:   Mat_Brdn   *lbb  = (Mat_Brdn *)lmvm->ctx;
110:   PetscInt    old_k, i;
111:   PetscScalar yty, yts;

113:   if (!lmvm->m) return 0;
114:   if (lmvm->prev_set) {
115:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
116:     VecAYPX(lmvm->Xprev, -1.0, X);
117:     VecAYPX(lmvm->Fprev, -1.0, F);
118:     /* Accept the update */
119:     lbb->needP = lbb->needQ = PETSC_TRUE;
120:     old_k                   = lmvm->k;
121:     MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
122:     /* If we hit the memory limit, shift the yty and yts arrays */
123:     if (old_k == lmvm->k) {
124:       for (i = 0; i <= lmvm->k - 1; ++i) {
125:         lbb->yty[i] = lbb->yty[i + 1];
126:         lbb->yts[i] = lbb->yts[i + 1];
127:       }
128:     }
129:     /* Accumulate the latest yTy and yTs dot products */
130:     VecDotBegin(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
131:     VecDotBegin(lmvm->Y[lmvm->k], lmvm->S[lmvm->k], &yts);
132:     VecDotEnd(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
133:     VecDotEnd(lmvm->Y[lmvm->k], lmvm->S[lmvm->k], &yts);
134:     lbb->yty[lmvm->k] = PetscRealPart(yty);
135:     lbb->yts[lmvm->k] = PetscRealPart(yts);
136:   }
137:   /* Save the solution and function to be used in the next update */
138:   VecCopy(X, lmvm->Xprev);
139:   VecCopy(F, lmvm->Fprev);
140:   lmvm->prev_set = PETSC_TRUE;
141:   return 0;
142: }

144: /*------------------------------------------------------------*/

146: static PetscErrorCode MatCopy_LMVMBadBrdn(Mat B, Mat M, MatStructure str)
147: {
148:   Mat_LMVM *bdata = (Mat_LMVM *)B->data;
149:   Mat_Brdn *bctx  = (Mat_Brdn *)bdata->ctx;
150:   Mat_LMVM *mdata = (Mat_LMVM *)M->data;
151:   Mat_Brdn *mctx  = (Mat_Brdn *)mdata->ctx;
152:   PetscInt  i;

154:   mctx->needP = bctx->needP;
155:   mctx->needQ = bctx->needQ;
156:   for (i = 0; i <= bdata->k; ++i) {
157:     mctx->yty[i] = bctx->yty[i];
158:     mctx->yts[i] = bctx->yts[i];
159:     VecCopy(bctx->P[i], mctx->P[i]);
160:     VecCopy(bctx->Q[i], mctx->Q[i]);
161:   }
162:   return 0;
163: }

165: /*------------------------------------------------------------*/

167: static PetscErrorCode MatReset_LMVMBadBrdn(Mat B, PetscBool destructive)
168: {
169:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
170:   Mat_Brdn *lbb  = (Mat_Brdn *)lmvm->ctx;

172:   lbb->needP = lbb->needQ = PETSC_TRUE;
173:   if (destructive && lbb->allocated) {
174:     PetscFree2(lbb->yty, lbb->yts);
175:     VecDestroyVecs(lmvm->m, &lbb->P);
176:     VecDestroyVecs(lmvm->m, &lbb->Q);
177:     lbb->allocated = PETSC_FALSE;
178:   }
179:   MatReset_LMVM(B, destructive);
180:   return 0;
181: }

183: /*------------------------------------------------------------*/

185: static PetscErrorCode MatAllocate_LMVMBadBrdn(Mat B, Vec X, Vec F)
186: {
187:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
188:   Mat_Brdn *lbb  = (Mat_Brdn *)lmvm->ctx;

190:   MatAllocate_LMVM(B, X, F);
191:   if (!lbb->allocated) {
192:     PetscMalloc2(lmvm->m, &lbb->yty, lmvm->m, &lbb->yts);
193:     if (lmvm->m > 0) {
194:       VecDuplicateVecs(X, lmvm->m, &lbb->P);
195:       VecDuplicateVecs(X, lmvm->m, &lbb->Q);
196:     }
197:     lbb->allocated = PETSC_TRUE;
198:   }
199:   return 0;
200: }

202: /*------------------------------------------------------------*/

204: static PetscErrorCode MatDestroy_LMVMBadBrdn(Mat B)
205: {
206:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
207:   Mat_Brdn *lbb  = (Mat_Brdn *)lmvm->ctx;

209:   if (lbb->allocated) {
210:     PetscFree2(lbb->yty, lbb->yts);
211:     VecDestroyVecs(lmvm->m, &lbb->P);
212:     VecDestroyVecs(lmvm->m, &lbb->Q);
213:     lbb->allocated = PETSC_FALSE;
214:   }
215:   PetscFree(lmvm->ctx);
216:   MatDestroy_LMVM(B);
217:   return 0;
218: }

220: /*------------------------------------------------------------*/

222: static PetscErrorCode MatSetUp_LMVMBadBrdn(Mat B)
223: {
224:   Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
225:   Mat_Brdn *lbb  = (Mat_Brdn *)lmvm->ctx;

227:   MatSetUp_LMVM(B);
228:   if (!lbb->allocated) {
229:     PetscMalloc2(lmvm->m, &lbb->yty, lmvm->m, &lbb->yts);
230:     if (lmvm->m > 0) {
231:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbb->P);
232:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbb->Q);
233:     }
234:     lbb->allocated = PETSC_TRUE;
235:   }
236:   return 0;
237: }

239: /*------------------------------------------------------------*/

241: PetscErrorCode MatCreate_LMVMBadBrdn(Mat B)
242: {
243:   Mat_LMVM *lmvm;
244:   Mat_Brdn *lbb;

246:   MatCreate_LMVM(B);
247:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMBADBROYDEN);
248:   B->ops->setup   = MatSetUp_LMVMBadBrdn;
249:   B->ops->destroy = MatDestroy_LMVMBadBrdn;
250:   B->ops->solve   = MatSolve_LMVMBadBrdn;

252:   lmvm                = (Mat_LMVM *)B->data;
253:   lmvm->square        = PETSC_TRUE;
254:   lmvm->ops->allocate = MatAllocate_LMVMBadBrdn;
255:   lmvm->ops->reset    = MatReset_LMVMBadBrdn;
256:   lmvm->ops->mult     = MatMult_LMVMBadBrdn;
257:   lmvm->ops->update   = MatUpdate_LMVMBadBrdn;
258:   lmvm->ops->copy     = MatCopy_LMVMBadBrdn;

260:   PetscNew(&lbb);
261:   lmvm->ctx      = (void *)lbb;
262:   lbb->allocated = PETSC_FALSE;
263:   lbb->needP = lbb->needQ = PETSC_TRUE;
264:   return 0;
265: }

267: /*------------------------------------------------------------*/

269: /*@
270:    MatCreateLMVMBadBroyden - Creates a limited-memory modified (aka "bad") Broyden-type
271:    approximation matrix used for a Jacobian. L-BadBrdn is not guaranteed to be
272:    symmetric or positive-definite.

274:    The provided local and global sizes must match the solution and function vectors
275:    used with `MatLMVMUpdate()` and `MatSolve()`. The resulting L-BadBrdn matrix will have
276:    storage vectors allocated with `VecCreateSeq()` in serial and `VecCreateMPI()` in
277:    parallel. To use the L-BadBrdn matrix with other vector types, the matrix must be
278:    created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
279:    This ensures that the internal storage and work vectors are duplicated from the
280:    correct type of vector.

282:    Collective

284:    Input Parameters:
285: +  comm - MPI communicator
286: .  n - number of local rows for storage vectors
287: -  N - global size of the storage vectors

289:    Output Parameter:
290: .  B - the matrix

292:    Options Database Keys:
293: +   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
294: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
295: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
296: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
297: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
298: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

300:    Level: intermediate

302:    Note:
303:    It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
304:    paradigm instead of this routine directly.

306: .seealso: [](chapter_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBADBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
307:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMSymBrdn()`
308: @*/
309: PetscErrorCode MatCreateLMVMBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
310: {
311:   MatCreate(comm, B);
312:   MatSetSizes(*B, n, n, N, N);
313:   MatSetType(*B, MATLMVMBADBROYDEN);
314:   MatSetUp(*B);
315:   return 0;
316: }