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:   PetscErrorCode    ierr;
35:   PetscInt          i, j;
36:   PetscScalar       yjtqi, sjtyi, ytx, stf, ytq;

41:   VecCheckSameSize(F, 2, dX, 3);
42:   VecCheckMatCompatible(B, dX, 3, F, 2);

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

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

77: /*------------------------------------------------------------*/

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

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

88:   work <- X

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

96:   Z <- J0 * work

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

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

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

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

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

137: /*------------------------------------------------------------*/

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

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

212:   /* Update the scaling */
213:   if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
214:     MatLMVMUpdate(ldfp->D, X, F);
215:   }

217:   if (ldfp->watchdog > ldfp->max_seq_rejects) {
218:     MatLMVMReset(B, PETSC_FALSE);
219:     if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
220:       MatLMVMReset(ldfp->D, PETSC_FALSE);
221:     }
222:   }

224:   /* Save the solution and function to be used in the next update */
225:   VecCopy(X, lmvm->Xprev);
226:   VecCopy(F, lmvm->Fprev);
227:   lmvm->prev_set = PETSC_TRUE;
228:   return(0);
229: }

231: /*------------------------------------------------------------*/

233: static PetscErrorCode MatCopy_LMVMDFP(Mat B, Mat M, MatStructure str)
234: {
235:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
236:   Mat_SymBrdn       *bctx = (Mat_SymBrdn*)bdata->ctx;
237:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
238:   Mat_SymBrdn       *mctx = (Mat_SymBrdn*)mdata->ctx;
239:   PetscErrorCode    ierr;
240:   PetscInt          i;

243:   mctx->needQ = bctx->needQ;
244:   for (i=0; i<=bdata->k; ++i) {
245:     mctx->ytq[i] = bctx->ytq[i];
246:     mctx->yts[i] = bctx->yts[i];
247:     VecCopy(bctx->Q[i], mctx->Q[i]);
248:   }
249:   mctx->scale_type      = bctx->scale_type;
250:   mctx->alpha           = bctx->alpha;
251:   mctx->beta            = bctx->beta;
252:   mctx->rho             = bctx->rho;
253:   mctx->sigma_hist      = bctx->sigma_hist;
254:   mctx->watchdog        = bctx->watchdog;
255:   mctx->max_seq_rejects = bctx->max_seq_rejects;
256:   switch (bctx->scale_type) {
257:   case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
258:     mctx->sigma = bctx->sigma;
259:     break;
260:   case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
261:     MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN);
262:     break;
263:   case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
264:     mctx->sigma = 1.0;
265:     break;
266:   default:
267:     break;
268:   }
269:   return(0);
270: }

272: /*------------------------------------------------------------*/

274: static PetscErrorCode MatReset_LMVMDFP(Mat B, PetscBool destructive)
275: {
276:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
277:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
278:   Mat_LMVM          *dbase;
279:   Mat_DiagBrdn      *dctx;
280:   PetscErrorCode    ierr;

283:   ldfp->watchdog = 0;
284:   ldfp->needQ = PETSC_TRUE;
285:   if (ldfp->allocated) {
286:     if (destructive) {
287:       VecDestroy(&ldfp->work);
288:       PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
289:       VecDestroyVecs(lmvm->m, &ldfp->Q);
290:       switch (ldfp->scale_type) {
291:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
292:         MatLMVMReset(ldfp->D, PETSC_TRUE);
293:         break;
294:       default:
295:         break;
296:       }
297:       ldfp->allocated = PETSC_FALSE;
298:     } else {
299:       switch (ldfp->scale_type) {
300:       case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
301:         ldfp->sigma = ldfp->delta;
302:         break;
303:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
304:         MatLMVMReset(ldfp->D, PETSC_FALSE);
305:         dbase = (Mat_LMVM*)ldfp->D->data;
306:         dctx = (Mat_DiagBrdn*)dbase->ctx;
307:         VecSet(dctx->invD, ldfp->delta);
308:         break;
309:       case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
310:         ldfp->sigma = 1.0;
311:         break;
312:       default:
313:         break;
314:       }
315:     }
316:   }
317:   MatReset_LMVM(B, destructive);
318:   return(0);
319: }

321: /*------------------------------------------------------------*/

323: static PetscErrorCode MatAllocate_LMVMDFP(Mat B, Vec X, Vec F)
324: {
325:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
326:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
327:   PetscErrorCode    ierr;

330:   MatAllocate_LMVM(B, X, F);
331:   if (!ldfp->allocated) {
332:     VecDuplicate(X, &ldfp->work);
333:     PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
334:     if (lmvm->m > 0) {
335:       VecDuplicateVecs(X, lmvm->m, &ldfp->Q);
336:     }
337:     switch (ldfp->scale_type) {
338:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
339:       MatLMVMAllocate(ldfp->D, X, F);
340:       break;
341:     default:
342:       break;
343:     }
344:     ldfp->allocated = PETSC_TRUE;
345:   }
346:   return(0);
347: }

349: /*------------------------------------------------------------*/

351: static PetscErrorCode MatDestroy_LMVMDFP(Mat B)
352: {
353:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
354:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
355:   PetscErrorCode    ierr;

358:   if (ldfp->allocated) {
359:     VecDestroy(&ldfp->work);
360:     PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
361:     VecDestroyVecs(lmvm->m, &ldfp->Q);
362:     ldfp->allocated = PETSC_FALSE;
363:   }
364:   MatDestroy(&ldfp->D);
365:   PetscFree(lmvm->ctx);
366:   MatDestroy_LMVM(B);
367:   return(0);
368: }

370: /*------------------------------------------------------------*/

372: static PetscErrorCode MatSetUp_LMVMDFP(Mat B)
373: {
374:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
375:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
376:   PetscErrorCode    ierr;
377:   PetscInt          n, N;

380:   MatSetUp_LMVM(B);
381:   if (!ldfp->allocated) {
382:     VecDuplicate(lmvm->Xprev, &ldfp->work);
383:     PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
384:     if (lmvm->m > 0) {
385:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &ldfp->Q);
386:     }
387:     switch (ldfp->scale_type) {
388:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
389:       MatGetLocalSize(B, &n, &n);
390:       MatGetSize(B, &N, &N);
391:       MatSetSizes(ldfp->D, n, n, N, N);
392:       MatSetUp(ldfp->D);
393:       break;
394:     default:
395:       break;
396:     }
397:     ldfp->allocated = PETSC_TRUE;
398:   }
399:   return(0);
400: }

402: /*------------------------------------------------------------*/

404: static PetscErrorCode MatSetFromOptions_LMVMDFP(PetscOptionItems *PetscOptionsObject, Mat B)
405: {
406:   PetscErrorCode               ierr;

409:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
410:   PetscOptionsHead(PetscOptionsObject,"DFP method for approximating SPD Jacobian actions (MATLMVMDFP)");
411:   MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
412:   PetscOptionsTail();
413:   return(0);
414: }

416: /*------------------------------------------------------------*/

418: PetscErrorCode MatCreate_LMVMDFP(Mat B)
419: {
420:   Mat_LMVM          *lmvm;
421:   Mat_SymBrdn       *ldfp;
422:   PetscErrorCode    ierr;

425:   MatCreate_LMVMSymBrdn(B);
426:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP);
427:   B->ops->setup = MatSetUp_LMVMDFP;
428:   B->ops->destroy = MatDestroy_LMVMDFP;
429:   B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
430:   B->ops->solve = MatSolve_LMVMDFP;

432:   lmvm = (Mat_LMVM*)B->data;
433:   lmvm->ops->allocate = MatAllocate_LMVMDFP;
434:   lmvm->ops->reset = MatReset_LMVMDFP;
435:   lmvm->ops->update = MatUpdate_LMVMDFP;
436:   lmvm->ops->mult = MatMult_LMVMDFP;
437:   lmvm->ops->copy = MatCopy_LMVMDFP;

439:   ldfp = (Mat_SymBrdn*)lmvm->ctx;
440:   ldfp->needP           = PETSC_FALSE;
441:   ldfp->phi             = 1.0;
442:   return(0);
443: }

445: /*------------------------------------------------------------*/

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

452:    The provided local and global sizes must match the solution and function vectors
453:    used with MatLMVMUpdate() and MatSolve(). The resulting L-DFP matrix will have
454:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
455:    parallel. To use the L-DFP matrix with other vector types, the matrix must be
456:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
457:    This ensures that the internal storage and work vectors are duplicated from the
458:    correct type of vector.

460:    Collective

462:    Input Parameters:
463: +  comm - MPI communicator, set to PETSC_COMM_SELF
464: .  n - number of local rows for storage vectors
465: -  N - global size of the storage vectors

467:    Output Parameter:
468: .  B - the matrix

470:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()

473:    Options Database Keys:
474: +   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
475: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
476: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
477: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
478: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
479: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
480: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

482:    Level: intermediate

484: .seealso: MatCreate(), MATLMVM, MATLMVMDFP, MatCreateLMVMBFGS(), MatCreateLMVMSR1(),