Actual source code: dfp.c

  1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
  2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>

  4: /*
  5:   Limited-memory Davidon-Fletcher-Powell method for approximating both
  6:   the forward product and inverse application of a Jacobian.
  7:  */

  9: /*------------------------------------------------------------*/

 11: /*
 12:   The solution method (approximate inverse Jacobian application) is
 13:   matrix-vector product version of the recursive formula given in
 14:   Equation (6.15) of Nocedal and Wright "Numerical Optimization" 2nd
 15:   edition, pg 139.

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

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

 23:   for i = 0,1,2,...,k
 24:     # Q[i] = (B_i)^{-1} * Y[i]
 25:     gamma = (S[i]^T F) / (Y[i]^T S[i])
 26:     zeta = (Y[i]^T dX) / (Y[i]^T Q[i])
 27:     dX <- dX + (gamma * S[i]) - (zeta * Q[i])
 28:   end
 29: */
 30: PetscErrorCode MatSolve_LMVMDFP(Mat B, Vec F, Vec dX)
 31: {
 32:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
 33:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
 34:   PetscInt     i, j;
 35:   PetscScalar  yjtqi, sjtyi, ytx, stf, ytq;

 39:   VecCheckSameSize(F, 2, dX, 3);
 40:   VecCheckMatCompatible(B, dX, 3, F, 2);

 42:   if (ldfp->needQ) {
 43:     /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
 44:     for (i = 0; i <= lmvm->k; ++i) {
 45:       MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], ldfp->Q[i]);
 46:       for (j = 0; j <= i - 1; ++j) {
 47:         /* Compute the necessary dot products */
 48:         VecDotBegin(lmvm->Y[j], ldfp->Q[i], &yjtqi);
 49:         VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
 50:         VecDotEnd(lmvm->Y[j], ldfp->Q[i], &yjtqi);
 51:         VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
 52:         /* Compute the pure DFP component of the inverse application*/
 53:         VecAXPBYPCZ(ldfp->Q[i], -PetscRealPart(yjtqi) / ldfp->ytq[j], PetscRealPart(sjtyi) / ldfp->yts[j], 1.0, ldfp->Q[j], lmvm->S[j]);
 54:       }
 55:       VecDot(lmvm->Y[i], ldfp->Q[i], &ytq);
 56:       ldfp->ytq[i] = PetscRealPart(ytq);
 57:     }
 58:     ldfp->needQ = PETSC_FALSE;
 59:   }

 61:   /* Start the outer loop (i) for the recursive formula */
 62:   MatSymBrdnApplyJ0Inv(B, F, dX);
 63:   for (i = 0; i <= lmvm->k; ++i) {
 64:     /* Get all the dot products we need */
 65:     VecDotBegin(lmvm->Y[i], dX, &ytx);
 66:     VecDotBegin(lmvm->S[i], F, &stf);
 67:     VecDotEnd(lmvm->Y[i], dX, &ytx);
 68:     VecDotEnd(lmvm->S[i], F, &stf);
 69:     /* Update dX_{i+1} = (B^{-1})_{i+1} * f */
 70:     VecAXPBYPCZ(dX, -PetscRealPart(ytx) / ldfp->ytq[i], PetscRealPart(stf) / ldfp->yts[i], 1.0, ldfp->Q[i], lmvm->S[i]);
 71:   }
 72:   return 0;
 73: }

 75: /*------------------------------------------------------------*/

 77: /*
 78:   The forward product for the approximate Jacobian is the matrix-free
 79:   implementation of the recursive formula given in Equation 6.13 of
 80:   Nocedal and Wright "Numerical Optimization" 2nd edition, pg 139.

 82:   This forward product has a two-loop form similar to the BFGS two-loop
 83:   formulation for the inverse Jacobian application. However, the S and
 84:   Y vectors have interchanged roles.

 86:   work <- X

 88:   for i = k,k-1,k-2,...,0
 89:     rho[i] = 1 / (Y[i]^T S[i])
 90:     alpha[i] = rho[i] * (Y[i]^T work)
 91:     work <- work - (alpha[i] * S[i])
 92:   end

 94:   Z <- J0 * work

 96:   for i = 0,1,2,...,k
 97:     beta = rho[i] * (S[i]^T Y)
 98:     Z <- Z + ((alpha[i] - beta) * Y[i])
 99:   end
100: */
101: PetscErrorCode MatMult_LMVMDFP(Mat B, Vec X, Vec Z)
102: {
103:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
104:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
105:   PetscInt     i;
106:   PetscReal   *alpha, beta;
107:   PetscScalar  ytx, stz;

109:   /* Copy the function into the work vector for the first loop */
110:   VecCopy(X, ldfp->work);

112:   /* Start the first loop */
113:   PetscMalloc1(lmvm->k + 1, &alpha);
114:   for (i = lmvm->k; i >= 0; --i) {
115:     VecDot(lmvm->Y[i], ldfp->work, &ytx);
116:     alpha[i] = PetscRealPart(ytx) / ldfp->yts[i];
117:     VecAXPY(ldfp->work, -alpha[i], lmvm->S[i]);
118:   }

120:   /* Apply the forward product with initial Jacobian */
121:   MatSymBrdnApplyJ0Fwd(B, ldfp->work, Z);

123:   /* Start the second loop */
124:   for (i = 0; i <= lmvm->k; ++i) {
125:     VecDot(lmvm->S[i], Z, &stz);
126:     beta = PetscRealPart(stz) / ldfp->yts[i];
127:     VecAXPY(Z, alpha[i] - beta, lmvm->Y[i]);
128:   }
129:   PetscFree(alpha);
130:   return 0;
131: }

133: /*------------------------------------------------------------*/

135: static PetscErrorCode MatUpdate_LMVMDFP(Mat B, Vec X, Vec F)
136: {
137:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
138:   Mat_SymBrdn  *ldfp = (Mat_SymBrdn *)lmvm->ctx;
139:   Mat_LMVM     *dbase;
140:   Mat_DiagBrdn *dctx;
141:   PetscInt      old_k, i;
142:   PetscReal     curvtol;
143:   PetscScalar   curvature, ytytmp, ststmp;

145:   if (!lmvm->m) return 0;
146:   if (lmvm->prev_set) {
147:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
148:     VecAYPX(lmvm->Xprev, -1.0, X);
149:     VecAYPX(lmvm->Fprev, -1.0, F);
150:     /* Test if the updates can be accepted */
151:     VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
152:     VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
153:     VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
154:     VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
155:     if (PetscRealPart(ststmp) < lmvm->eps) {
156:       curvtol = 0.0;
157:     } else {
158:       curvtol = lmvm->eps * PetscRealPart(ststmp);
159:     }
160:     if (PetscRealPart(curvature) > curvtol) {
161:       /* Update is good, accept it */
162:       ldfp->watchdog = 0;
163:       ldfp->needQ    = PETSC_TRUE;
164:       old_k          = lmvm->k;
165:       MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
166:       /* If we hit the memory limit, shift the yts, yty and sts arrays */
167:       if (old_k == lmvm->k) {
168:         for (i = 0; i <= lmvm->k - 1; ++i) {
169:           ldfp->yts[i] = ldfp->yts[i + 1];
170:           ldfp->yty[i] = ldfp->yty[i + 1];
171:           ldfp->sts[i] = ldfp->sts[i + 1];
172:         }
173:       }
174:       /* Update history of useful scalars */
175:       VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);
176:       ldfp->yts[lmvm->k] = PetscRealPart(curvature);
177:       ldfp->yty[lmvm->k] = PetscRealPart(ytytmp);
178:       ldfp->sts[lmvm->k] = PetscRealPart(ststmp);
179:       /* Compute the scalar scale if necessary */
180:       if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) MatSymBrdnComputeJ0Scalar(B);
181:     } else {
182:       /* Update is bad, skip it */
183:       ++lmvm->nrejects;
184:       ++ldfp->watchdog;
185:     }
186:   } else {
187:     switch (ldfp->scale_type) {
188:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
189:       dbase = (Mat_LMVM *)ldfp->D->data;
190:       dctx  = (Mat_DiagBrdn *)dbase->ctx;
191:       VecSet(dctx->invD, ldfp->delta);
192:       break;
193:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
194:       ldfp->sigma = ldfp->delta;
195:       break;
196:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
197:       ldfp->sigma = 1.0;
198:       break;
199:     default:
200:       break;
201:     }
202:   }

204:   /* Update the scaling */
205:   if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) MatLMVMUpdate(ldfp->D, X, F);

207:   if (ldfp->watchdog > ldfp->max_seq_rejects) {
208:     MatLMVMReset(B, PETSC_FALSE);
209:     if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) MatLMVMReset(ldfp->D, PETSC_FALSE);
210:   }

212:   /* Save the solution and function to be used in the next update */
213:   VecCopy(X, lmvm->Xprev);
214:   VecCopy(F, lmvm->Fprev);
215:   lmvm->prev_set = PETSC_TRUE;
216:   return 0;
217: }

219: /*------------------------------------------------------------*/

221: static PetscErrorCode MatCopy_LMVMDFP(Mat B, Mat M, MatStructure str)
222: {
223:   Mat_LMVM    *bdata = (Mat_LMVM *)B->data;
224:   Mat_SymBrdn *bctx  = (Mat_SymBrdn *)bdata->ctx;
225:   Mat_LMVM    *mdata = (Mat_LMVM *)M->data;
226:   Mat_SymBrdn *mctx  = (Mat_SymBrdn *)mdata->ctx;
227:   PetscInt     i;

229:   mctx->needQ = bctx->needQ;
230:   for (i = 0; i <= bdata->k; ++i) {
231:     mctx->ytq[i] = bctx->ytq[i];
232:     mctx->yts[i] = bctx->yts[i];
233:     VecCopy(bctx->Q[i], mctx->Q[i]);
234:   }
235:   mctx->scale_type      = bctx->scale_type;
236:   mctx->alpha           = bctx->alpha;
237:   mctx->beta            = bctx->beta;
238:   mctx->rho             = bctx->rho;
239:   mctx->sigma_hist      = bctx->sigma_hist;
240:   mctx->watchdog        = bctx->watchdog;
241:   mctx->max_seq_rejects = bctx->max_seq_rejects;
242:   switch (bctx->scale_type) {
243:   case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
244:     mctx->sigma = bctx->sigma;
245:     break;
246:   case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
247:     MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN);
248:     break;
249:   case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
250:     mctx->sigma = 1.0;
251:     break;
252:   default:
253:     break;
254:   }
255:   return 0;
256: }

258: /*------------------------------------------------------------*/

260: static PetscErrorCode MatReset_LMVMDFP(Mat B, PetscBool destructive)
261: {
262:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
263:   Mat_SymBrdn  *ldfp = (Mat_SymBrdn *)lmvm->ctx;
264:   Mat_LMVM     *dbase;
265:   Mat_DiagBrdn *dctx;

267:   ldfp->watchdog = 0;
268:   ldfp->needQ    = PETSC_TRUE;
269:   if (ldfp->allocated) {
270:     if (destructive) {
271:       VecDestroy(&ldfp->work);
272:       PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
273:       VecDestroyVecs(lmvm->m, &ldfp->Q);
274:       switch (ldfp->scale_type) {
275:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
276:         MatLMVMReset(ldfp->D, PETSC_TRUE);
277:         break;
278:       default:
279:         break;
280:       }
281:       ldfp->allocated = PETSC_FALSE;
282:     } else {
283:       switch (ldfp->scale_type) {
284:       case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
285:         ldfp->sigma = ldfp->delta;
286:         break;
287:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
288:         MatLMVMReset(ldfp->D, PETSC_FALSE);
289:         dbase = (Mat_LMVM *)ldfp->D->data;
290:         dctx  = (Mat_DiagBrdn *)dbase->ctx;
291:         VecSet(dctx->invD, ldfp->delta);
292:         break;
293:       case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
294:         ldfp->sigma = 1.0;
295:         break;
296:       default:
297:         break;
298:       }
299:     }
300:   }
301:   MatReset_LMVM(B, destructive);
302:   return 0;
303: }

305: /*------------------------------------------------------------*/

307: static PetscErrorCode MatAllocate_LMVMDFP(Mat B, Vec X, Vec F)
308: {
309:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
310:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;

312:   MatAllocate_LMVM(B, X, F);
313:   if (!ldfp->allocated) {
314:     VecDuplicate(X, &ldfp->work);
315:     PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
316:     if (lmvm->m > 0) VecDuplicateVecs(X, lmvm->m, &ldfp->Q);
317:     switch (ldfp->scale_type) {
318:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
319:       MatLMVMAllocate(ldfp->D, X, F);
320:       break;
321:     default:
322:       break;
323:     }
324:     ldfp->allocated = PETSC_TRUE;
325:   }
326:   return 0;
327: }

329: /*------------------------------------------------------------*/

331: static PetscErrorCode MatDestroy_LMVMDFP(Mat B)
332: {
333:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
334:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;

336:   if (ldfp->allocated) {
337:     VecDestroy(&ldfp->work);
338:     PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
339:     VecDestroyVecs(lmvm->m, &ldfp->Q);
340:     ldfp->allocated = PETSC_FALSE;
341:   }
342:   MatDestroy(&ldfp->D);
343:   PetscFree(lmvm->ctx);
344:   MatDestroy_LMVM(B);
345:   return 0;
346: }

348: /*------------------------------------------------------------*/

350: static PetscErrorCode MatSetUp_LMVMDFP(Mat B)
351: {
352:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
353:   Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
354:   PetscInt     n, N;

356:   MatSetUp_LMVM(B);
357:   if (!ldfp->allocated) {
358:     VecDuplicate(lmvm->Xprev, &ldfp->work);
359:     PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
360:     if (lmvm->m > 0) VecDuplicateVecs(lmvm->Xprev, lmvm->m, &ldfp->Q);
361:     switch (ldfp->scale_type) {
362:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
363:       MatGetLocalSize(B, &n, &n);
364:       MatGetSize(B, &N, &N);
365:       MatSetSizes(ldfp->D, n, n, N, N);
366:       MatSetUp(ldfp->D);
367:       break;
368:     default:
369:       break;
370:     }
371:     ldfp->allocated = PETSC_TRUE;
372:   }
373:   return 0;
374: }

376: /*------------------------------------------------------------*/

378: static PetscErrorCode MatSetFromOptions_LMVMDFP(Mat B, PetscOptionItems *PetscOptionsObject)
379: {
380:   MatSetFromOptions_LMVM(B, PetscOptionsObject);
381:   PetscOptionsHeadBegin(PetscOptionsObject, "DFP method for approximating SPD Jacobian actions (MATLMVMDFP)");
382:   MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject);
383:   PetscOptionsHeadEnd();
384:   return 0;
385: }

387: /*------------------------------------------------------------*/

389: PetscErrorCode MatCreate_LMVMDFP(Mat B)
390: {
391:   Mat_LMVM    *lmvm;
392:   Mat_SymBrdn *ldfp;

394:   MatCreate_LMVMSymBrdn(B);
395:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP);
396:   B->ops->setup          = MatSetUp_LMVMDFP;
397:   B->ops->destroy        = MatDestroy_LMVMDFP;
398:   B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
399:   B->ops->solve          = MatSolve_LMVMDFP;

401:   lmvm                = (Mat_LMVM *)B->data;
402:   lmvm->ops->allocate = MatAllocate_LMVMDFP;
403:   lmvm->ops->reset    = MatReset_LMVMDFP;
404:   lmvm->ops->update   = MatUpdate_LMVMDFP;
405:   lmvm->ops->mult     = MatMult_LMVMDFP;
406:   lmvm->ops->copy     = MatCopy_LMVMDFP;

408:   ldfp        = (Mat_SymBrdn *)lmvm->ctx;
409:   ldfp->needP = PETSC_FALSE;
410:   ldfp->phi   = 1.0;
411:   return 0;
412: }

414: /*------------------------------------------------------------*/

416: /*@
417:    MatCreateLMVMDFP - Creates a limited-memory Davidon-Fletcher-Powell (DFP) matrix
418:    used for approximating Jacobians. L-DFP is symmetric positive-definite by
419:    construction, and is the dual of L-BFGS where Y and S vectors swap roles.

421:    The provided local and global sizes must match the solution and function vectors
422:    used with `MatLMVMUpdate()` and `MatSolve()`. The resulting L-DFP matrix will have
423:    storage vectors allocated with `VecCreateSeq()` in serial and `VecCreateMPI()` in
424:    parallel. To use the L-DFP matrix with other vector types, the matrix must be
425:    created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
426:    This ensures that the internal storage and work vectors are duplicated from the
427:    correct type of vector.

429:    Collective

431:    Input Parameters:
432: +  comm - MPI communicator
433: .  n - number of local rows for storage vectors
434: -  N - global size of the storage vectors

436:    Output Parameter:
437: .  B - the matrix

439:    Options Database Keys:
440: +   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
441: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
442: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
443: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
444: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
445: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

447:    Level: intermediate

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

453: .seealso: [](chapter_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDFP`, `MatCreateLMVMBFGS()`, `MatCreateLMVMSR1()`,
454:           `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`, `MatCreateLMVMSymBrdn()`
455: @*/
456: PetscErrorCode MatCreateLMVMDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
457: {
458:   MatCreate(comm, B);
459:   MatSetSizes(*B, n, n, N, N);
460:   MatSetType(*B, MATLMVMDFP);
461:   MatSetUp(*B);
462:   return 0;
463: }