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) VecAXPBY(ldb->invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew);
319:     } else {
320:       MatLMVMReset(B, PETSC_FALSE);
321:     }
322:     /* End DiagBrdn update */
323:   }
324:   /* Save the solution and function to be used in the next update */
325:   VecCopy(X, lmvm->Xprev);
326:   VecCopy(F, lmvm->Fprev);
327:   lmvm->prev_set = PETSC_TRUE;
328:   return 0;
329: }

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

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

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

361: /*------------------------------------------------------------*/

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

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

379: /*------------------------------------------------------------*/

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

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

403: /*------------------------------------------------------------*/

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

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

426: /*------------------------------------------------------------*/

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

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

448: /*------------------------------------------------------------*/

450: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
451: {
452:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
453:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

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

471: /*------------------------------------------------------------*/

473: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
474: {
475:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
476:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

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

493: /*------------------------------------------------------------*/

495: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
496: {
497:   Mat_LMVM     *lmvm;
498:   Mat_DiagBrdn *ldb;

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

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

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

534: /*------------------------------------------------------------*/

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

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

549:    In order to use the DiagBrdn matrix with other vector types, i.e. doing matrix-vector products
550:    and matrix solves, the matrix must first be created using `MatCreate()` and `MatSetType()`,
551:    followed by `MatLMVMAllocate()`. Then it will be available for updating
552:    (via `MatLMVMUpdate()`) in one's favored solver implementation.

554:    Collective

556:    Input Parameters:
557: +  comm - MPI communicator
558: .  n - number of local rows for storage vectors
559: -  N - global size of the storage vectors

561:    Output Parameter:
562: .  B - the matrix

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

573:    Level: intermediate

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

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