Actual source code: bfgs.c

petsc-master 2019-06-15
Report Typos and Errors
  1:  #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
  2:  #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>

  4: /*
  5:   Limited-memory Broyden-Fletcher-Goldfarb-Shano 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 adapted 
 13:    from Algorithm 7.4 on page 178 of Nocedal and Wright "Numerical Optimization" 
 14:    2nd edition (https://doi.org/10.1007/978-0-387-40065-5). The initial inverse 
 15:    Jacobian application falls back onto the gamma scaling recommended in equation 
 16:    (7.20) if the user has not provided any estimation of the initial Jacobian or 
 17:    its inverse.

 19:    work <- F

 21:    for i = k,k-1,k-2,...,0
 22:      rho[i] = 1 / (Y[i]^T S[i])
 23:      alpha[i] = rho[i] * (S[i]^T work)
 24:      Fwork <- work - (alpha[i] * Y[i])
 25:    end

 27:    dX <- J0^{-1} * work

 29:    for i = 0,1,2,...,k
 30:      beta = rho[i] * (Y[i]^T dX)
 31:      dX <- dX + ((alpha[i] - beta) * S[i])
 32:    end
 33: */
 34: PetscErrorCode MatSolve_LMVMBFGS(Mat B, Vec F, Vec dX)
 35: {
 36:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
 37:   Mat_SymBrdn       *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
 38:   PetscErrorCode    ierr;
 39:   PetscInt          i;
 40:   PetscReal         *alpha, beta;
 41:   PetscScalar       stf, ytx;
 42: 
 44:   VecCheckSameSize(F, 2, dX, 3);
 45:   VecCheckMatCompatible(B, dX, 3, F, 2);
 46: 
 47:   /* Copy the function into the work vector for the first loop */
 48:   VecCopy(F, lbfgs->work);
 49: 
 50:   /* Start the first loop */
 51:   PetscMalloc1(lmvm->k+1, &alpha);
 52:   for (i = lmvm->k; i >= 0; --i) {
 53:     VecDot(lmvm->S[i], lbfgs->work, &stf);
 54:     alpha[i] = PetscRealPart(stf)/lbfgs->yts[i];
 55:     VecAXPY(lbfgs->work, -alpha[i], lmvm->Y[i]);
 56:   }
 57: 
 58:   /* Invert the initial Jacobian onto the work vector (or apply scaling) */
 59:   MatSymBrdnApplyJ0Inv(B, lbfgs->work, dX);
 60: 
 61:   /* Start the second loop */
 62:   for (i = 0; i <= lmvm->k; ++i) {
 63:     VecDot(lmvm->Y[i], dX, &ytx);
 64:     beta = PetscRealPart(ytx)/lbfgs->yts[i];
 65:     VecAXPY(dX, alpha[i]-beta, lmvm->S[i]);
 66:   }
 67:   PetscFree(alpha);
 68:   return(0);
 69: }

 71: /*------------------------------------------------------------*/

 73: /*
 74:   The forward product for the approximate Jacobian is the matrix-free 
 75:   implementation of Equation (6.19) in Nocedal and Wright "Numerical 
 76:   Optimization" 2nd Edition, pg 140.
 77:   
 78:   This forward product has the same structure as the inverse Jacobian 
 79:   application in the DFP formulation, except with S and Y exchanging 
 80:   roles.
 81:   
 82:   Note: P[i] = (B_i)*S[i] terms are computed ahead of time whenever 
 83:   the matrix is updated with a new (S[i], Y[i]) pair. This allows 
 84:   repeated calls of MatMult inside KSP solvers without unnecessarily 
 85:   recomputing P[i] terms in expensive nested-loops.

 87:   Z <- J0 * X

 89:   for i = 0,1,2,...,k
 90:     P[i] <- J0 * S[i]
 91:     for j = 0,1,2,...,(i-1)
 92:       gamma = (Y[j]^T S[i]) / (Y[j]^T S[j])
 93:       zeta = (S[j]^ P[i]) / (S[j]^T P[j])
 94:       P[i] <- P[i] - (zeta * P[j]) + (gamma * Y[j])
 95:     end
 96:     gamma = (Y[i]^T X) / (Y[i]^T S[i])
 97:     zeta = (S[i]^T Z) / (S[i]^T P[i])
 98:     Z <- Z - (zeta * P[i]) + (gamma * Y[i])
 99:   end
100: */
101: PetscErrorCode MatMult_LMVMBFGS(Mat B, Vec X, Vec Z)
102: {
103:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
104:   Mat_SymBrdn       *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
105:   PetscErrorCode    ierr;
106:   PetscInt          i, j;
107:   PetscScalar       sjtpi, yjtsi, ytx, stz, stp;
108: 
110:   VecCheckSameSize(X, 2, Z, 3);
111:   VecCheckMatCompatible(B, X, 2, Z, 3);
112: 
113:   if (lbfgs->needP) {
114:     /* Pre-compute (P[i] = B_i * S[i]) */
115:     for (i = 0; i <= lmvm->k; ++i) {
116:       MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lbfgs->P[i]);
117:       for (j = 0; j <= i-1; ++j) {
118:         /* Compute the necessary dot products */
119:         VecDotBegin(lmvm->S[j], lbfgs->P[i], &sjtpi);
120:         VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
121:         VecDotEnd(lmvm->S[j], lbfgs->P[i], &sjtpi);
122:         VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
123:         /* Compute the pure BFGS component of the forward product */
124:         VecAXPBYPCZ(lbfgs->P[i], -PetscRealPart(sjtpi)/lbfgs->stp[j], PetscRealPart(yjtsi)/lbfgs->yts[j], 1.0, lbfgs->P[j], lmvm->Y[j]);
125:       }
126:       VecDot(lmvm->S[i], lbfgs->P[i], &stp);
127:       lbfgs->stp[i] = PetscRealPart(stp);
128:     }
129:     lbfgs->needP = PETSC_FALSE;
130:   }
131: 
132:   /* Start the outer loop (i) for the recursive formula */
133:   MatSymBrdnApplyJ0Fwd(B, X, Z);
134:   for (i = 0; i <= lmvm->k; ++i) {
135:     /* Get all the dot products we need */
136:     VecDotBegin(lmvm->S[i], Z, &stz);
137:     VecDotBegin(lmvm->Y[i], X, &ytx);
138:     VecDotEnd(lmvm->S[i], Z, &stz);
139:     VecDotEnd(lmvm->Y[i], X, &ytx);
140:     /* Update Z_{i+1} = B_{i+1} * X */
141:     VecAXPBYPCZ(Z, -PetscRealPart(stz)/lbfgs->stp[i], PetscRealPart(ytx)/lbfgs->yts[i], 1.0, lbfgs->P[i], lmvm->Y[i]);
142:   }
143:   return(0);
144: }

146: /*------------------------------------------------------------*/

148: static PetscErrorCode MatUpdate_LMVMBFGS(Mat B, Vec X, Vec F)
149: {
150:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
151:   Mat_SymBrdn       *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
152:   Mat_LMVM          *dbase;
153:   Mat_DiagBrdn      *dctx;
154:   PetscErrorCode    ierr;
155:   PetscInt          old_k, i;
156:   PetscReal         curvtol;
157:   PetscScalar       curvature, ytytmp, ststmp;

160:   if (!lmvm->m) return(0);
161:   if (lmvm->prev_set) {
162:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
163:     VecAYPX(lmvm->Xprev, -1.0, X);
164:     VecAYPX(lmvm->Fprev, -1.0, F);
165:     /* Test if the updates can be accepted */
166:     VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
167:     VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
168:     VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
169:     VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
170:     if (PetscRealPart(ststmp) < lmvm->eps) {
171:       curvtol = 0.0;
172:     } else {
173:       curvtol = lmvm->eps * PetscRealPart(ststmp);
174:     }
175:     if (PetscRealPart(curvature) > curvtol) {
176:       /* Update is good, accept it */
177:       lbfgs->watchdog = 0;
178:       lbfgs->needP = PETSC_TRUE;
179:       old_k = lmvm->k;
180:       MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
181:       /* If we hit the memory limit, shift the yts, yty and sts arrays */
182:       if (old_k == lmvm->k) {
183:         for (i = 0; i <= lmvm->k-1; ++i) {
184:           lbfgs->yts[i] = lbfgs->yts[i+1];
185:           lbfgs->yty[i] = lbfgs->yty[i+1];
186:           lbfgs->sts[i] = lbfgs->sts[i+1];
187:         }
188:       }
189:       /* Update history of useful scalars */
190:       VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);
191:       lbfgs->yts[lmvm->k] = PetscRealPart(curvature);
192:       lbfgs->yty[lmvm->k] = PetscRealPart(ytytmp);
193:       lbfgs->sts[lmvm->k] = PetscRealPart(ststmp);
194:       /* Compute the scalar scale if necessary */
195:       if (lbfgs->scale_type == SYMBRDN_SCALE_SCALAR) {
196:         MatSymBrdnComputeJ0Scalar(B);
197:       }
198:     } else {
199:       /* Update is bad, skip it */
200:       ++lmvm->nrejects;
201:       ++lbfgs->watchdog;
202:     }
203:   } else {
204:     switch (lbfgs->scale_type) {
205:     case SYMBRDN_SCALE_DIAG:
206:       dbase = (Mat_LMVM*)lbfgs->D->data;
207:       dctx = (Mat_DiagBrdn*)dbase->ctx;
208:       VecSet(dctx->invD, lbfgs->delta);
209:       break;
210:     case SYMBRDN_SCALE_SCALAR:
211:       lbfgs->sigma = lbfgs->delta;
212:       break;
213:     case SYMBRDN_SCALE_NONE:
214:       lbfgs->sigma = 1.0;
215:       break;
216:     default:
217:       break;
218:     }
219:   }
220: 
221:   /* Update the scaling */
222:   if (lbfgs->scale_type == SYMBRDN_SCALE_DIAG) {
223:     MatLMVMUpdate(lbfgs->D, X, F);
224:   }
225: 
226:   if (lbfgs->watchdog > lbfgs->max_seq_rejects) {
227:     MatLMVMReset(B, PETSC_FALSE);
228:     if (lbfgs->scale_type == SYMBRDN_SCALE_DIAG) {
229:       MatLMVMReset(lbfgs->D, PETSC_FALSE);
230:     }
231:   }

233:   /* Save the solution and function to be used in the next update */
234:   VecCopy(X, lmvm->Xprev);
235:   VecCopy(F, lmvm->Fprev);
236:   lmvm->prev_set = PETSC_TRUE;
237:   return(0);
238: }

240: /*------------------------------------------------------------*/

242: static PetscErrorCode MatCopy_LMVMBFGS(Mat B, Mat M, MatStructure str)
243: {
244:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
245:   Mat_SymBrdn       *bctx = (Mat_SymBrdn*)bdata->ctx;
246:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
247:   Mat_SymBrdn       *mctx = (Mat_SymBrdn*)mdata->ctx;
248:   PetscErrorCode    ierr;
249:   PetscInt          i;

252:   mctx->needP = bctx->needP;
253:   for (i=0; i<=bdata->k; ++i) {
254:     mctx->stp[i] = bctx->stp[i];
255:     mctx->yts[i] = bctx->yts[i];
256:     VecCopy(bctx->P[i], mctx->P[i]);
257:   }
258:   mctx->scale_type      = bctx->scale_type;
259:   mctx->alpha           = bctx->alpha;
260:   mctx->beta            = bctx->beta;
261:   mctx->rho             = bctx->rho;
262:   mctx->delta           = bctx->delta;
263:   mctx->sigma_hist      = bctx->sigma_hist;
264:   mctx->watchdog        = bctx->watchdog;
265:   mctx->max_seq_rejects = bctx->max_seq_rejects;
266:   switch (bctx->scale_type) {
267:   case SYMBRDN_SCALE_SCALAR:
268:     mctx->sigma = bctx->sigma;
269:     break;
270:   case SYMBRDN_SCALE_DIAG:
271:     MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN);
272:     break;
273:   case SYMBRDN_SCALE_NONE:
274:     mctx->sigma = 1.0;
275:     break;
276:   default:
277:     break;
278:   }
279:   return(0);
280: }

282: /*------------------------------------------------------------*/

284: static PetscErrorCode MatReset_LMVMBFGS(Mat B, PetscBool destructive)
285: {
286:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
287:   Mat_SymBrdn       *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
288:   Mat_LMVM          *dbase;
289:   Mat_DiagBrdn      *dctx;
290:   PetscErrorCode    ierr;
291: 
293:   lbfgs->watchdog = 0;
294:   lbfgs->needP = PETSC_TRUE;
295:   if (lbfgs->allocated) {
296:     if (destructive) {
297:       VecDestroy(&lbfgs->work);
298:       PetscFree4(lbfgs->stp, lbfgs->yts, lbfgs->yty, lbfgs->sts);
299:       VecDestroyVecs(lmvm->m, &lbfgs->P);
300:       switch (lbfgs->scale_type) {
301:       case SYMBRDN_SCALE_DIAG:
302:         MatLMVMReset(lbfgs->D, PETSC_TRUE);
303:         break;
304:       default:
305:         break;
306:       }
307:       lbfgs->allocated = PETSC_FALSE;
308:     } else {
309:       switch (lbfgs->scale_type) {
310:       case SYMBRDN_SCALE_SCALAR:
311:         lbfgs->sigma = lbfgs->delta;
312:         break;
313:       case SYMBRDN_SCALE_DIAG:
314:         MatLMVMReset(lbfgs->D, PETSC_FALSE);
315:         dbase = (Mat_LMVM*)lbfgs->D->data;
316:         dctx = (Mat_DiagBrdn*)dbase->ctx;
317:         VecSet(dctx->invD, lbfgs->delta);
318:         break;
319:       case SYMBRDN_SCALE_NONE:
320:         lbfgs->sigma = 1.0;
321:         break;
322:       default:
323:         break;
324:       }
325:     }
326:   }
327:   MatReset_LMVM(B, destructive);
328:   return(0);
329: }

331: /*------------------------------------------------------------*/

333: static PetscErrorCode MatAllocate_LMVMBFGS(Mat B, Vec X, Vec F)
334: {
335:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
336:   Mat_SymBrdn       *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
337:   PetscErrorCode    ierr;
338: 
340:   MatAllocate_LMVM(B, X, F);
341:   if (!lbfgs->allocated) {
342:     VecDuplicate(X, &lbfgs->work);
343:     PetscMalloc4(lmvm->m, &lbfgs->stp, lmvm->m, &lbfgs->yts, lmvm->m, &lbfgs->yty, lmvm->m, &lbfgs->sts);
344:     if (lmvm->m > 0) {
345:       VecDuplicateVecs(X, lmvm->m, &lbfgs->P);
346:     }
347:     switch (lbfgs->scale_type) {
348:     case SYMBRDN_SCALE_DIAG:
349:       MatLMVMAllocate(lbfgs->D, X, F);
350:       break;
351:     default:
352:       break;
353:     }
354:     lbfgs->allocated = PETSC_TRUE;
355:   }
356:   return(0);
357: }

359: /*------------------------------------------------------------*/

361: static PetscErrorCode MatDestroy_LMVMBFGS(Mat B)
362: {
363:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
364:   Mat_SymBrdn       *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
365:   PetscErrorCode    ierr;

368:   if (lbfgs->allocated) {
369:     VecDestroy(&lbfgs->work);
370:     PetscFree4(lbfgs->stp, lbfgs->yts, lbfgs->yty, lbfgs->sts);
371:     VecDestroyVecs(lmvm->m, &lbfgs->P);
372:     lbfgs->allocated = PETSC_FALSE;
373:   }
374:   MatDestroy(&lbfgs->D);
375:   PetscFree(lmvm->ctx);
376:   MatDestroy_LMVM(B);
377:   return(0);
378: }

380: /*------------------------------------------------------------*/

382: static PetscErrorCode MatSetUp_LMVMBFGS(Mat B)
383: {
384:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
385:   Mat_SymBrdn       *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
386:   PetscErrorCode    ierr;
387:   PetscInt          n, N;
388: 
390:   MatSetUp_LMVM(B);
391:   if (!lbfgs->allocated) {
392:     VecDuplicate(lmvm->Xprev, &lbfgs->work);
393:     PetscMalloc4(lmvm->m, &lbfgs->stp, lmvm->m, &lbfgs->yts, lmvm->m, &lbfgs->yty, lmvm->m, &lbfgs->sts);
394:     if (lmvm->m > 0) {
395:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbfgs->P);
396:     }
397:     switch (lbfgs->scale_type) {
398:     case SYMBRDN_SCALE_DIAG:
399:       MatGetLocalSize(B, &n, &n);
400:       MatGetSize(B, &N, &N);
401:       MatSetSizes(lbfgs->D, n, n, N, N);
402:       MatSetUp(lbfgs->D);
403:       break;
404:     default:
405:       break;
406:     }
407:     lbfgs->allocated = PETSC_TRUE;
408:   }
409:   return(0);
410: }

412: /*------------------------------------------------------------*/

414: static PetscErrorCode MatSetFromOptions_LMVMBFGS(PetscOptionItems *PetscOptionsObject, Mat B)
415: {
416:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
417:   Mat_SymBrdn       *lbfgs = (Mat_SymBrdn*)lmvm->ctx;
418:   Mat_LMVM          *dbase;
419:   Mat_DiagBrdn      *dctx;
420:   PetscErrorCode    ierr;

423:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
424:   PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
425:   PetscOptionsEList("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "", Scale_Table, SYMBRDN_SCALE_SIZE, Scale_Table[lbfgs->scale_type], &lbfgs->scale_type,NULL);
426:   PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",lbfgs->theta,&lbfgs->theta,NULL);
427:   PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",lbfgs->rho,&lbfgs->rho,NULL);
428:   PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",lbfgs->alpha,&lbfgs->alpha,NULL);
429:   PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",lbfgs->beta,&lbfgs->beta,NULL);
430:   PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",lbfgs->sigma_hist,&lbfgs->sigma_hist,NULL);
431:   PetscOptionsTail();
432:   if ((lbfgs->theta < 0.0) || (lbfgs->theta > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
433:   if ((lbfgs->alpha < 0.0) || (lbfgs->alpha > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
434:   if ((lbfgs->rho < 0.0) || (lbfgs->rho > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "update limiter in the J0 scaling cannot be outside the range of [0, 1]");
435:   if (lbfgs->sigma_hist < 0) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
436:   if (lbfgs->scale_type == SYMBRDN_SCALE_DIAG) {
437:     MatSetFromOptions(lbfgs->D);
438:     dbase = (Mat_LMVM*)lbfgs->D->data;
439:     dctx = (Mat_DiagBrdn*)dbase->ctx;
440:     dctx->delta_min  = lbfgs->delta_min;
441:     dctx->delta_max  = lbfgs->delta_max;
442:     dctx->theta      = lbfgs->theta;
443:     dctx->rho        = lbfgs->rho;
444:     dctx->alpha      = lbfgs->alpha;
445:     dctx->beta       = lbfgs->beta;
446:     dctx->sigma_hist = lbfgs->sigma_hist;
447:   }
448:   return(0);
449: }

451: /*------------------------------------------------------------*/

453: PetscErrorCode MatCreate_LMVMBFGS(Mat B)
454: {
455:   Mat_LMVM          *lmvm;
456:   Mat_SymBrdn       *lbfgs;
457:   PetscErrorCode    ierr;

460:   MatCreate_LMVM(B);
461:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMBFGS);
462:   MatSetOption(B, MAT_SPD, PETSC_TRUE);
463:   B->ops->view = MatView_LMVMSymBrdn;
464:   B->ops->setup = MatSetUp_LMVMBFGS;
465:   B->ops->destroy = MatDestroy_LMVMBFGS;
466:   B->ops->setfromoptions = MatSetFromOptions_LMVMBFGS;
467:   B->ops->solve = MatSolve_LMVMBFGS;

469:   lmvm = (Mat_LMVM*)B->data;
470:   lmvm->square = PETSC_TRUE;
471:   lmvm->ops->allocate = MatAllocate_LMVMBFGS;
472:   lmvm->ops->reset = MatReset_LMVMBFGS;
473:   lmvm->ops->update = MatUpdate_LMVMBFGS;
474:   lmvm->ops->mult = MatMult_LMVMBFGS;
475:   lmvm->ops->copy = MatCopy_LMVMBFGS;

477:   PetscNewLog(B, &lbfgs);
478:   lmvm->ctx = (void*)lbfgs;
479:   lbfgs->allocated       = PETSC_FALSE;
480:   lbfgs->needP           = PETSC_TRUE;
481:   lbfgs->phi             = 0.0;
482:   lbfgs->theta           = 0.125;
483:   lbfgs->alpha           = 1.0;
484:   lbfgs->rho             = 1.0;
485:   lbfgs->beta            = 0.5;
486:   lbfgs->sigma           = 1.0;
487:   lbfgs->delta           = 1.0;
488:   lbfgs->delta_min       = 1e-7;
489:   lbfgs->delta_max       = 100.0;
490:   lbfgs->sigma_hist      = 1;
491:   lbfgs->scale_type      = SYMBRDN_SCALE_DIAG;
492:   lbfgs->watchdog        = 0;
493:   lbfgs->max_seq_rejects = lmvm->m/2;
494: 
495:   MatCreate(PetscObjectComm((PetscObject)B), &lbfgs->D);
496:   MatSetType(lbfgs->D, MATLMVMDIAGBRDN);
497:   MatSetOptionsPrefix(lbfgs->D, "J0_");
498:   return(0);
499: }

501: /*------------------------------------------------------------*/

503: /*@
504:    MatCreateLMVMBFGS - Creates a limited-memory Broyden-Fletcher-Goldfarb-Shano (BFGS)
505:    matrix used for approximating Jacobians. L-BFGS is symmetric positive-definite by 
506:    construction, and is commonly used to approximate Hessians in optimization 
507:    problems.
508:    
509:    The provided local and global sizes must match the solution and function vectors 
510:    used with MatLMVMUpdate() and MatSolve(). The resulting L-BFGS matrix will have 
511:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in 
512:    parallel. To use the L-BFGS matrix with other vector types, the matrix must be 
513:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate(). 
514:    This ensures that the internal storage and work vectors are duplicated from the 
515:    correct type of vector.

517:    Collective

519:    Input Parameters:
520: +  comm - MPI communicator, set to PETSC_COMM_SELF
521: .  n - number of local rows for storage vectors
522: -  N - global size of the storage vectors

524:    Output Parameter:
525: .  B - the matrix

527:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
528:    paradigm instead of this routine directly.

530:    Options Database Keys:
531: .   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
532: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
533: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
534: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
535: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
536: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
537: .   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

539:    Level: intermediate

541: .seealso: MatCreate(), MATLMVM, MATLMVMBFGS, MatCreateLMVMDFP(), MatCreateLMVMSR1(), 
542:           MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn(), MatCreateLMVMSymBrdn()
543: @*/
544: PetscErrorCode MatCreateLMVMBFGS(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
545: {
546:   PetscErrorCode    ierr;
547: 
549:   MatCreate(comm, B);
550:   MatSetSizes(*B, n, n, N, N);
551:   MatSetType(*B, MATLMVMBFGS);
552:   MatSetUp(*B);
553:   return(0);
554: }