Actual source code: dfp.c

petsc-3.11.3 2019-06-26
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 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.
 16:   
 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;
 37: 
 41:   VecCheckSameSize(F, 2, dX, 3);
 42:   VecCheckMatCompatible(B, dX, 3, F, 2);
 43: 
 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:   }
 62: 
 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.
 83:   
 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;
111: 
113:   /* Copy the function into the work vector for the first loop */
114:   VecCopy(X, ldfp->work);
115: 
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:   }
123: 
124:   /* Apply the forward product with initial Jacobian */
125:   MatSymBrdnApplyJ0Fwd(B, ldfp->work, Z);
126: 
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 == SYMBRDN_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 SYMBRDN_SCALE_DIAG:
197:       dbase = (Mat_LMVM*)ldfp->D->data;
198:       dctx = (Mat_DiagBrdn*)dbase->ctx;
199:       VecSet(dctx->invD, ldfp->delta);
200:       break;
201:     case SYMBRDN_SCALE_SCALAR:
202:       ldfp->sigma = ldfp->delta;
203:       break;
204:     case SYMBRDN_SCALE_NONE:
205:       ldfp->sigma = 1.0;
206:       break;
207:     default:
208:       break;
209:     }
210:   }
211: 
212:   /* Update the scaling */
213:   if (ldfp->scale_type == SYMBRDN_SCALE_DIAG) {
214:     MatLMVMUpdate(ldfp->D, X, F);
215:   }
216: 
217:   if (ldfp->watchdog > ldfp->max_seq_rejects) {
218:     MatLMVMReset(B, PETSC_FALSE);
219:     if (ldfp->scale_type == SYMBRDN_SCALE_DIAG) {
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 SYMBRDN_SCALE_SCALAR:
258:     mctx->sigma = bctx->sigma;
259:     break;
260:   case SYMBRDN_SCALE_DIAG:
261:     MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN);
262:     break;
263:   case SYMBRDN_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;
281: 
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 SYMBRDN_SCALE_DIAG:
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 SYMBRDN_SCALE_SCALAR:
301:         ldfp->sigma = ldfp->delta;
302:         break;
303:       case SYMBRDN_SCALE_DIAG:
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 SYMBRDN_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;
328: 
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 SYMBRDN_SCALE_DIAG:
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;
378: 
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 SYMBRDN_SCALE_DIAG:
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:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
407:   Mat_SymBrdn       *ldfp = (Mat_SymBrdn*)lmvm->ctx;
408:   Mat_LMVM          *dbase;
409:   Mat_DiagBrdn      *dctx;
410:   PetscErrorCode    ierr;

413:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
414:   PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
415:   PetscOptionsEList("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "", Scale_Table, SYMBRDN_SCALE_SIZE, Scale_Table[ldfp->scale_type], &ldfp->scale_type,NULL);
416:   PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",ldfp->theta,&ldfp->theta,NULL);
417:   PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",ldfp->rho,&ldfp->rho,NULL);
418:   PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",ldfp->alpha,&ldfp->alpha,NULL);
419:   PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",ldfp->beta,&ldfp->beta,NULL);
420:   PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",ldfp->sigma_hist,&ldfp->sigma_hist,NULL);
421:   PetscOptionsTail();
422:   if ((ldfp->theta < 0.0) || (ldfp->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]");
423:   if ((ldfp->alpha < 0.0) || (ldfp->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]");
424:   if ((ldfp->rho < 0.0) || (ldfp->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]");
425:   if (ldfp->sigma_hist < 0) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
426:   if (ldfp->scale_type == SYMBRDN_SCALE_DIAG) {
427:     MatSetFromOptions(ldfp->D);
428:     dbase = (Mat_LMVM*)ldfp->D->data;
429:     dctx = (Mat_DiagBrdn*)dbase->ctx;
430:     dctx->delta_min  = ldfp->delta_min;
431:     dctx->delta_max  = ldfp->delta_max;
432:     dctx->theta      = ldfp->theta;
433:     dctx->rho        = ldfp->rho;
434:     dctx->alpha      = ldfp->alpha;
435:     dctx->beta       = ldfp->beta;
436:     dctx->sigma_hist = ldfp->sigma_hist;
437:   }
438:   return(0);
439: }

441: /*------------------------------------------------------------*/

443: PetscErrorCode MatCreate_LMVMDFP(Mat B)
444: {
445:   Mat_LMVM          *lmvm;
446:   Mat_SymBrdn       *ldfp;
447:   PetscErrorCode    ierr;

450:   MatCreate_LMVM(B);
451:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP);
452:   MatSetOption(B, MAT_SPD, PETSC_TRUE);
453:   B->ops->view = MatView_LMVMSymBrdn;
454:   B->ops->setup = MatSetUp_LMVMDFP;
455:   B->ops->destroy = MatDestroy_LMVMDFP;
456:   B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
457:   B->ops->solve = MatSolve_LMVMDFP;

459:   lmvm = (Mat_LMVM*)B->data;
460:   lmvm->square = PETSC_TRUE;
461:   lmvm->ops->allocate = MatAllocate_LMVMDFP;
462:   lmvm->ops->reset = MatReset_LMVMDFP;
463:   lmvm->ops->update = MatUpdate_LMVMDFP;
464:   lmvm->ops->mult = MatMult_LMVMDFP;
465:   lmvm->ops->copy = MatCopy_LMVMDFP;

467:   PetscNewLog(B, &ldfp);
468:   lmvm->ctx = (void*)ldfp;
469:   ldfp->allocated       = PETSC_FALSE;
470:   ldfp->needQ           = PETSC_TRUE;
471:   ldfp->phi             = 1.0;
472:   ldfp->theta           = 0.125;
473:   ldfp->alpha           = 1.0;
474:   ldfp->rho             = 1.0;
475:   ldfp->beta            = 0.5;
476:   ldfp->sigma           = 1.0;
477:   ldfp->delta           = 1.0;
478:   ldfp->delta_min       = 1e-7;
479:   ldfp->delta_max       = 100.0;
480:   ldfp->sigma_hist      = 1;
481:   ldfp->scale_type      = SYMBRDN_SCALE_DIAG;
482:   ldfp->watchdog        = 0;
483:   ldfp->max_seq_rejects = lmvm->m/2;
484: 
485:   MatCreate(PetscObjectComm((PetscObject)B), &ldfp->D);
486:   MatSetType(ldfp->D, MATLMVMDIAGBRDN);
487:   MatSetOptionsPrefix(ldfp->D, "J0_");
488:   return(0);
489: }

491: /*------------------------------------------------------------*/

493: /*@
494:    MatCreateLMVMDFP - Creates a limited-memory Davidon-Fletcher-Powell (DFP) matrix 
495:    used for approximating Jacobians. L-DFP is symmetric positive-definite by 
496:    construction, and is the dual of L-BFGS where Y and S vectors swap roles.
497:    
498:    The provided local and global sizes must match the solution and function vectors 
499:    used with MatLMVMUpdate() and MatSolve(). The resulting L-DFP matrix will have 
500:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in 
501:    parallel. To use the L-DFP matrix with other vector types, the matrix must be 
502:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate(). 
503:    This ensures that the internal storage and work vectors are duplicated from the 
504:    correct type of vector.

506:    Collective on MPI_Comm

508:    Input Parameters:
509: +  comm - MPI communicator, set to PETSC_COMM_SELF
510: .  n - number of local rows for storage vectors
511: -  N - global size of the storage vectors

513:    Output Parameter:
514: .  B - the matrix

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

519:    Options Database Keys:
520: .   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
521: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
522: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
523: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
524: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
525: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
526: .   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

528:    Level: intermediate

530: .seealso: MatCreate(), MATLMVM, MATLMVMDFP, MatCreateLMVMBFGS(), MatCreateLMVMSR1(), 
531:            MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn(), MatCreateLMVMSymBrdn()
532: @*/
533: PetscErrorCode MatCreateLMVMDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
534: {
535:   PetscErrorCode    ierr;
536: 
538:   MatCreate(comm, B);
539:   MatSetSizes(*B, n, n, N, N);
540:   MatSetType(*B, MATLMVMDFP);
541:   MatSetUp(*B);
542:   return(0);
543: }