Actual source code: diagbrdn.c

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

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

  5: static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
  6: {
  7:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
  8:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;

 10:   VecCheckSameSize(F, 2, dX, 3);
 11:   VecCheckMatCompatible(B, dX, 3, F, 2);
 12:   VecPointwiseMult(dX, ldb->invD, F);
 13:   return 0;
 14: }

 16: /*------------------------------------------------------------*/

 18: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
 19: {
 20:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
 21:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;

 23:   VecCheckSameSize(X, 2, Z, 3);
 24:   VecCheckMatCompatible(B, X, 2, Z, 3);
 25:   VecPointwiseDivide(Z, X, ldb->invD);
 26:   return 0;
 27: }

 29: /*------------------------------------------------------------*/

 31: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
 32: {
 33:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
 34:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
 35:   PetscInt          old_k, i, start;
 36:   PetscScalar       yty, ststmp, curvature, ytDy, stDs, ytDs;
 37:   PetscReal         curvtol, sigma, yy_sum, ss_sum, ys_sum, denom;

 39:   if (!lmvm->m) return 0;
 40:   if (lmvm->prev_set) {
 41:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
 42:     VecAYPX(lmvm->Xprev, -1.0, X);
 43:     VecAYPX(lmvm->Fprev, -1.0, F);
 44:     /* Compute tolerance for accepting the update */
 45:     VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
 46:     VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
 47:     VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
 48:     VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
 49:     if (PetscRealPart(ststmp) < lmvm->eps) {
 50:       curvtol = 0.0;
 51:     } else {
 52:       curvtol = lmvm->eps * PetscRealPart(ststmp);
 53:     }
 54:     /* Test the curvature for the update */
 55:     if (PetscRealPart(curvature) > curvtol) {
 56:       /* Update is good so we accept it */
 57:       old_k = lmvm->k;
 58:       MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
 59:       /* If we hit the memory limit, shift the yty and yts arrays */
 60:       if (old_k == lmvm->k) {
 61:         for (i = 0; i <= lmvm->k-1; ++i) {
 62:           ldb->yty[i] = ldb->yty[i+1];
 63:           ldb->yts[i] = ldb->yts[i+1];
 64:           ldb->sts[i] = ldb->sts[i+1];
 65:         }
 66:       }
 67:       /* Accept dot products into the history */
 68:       VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
 69:       ldb->yty[lmvm->k] = PetscRealPart(yty);
 70:       ldb->yts[lmvm->k] = PetscRealPart(curvature);
 71:       ldb->sts[lmvm->k] = PetscRealPart(ststmp);
 72:       if (ldb->forward) {
 73:         /* We are doing diagonal scaling of the forward Hessian B */
 74:         /*  BFGS = DFP = inv(D); */
 75:         VecCopy(ldb->invD, ldb->invDnew);
 76:         VecReciprocal(ldb->invDnew);

 78:         /*  V = y*y */
 79:         VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]);

 81:         /*  W = inv(D)*s */
 82:         VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]);
 83:         VecDot(ldb->W, lmvm->S[lmvm->k], &stDs);

 85:         /*  Safeguard stDs */
 86:         stDs = PetscMax(PetscRealPart(stDs), ldb->tol);

 88:         if (1.0 != ldb->theta) {
 89:           /*  BFGS portion of the update */
 90:           /*  U = (inv(D)*s)*(inv(D)*s) */
 91:           VecPointwiseMult(ldb->U, ldb->W, ldb->W);

 93:           /*  Assemble */
 94:           VecAXPBY(ldb->BFGS, -1.0/stDs, 0.0, ldb->U);
 95:         }
 96:         if (0.0 != ldb->theta) {
 97:           /*  DFP portion of the update */
 98:           /*  U = inv(D)*s*y */
 99:           VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]);

101:           /*  Assemble */
102:           VecAXPBY(ldb->DFP, stDs/ldb->yts[lmvm->k], 0.0, ldb->V);
103:           VecAXPY(ldb->DFP, -2.0, ldb->U);
104:         }

106:         if (0.0 == ldb->theta) {
107:           VecAXPY(ldb->invDnew, 1.0, ldb->BFGS);
108:         } else if (1.0 == ldb->theta) {
109:           VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->DFP);
110:         } else {
111:           /*  Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
112:           VecAXPBYPCZ(ldb->invDnew, 1.0-ldb->theta, (ldb->theta)/ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP);
113:         }

115:         VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
116:         /*  Obtain inverse and ensure positive definite */
117:         VecReciprocal(ldb->invDnew);
118:         VecAbs(ldb->invDnew);

120:       } else {
121:         /* Inverse Hessian update instead. */
122:         VecCopy(ldb->invD, ldb->invDnew);

124:         /*  V = s*s */
125:         VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]);

127:         /*  W = D*y */
128:         VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]);
129:         VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy);

131:         /*  Safeguard ytDy */
132:         ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);

134:         if (1.0 != ldb->theta) {
135:           /*  BFGS portion of the update */
136:           /*  U = s*Dy */
137:           VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]);

139:           /*  Assemble */
140:           VecAXPBY(ldb->BFGS, ytDy/ldb->yts[lmvm->k], 0.0, ldb->V);
141:           VecAXPY(ldb->BFGS, -2.0, ldb->U);
142:         }
143:         if (0.0 != ldb->theta) {
144:           /*  DFP portion of the update */

146:           /*  U = (inv(D)*y)*(inv(D)*y) */
147:           VecPointwiseMult(ldb->U, ldb->W, ldb->W);

149:           /*  Assemble */
150:           VecAXPBY(ldb->DFP, -1.0/ytDy, 0.0, ldb->U);
151:         }

153:         if (0.0 == ldb->theta) {
154:           VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->BFGS);
155:         } else if (1.0 == ldb->theta) {
156:           VecAXPY(ldb->invDnew, 1.0, ldb->DFP);
157:         } else {
158:           /*  Broyden update U=(1-theta)*P + theta*Q */
159:           VecAXPBYPCZ(ldb->invDnew, (1.0-ldb->theta)/ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP);
160:         }
161:         VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
162:         /*  Ensure positive definite */
163:         VecAbs(ldb->invDnew);
164:       }
165:       if (ldb->sigma_hist > 0) {
166:         /*  Start with re-scaling on the newly computed diagonal */
167:         if (0.5 == ldb->beta) {
168:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
169:             VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew);
170:             VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);

172:             VecDotBegin(ldb->V, lmvm->Y[0], &ytDy);
173:             VecDotBegin(ldb->W, lmvm->S[0], &stDs);
174:             VecDotEnd(ldb->V, lmvm->Y[0], &ytDy);
175:             VecDotEnd(ldb->W, lmvm->S[0], &stDs);

177:             ss_sum = PetscRealPart(stDs);
178:             yy_sum = PetscRealPart(ytDy);
179:             ys_sum = ldb->yts[0];
180:           } else {
181:             VecCopy(ldb->invDnew, ldb->U);
182:             VecReciprocal(ldb->U);

184:             /*  Compute summations for scalar scaling */
185:             yy_sum = 0;       /*  No safeguard required */
186:             ys_sum = 0;       /*  No safeguard required */
187:             ss_sum = 0;       /*  No safeguard required */
188:             start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
189:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
190:               VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U);
191:               VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);

193:               VecDotBegin(ldb->W, lmvm->S[i], &stDs);
194:               VecDotBegin(ldb->V, lmvm->Y[i], &ytDy);
195:               VecDotEnd(ldb->W, lmvm->S[i], &stDs);
196:               VecDotEnd(ldb->V, lmvm->Y[i], &ytDy);

198:               ss_sum += PetscRealPart(stDs);
199:               ys_sum += ldb->yts[i];
200:               yy_sum += PetscRealPart(ytDy);
201:             }
202:           }
203:         } else if (0.0 == ldb->beta) {
204:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
205:             /*  Compute summations for scalar scaling */
206:             VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);

208:             VecDotBegin(ldb->W, lmvm->Y[0], &ytDs);
209:             VecDotBegin(ldb->W, ldb->W, &stDs);
210:             VecDotEnd(ldb->W, lmvm->Y[0], &ytDs);
211:             VecDotEnd(ldb->W, ldb->W, &stDs);

213:             ys_sum = PetscRealPart(ytDs);
214:             ss_sum = PetscRealPart(stDs);
215:             yy_sum = ldb->yty[0];
216:           } else {
217:             VecCopy(ldb->invDnew, ldb->U);
218:             VecReciprocal(ldb->U);

220:             /*  Compute summations for scalar scaling */
221:             yy_sum = 0;       /*  No safeguard required */
222:             ys_sum = 0;       /*  No safeguard required */
223:             ss_sum = 0;       /*  No safeguard required */
224:             start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
225:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
226:               VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);

228:               VecDotBegin(ldb->W, lmvm->Y[i], &ytDs);
229:               VecDotBegin(ldb->W, ldb->W, &stDs);
230:               VecDotEnd(ldb->W, lmvm->Y[i], &ytDs);
231:               VecDotEnd(ldb->W, ldb->W, &stDs);

233:               ss_sum += PetscRealPart(stDs);
234:               ys_sum += PetscRealPart(ytDs);
235:               yy_sum += ldb->yty[i];
236:             }
237:           }
238:         } else if (1.0 == ldb->beta) {
239:           /*  Compute summations for scalar scaling */
240:           yy_sum = 0; /*  No safeguard required */
241:           ys_sum = 0; /*  No safeguard required */
242:           ss_sum = 0; /*  No safeguard required */
243:           start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
244:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
245:             VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew);

247:             VecDotBegin(ldb->V, lmvm->S[i], &ytDs);
248:             VecDotBegin(ldb->V, ldb->V, &ytDy);
249:             VecDotEnd(ldb->V, lmvm->S[i], &ytDs);
250:             VecDotEnd(ldb->V, ldb->V, &ytDy);

252:             yy_sum += PetscRealPart(ytDy);
253:             ys_sum += PetscRealPart(ytDs);
254:             ss_sum += ldb->sts[i];
255:           }
256:         } else {
257:           VecCopy(ldb->invDnew, ldb->U);
258:           VecPow(ldb->U, ldb->beta-1);

260:           /*  Compute summations for scalar scaling */
261:           yy_sum = 0; /*  No safeguard required */
262:           ys_sum = 0; /*  No safeguard required */
263:           ss_sum = 0; /*  No safeguard required */
264:           start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
265:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
266:             VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]);
267:             VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]);

269:             VecDotBegin(ldb->V, ldb->W, &ytDs);
270:             VecDotBegin(ldb->V, ldb->V, &ytDy);
271:             VecDotBegin(ldb->W, ldb->W, &stDs);
272:             VecDotEnd(ldb->V, ldb->W, &ytDs);
273:             VecDotEnd(ldb->V, ldb->V, &ytDy);
274:             VecDotEnd(ldb->W, ldb->W, &stDs);

276:             yy_sum += PetscRealPart(ytDy);
277:             ys_sum += PetscRealPart(ytDs);
278:             ss_sum += PetscRealPart(stDs);
279:           }
280:         }

282:         if (0.0 == ldb->alpha) {
283:           /*  Safeguard ys_sum  */
284:           ys_sum = PetscMax(ldb->tol, ys_sum);

286:           sigma = ss_sum / ys_sum;
287:         } else if (1.0 == ldb->alpha) {
288:           /* yy_sum is never 0; if it were, we'd be at the minimum */
289:           sigma = ys_sum / yy_sum;
290:         } else {
291:           denom = 2.0*ldb->alpha*yy_sum;

293:           /*  Safeguard denom */
294:           denom = PetscMax(ldb->tol, denom);

296:           sigma = ((2.0*ldb->alpha-1)*ys_sum + PetscSqrtReal((2.0*ldb->alpha-1)*(2.0*ldb->alpha-1)*ys_sum*ys_sum - 4.0*ldb->alpha*(ldb->alpha-1)*yy_sum*ss_sum)) / denom;
297:         }
298:       } else {
299:         sigma = 1.0;
300:       }
301:       /*  If Q has small values, then Q^(r_beta - 1)
302:        can have very large values.  Hence, ys_sum
303:        and ss_sum can be infinity.  In this case,
304:        sigma can either be not-a-number or infinity. */

306:       if (PetscIsInfOrNanScalar(sigma)) {
307:         /*  sigma is not-a-number; skip rescaling */
308:       } else if (0.0 == sigma) {
309:         /*  sigma is zero; this is a bad case; skip rescaling */
310:       } else {
311:         /*  sigma is positive */
312:         VecScale(ldb->invDnew, sigma);
313:       }

315:       /* Combine the old diagonal and the new diagonal using a convex limiter */
316:       if (1.0 == ldb->rho) {
317:         VecCopy(ldb->invDnew, ldb->invD);
318:       } else if (ldb->rho) {
319:         VecAXPBY(ldb->invD, 1.0-ldb->rho, ldb->rho, ldb->invDnew);
320:       }
321:     } else {
322:       MatLMVMReset(B, PETSC_FALSE);
323:     }
324:     /* End DiagBrdn update */

326:   }
327:   /* Save the solution and function to be used in the next update */
328:   VecCopy(X, lmvm->Xprev);
329:   VecCopy(F, lmvm->Fprev);
330:   lmvm->prev_set = PETSC_TRUE;
331:   return 0;
332: }

334: /*------------------------------------------------------------*/

336: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
337: {
338:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
339:   Mat_DiagBrdn      *bctx = (Mat_DiagBrdn*)bdata->ctx;
340:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
341:   Mat_DiagBrdn      *mctx = (Mat_DiagBrdn*)mdata->ctx;
342:   PetscInt          i;

344:   mctx->theta = bctx->theta;
345:   mctx->alpha = bctx->alpha;
346:   mctx->beta = bctx->beta;
347:   mctx->rho = bctx->rho;
348:   mctx->delta = bctx->delta;
349:   mctx->delta_min = bctx->delta_min;
350:   mctx->delta_max = bctx->delta_max;
351:   mctx->tol = bctx->tol;
352:   mctx->sigma = bctx->sigma;
353:   mctx->sigma_hist = bctx->sigma_hist;
354:   mctx->forward = bctx->forward;
355:   VecCopy(bctx->invD, mctx->invD);
356:   for (i=0; i<=bdata->k; ++i) {
357:     mctx->yty[i] = bctx->yty[i];
358:     mctx->yts[i] = bctx->yts[i];
359:     mctx->sts[i] = bctx->sts[i];
360:   }
361:   return 0;
362: }

364: /*------------------------------------------------------------*/

366: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
367: {
368:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
369:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
370:   PetscBool         isascii;

372:   PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
373:   if (isascii) {
374:     PetscViewerASCIIPrintf(pv,"Scale history: %d\n",ldb->sigma_hist);
375:     PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)ldb->alpha, (double)ldb->beta, (double)ldb->rho);
376:     PetscViewerASCIIPrintf(pv,"Convex factor: theta=%g\n", (double)ldb->theta);
377:   }
378:   MatView_LMVM(B, pv);
379:   return 0;
380: }

382: /*------------------------------------------------------------*/

384: static PetscErrorCode MatSetFromOptions_DiagBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
385: {
386:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
387:   Mat_DiagBrdn       *ldb = (Mat_DiagBrdn*)lmvm->ctx;

389:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
390:   PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
391:   PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",ldb->theta,&ldb->theta,NULL);
392:   PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",ldb->rho,&ldb->rho,NULL);
393:   PetscOptionsReal("-mat_lmvm_tol","(developer) tolerance for bounding rescaling denominator","",ldb->tol,&ldb->tol,NULL);
394:   PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",ldb->alpha,&ldb->alpha,NULL);
395:   PetscOptionsBool("-mat_lmvm_forward","Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.","",ldb->forward,&ldb->forward,NULL);
396:   PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",ldb->beta,&ldb->beta,NULL);
397:   PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",ldb->sigma_hist,&ldb->sigma_hist,NULL);
398:   PetscOptionsTail();
403:   return 0;
404: }

406: /*------------------------------------------------------------*/

408: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
409: {
410:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
411:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;

413:   VecSet(ldb->invD, ldb->delta);
414:   if (destructive && ldb->allocated) {
415:     PetscFree3(ldb->yty, ldb->yts, ldb->sts);
416:     VecDestroy(&ldb->invDnew);
417:     VecDestroy(&ldb->invD);
418:     VecDestroy(&ldb->BFGS);
419:     VecDestroy(&ldb->DFP);
420:     VecDestroy(&ldb->U);
421:     VecDestroy(&ldb->V);
422:     VecDestroy(&ldb->W);
423:     ldb->allocated = PETSC_FALSE;
424:   }
425:   MatReset_LMVM(B, destructive);
426:   return 0;
427: }

429: /*------------------------------------------------------------*/

431: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
432: {
433:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
434:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;

436:   MatAllocate_LMVM(B, X, F);
437:   if (!ldb->allocated) {
438:     PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
439:     VecDuplicate(lmvm->Xprev, &ldb->invDnew);
440:     VecDuplicate(lmvm->Xprev, &ldb->invD);
441:     VecDuplicate(lmvm->Xprev, &ldb->BFGS);
442:     VecDuplicate(lmvm->Xprev, &ldb->DFP);
443:     VecDuplicate(lmvm->Xprev, &ldb->U);
444:     VecDuplicate(lmvm->Xprev, &ldb->V);
445:     VecDuplicate(lmvm->Xprev, &ldb->W);
446:     ldb->allocated = PETSC_TRUE;
447:   }
448:   return 0;
449: }

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

453: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
454: {
455:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
456:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;

458:   if (ldb->allocated) {
459:     PetscFree3(ldb->yty, ldb->yts, ldb->sts);
460:     VecDestroy(&ldb->invDnew);
461:     VecDestroy(&ldb->invD);
462:     VecDestroy(&ldb->BFGS);
463:     VecDestroy(&ldb->DFP);
464:     VecDestroy(&ldb->U);
465:     VecDestroy(&ldb->V);
466:     VecDestroy(&ldb->W);
467:     ldb->allocated = PETSC_FALSE;
468:   }
469:   PetscFree(lmvm->ctx);
470:   MatDestroy_LMVM(B);
471:   return 0;
472: }

474: /*------------------------------------------------------------*/

476: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
477: {
478:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
479:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;

481:   MatSetUp_LMVM(B);
482:   if (!ldb->allocated) {
483:     PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
484:     VecDuplicate(lmvm->Xprev, &ldb->invDnew);
485:     VecDuplicate(lmvm->Xprev, &ldb->invD);
486:     VecDuplicate(lmvm->Xprev, &ldb->BFGS);
487:     VecDuplicate(lmvm->Xprev, &ldb->DFP);
488:     VecDuplicate(lmvm->Xprev, &ldb->U);
489:     VecDuplicate(lmvm->Xprev, &ldb->V);
490:     VecDuplicate(lmvm->Xprev, &ldb->W);
491:     ldb->allocated = PETSC_TRUE;
492:   }
493:   return 0;
494: }

496: /*------------------------------------------------------------*/

498: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
499: {
500:   Mat_LMVM          *lmvm;
501:   Mat_DiagBrdn      *ldb;

503:   MatCreate_LMVM(B);
504:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN);
505:   B->ops->setup = MatSetUp_DiagBrdn;
506:   B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
507:   B->ops->destroy = MatDestroy_DiagBrdn;
508:   B->ops->solve = MatSolve_DiagBrdn;
509:   B->ops->view = MatView_DiagBrdn;

511:   lmvm = (Mat_LMVM*)B->data;
512:   lmvm->square = PETSC_TRUE;
513:   lmvm->m = 1;
514:   lmvm->ops->allocate = MatAllocate_DiagBrdn;
515:   lmvm->ops->reset = MatReset_DiagBrdn;
516:   lmvm->ops->mult = MatMult_DiagBrdn;
517:   lmvm->ops->update = MatUpdate_DiagBrdn;
518:   lmvm->ops->copy = MatCopy_DiagBrdn;

520:   PetscNewLog(B, &ldb);
521:   lmvm->ctx = (void*)ldb;
522:   ldb->theta = 0.0;
523:   ldb->alpha = 1.0;
524:   ldb->rho = 1.0;
525:   ldb->forward = PETSC_TRUE;
526:   ldb->beta = 0.5;
527:   ldb->sigma = 1.0;
528:   ldb->delta = 1.0;
529:   ldb->delta_min = 1e-7;
530:   ldb->delta_max = 100.0;
531:   ldb->tol = 1e-8;
532:   ldb->sigma_hist = 1;
533:   ldb->allocated = PETSC_FALSE;
534:   return 0;
535: }

537: /*------------------------------------------------------------*/

539: /*@
540:    MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
541:    for approximating Hessians. It consists of a convex combination of DFP and BFGS
542:    diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
543:    To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
544:    We also ensure positive definiteness by taking the VecAbs() of the final vector.

546:    There are two ways of approximating the diagonal: using the forward (B) update
547:    schemes for BFGS and DFP and then taking the inverse, or directly working with
548:    the inverse (H) update schemes for the BFGS and DFP updates, derived using the
549:    Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a
550:    parameter below.

552:    In order to use the DiagBrdn matrix with other vector types, i.e. doing MatMults
553:    and MatSolves, the matrix must first be created using MatCreate() and MatSetType(),
554:    followed by MatLMVMAllocate(). Then it will be available for updating
555:    (via MatLMVMUpdate) in one's favored solver implementation.
556:    This allows for MPI compatibility.

558:    Collective

560:    Input Parameters:
561: +  comm - MPI communicator, set to PETSC_COMM_SELF
562: .  n - number of local rows for storage vectors
563: -  N - global size of the storage vectors

565:    Output Parameter:
566: .  B - the matrix

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

571:    Options Database Keys:
572: +   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
573: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
574: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
575: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
576: .   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
577: .   -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
578: -   -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal

580:    Level: intermediate

582: .seealso: MatCreate(), MATLMVM, MATLMVMDIAGBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
583:           MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMSymBrdn()
584: @*/
585: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
586: {
587:   MatCreate(comm, B);
588:   MatSetSizes(*B, n, n, N, N);
589:   MatSetType(*B, MATLMVMDIAGBROYDEN);
590:   MatSetUp(*B);
591:   return 0;
592: }