Actual source code: diagbrdn.c

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

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

  8:   PetscFunctionBegin;
  9:   PetscCall(VecPointwiseMult(dX, ldb->invD, F));
 10:   PetscFunctionReturn(PETSC_SUCCESS);
 11: }

 13: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
 14: {
 15:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
 16:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

 18:   PetscFunctionBegin;
 19:   PetscCall(VecPointwiseDivide(Z, X, ldb->invD));
 20:   PetscFunctionReturn(PETSC_SUCCESS);
 21: }

 23: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
 24: {
 25:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
 26:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;
 27:   PetscInt      old_k, i, start;
 28:   PetscScalar   yty, curvature, ytDy, stDs, ytDs;
 29:   PetscReal     curvtol, sigma, yy_sum, ss_sum, ys_sum, denom, ststmp;
 30:   PetscReal     stDsr, ytDyr;

 32:   PetscFunctionBegin;
 33:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
 34:   if (lmvm->prev_set) {
 35:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
 36:     PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
 37:     PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));

 39:     /* Test if the updates can be accepted */
 40:     PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ststmp));
 41:     if (ststmp < lmvm->eps) curvtol = 0.0;
 42:     else curvtol = lmvm->eps * ststmp;

 44:     /* Test the curvature for the update */
 45:     if (PetscRealPart(curvature) > curvtol) {
 46:       /* Update is good so we accept it */
 47:       old_k = lmvm->k;
 48:       PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
 49:       /* If we hit the memory limit, shift the yty and yts arrays */
 50:       if (old_k == lmvm->k) {
 51:         for (i = 0; i <= lmvm->k - 1; ++i) {
 52:           ldb->yty[i] = ldb->yty[i + 1];
 53:           ldb->yts[i] = ldb->yts[i + 1];
 54:           ldb->sts[i] = ldb->sts[i + 1];
 55:         }
 56:       }
 57:       /* Accept dot products into the history */
 58:       PetscCall(VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty));
 59:       ldb->yty[lmvm->k] = PetscRealPart(yty);
 60:       ldb->yts[lmvm->k] = PetscRealPart(curvature);
 61:       ldb->sts[lmvm->k] = ststmp;
 62:       if (ldb->forward) {
 63:         /* We are doing diagonal scaling of the forward Hessian B */
 64:         /*  BFGS = DFP = inv(D); */
 65:         PetscCall(VecCopy(ldb->invD, ldb->invDnew));
 66:         PetscCall(VecReciprocal(ldb->invDnew));

 68:         /*  V = y*y */
 69:         PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]));

 71:         /*  W = inv(D)*s */
 72:         PetscCall(VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]));
 73:         PetscCall(VecDot(ldb->W, lmvm->S[lmvm->k], &stDs));

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

 78:         if (1.0 != ldb->theta) {
 79:           /*  BFGS portion of the update */
 80:           /*  U = (inv(D)*s)*(inv(D)*s) */
 81:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));

 83:           /*  Assemble */
 84:           PetscCall(VecAXPBY(ldb->BFGS, -1.0 / stDs, 0.0, ldb->U));
 85:         }
 86:         if (0.0 != ldb->theta) {
 87:           /*  DFP portion of the update */
 88:           /*  U = inv(D)*s*y */
 89:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]));

 91:           /*  Assemble */
 92:           PetscCall(VecAXPBY(ldb->DFP, stDs / ldb->yts[lmvm->k], 0.0, ldb->V));
 93:           PetscCall(VecAXPY(ldb->DFP, -2.0, ldb->U));
 94:         }

 96:         if (0.0 == ldb->theta) {
 97:           PetscCall(VecAXPY(ldb->invDnew, 1.0, ldb->BFGS));
 98:         } else if (1.0 == ldb->theta) {
 99:           PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->DFP));
100:         } else {
101:           /*  Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
102:           PetscCall(VecAXPBYPCZ(ldb->invDnew, 1.0 - ldb->theta, (ldb->theta) / ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP));
103:         }

105:         PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
106:         /*  Obtain inverse and ensure positive definite */
107:         PetscCall(VecReciprocal(ldb->invDnew));
108:         PetscCall(VecAbs(ldb->invDnew));

110:       } else {
111:         /* Inverse Hessian update instead. */
112:         PetscCall(VecCopy(ldb->invD, ldb->invDnew));

114:         /*  V = s*s */
115:         PetscCall(VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]));

117:         /*  W = D*y */
118:         PetscCall(VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]));
119:         PetscCall(VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy));

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

124:         if (1.0 != ldb->theta) {
125:           /*  BFGS portion of the update */
126:           /*  U = s*Dy */
127:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]));

129:           /*  Assemble */
130:           PetscCall(VecAXPBY(ldb->BFGS, ytDy / ldb->yts[lmvm->k], 0.0, ldb->V));
131:           PetscCall(VecAXPY(ldb->BFGS, -2.0, ldb->U));
132:         }
133:         if (0.0 != ldb->theta) {
134:           /*  DFP portion of the update */

136:           /*  U = (inv(D)*y)*(inv(D)*y) */
137:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));

139:           /*  Assemble */
140:           PetscCall(VecAXPBY(ldb->DFP, -1.0 / ytDy, 0.0, ldb->U));
141:         }

143:         if (0.0 == ldb->theta) {
144:           PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->BFGS));
145:         } else if (1.0 == ldb->theta) {
146:           PetscCall(VecAXPY(ldb->invDnew, 1.0, ldb->DFP));
147:         } else {
148:           /*  Broyden update U=(1-theta)*P + theta*Q */
149:           PetscCall(VecAXPBYPCZ(ldb->invDnew, (1.0 - ldb->theta) / ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP));
150:         }
151:         PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
152:         /*  Ensure positive definite */
153:         PetscCall(VecAbs(ldb->invDnew));
154:       }
155:       if (ldb->sigma_hist > 0) {
156:         /*  Start with re-scaling on the newly computed diagonal */
157:         if (0.5 == ldb->beta) {
158:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
159:             PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew));
160:             PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));

162:             PetscCall(VecDot(ldb->V, lmvm->Y[0], &ytDy));
163:             PetscCall(VecDot(ldb->W, lmvm->S[0], &stDs));

165:             ss_sum = PetscRealPart(stDs);
166:             yy_sum = PetscRealPart(ytDy);
167:             ys_sum = ldb->yts[0];
168:           } else {
169:             PetscCall(VecCopy(ldb->invDnew, ldb->U));
170:             PetscCall(VecReciprocal(ldb->U));

172:             /*  Compute summations for scalar scaling */
173:             yy_sum = 0; /*  No safeguard required */
174:             ys_sum = 0; /*  No safeguard required */
175:             ss_sum = 0; /*  No safeguard required */
176:             start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
177:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
178:               PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U));
179:               PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));

181:               PetscCall(VecDot(ldb->W, lmvm->S[i], &stDs));
182:               PetscCall(VecDot(ldb->V, lmvm->Y[i], &ytDy));

184:               ss_sum += PetscRealPart(stDs);
185:               ys_sum += ldb->yts[i];
186:               yy_sum += PetscRealPart(ytDy);
187:             }
188:           }
189:         } else if (0.0 == ldb->beta) {
190:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
191:             /*  Compute summations for scalar scaling */
192:             PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));

194:             PetscCall(VecDotNorm2(lmvm->Y[0], ldb->W, &ytDs, &stDsr));

196:             ys_sum = PetscRealPart(ytDs);
197:             ss_sum = stDsr;
198:             yy_sum = ldb->yty[0];
199:           } else {
200:             PetscCall(VecCopy(ldb->invDnew, ldb->U));
201:             PetscCall(VecReciprocal(ldb->U));

203:             /*  Compute summations for scalar scaling */
204:             yy_sum = 0; /*  No safeguard required */
205:             ys_sum = 0; /*  No safeguard required */
206:             ss_sum = 0; /*  No safeguard required */
207:             start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
208:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
209:               PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));

211:               PetscCall(VecDotNorm2(lmvm->Y[i], ldb->W, &ytDs, &stDsr));

213:               ss_sum += stDsr;
214:               ys_sum += PetscRealPart(ytDs);
215:               yy_sum += ldb->yty[i];
216:             }
217:           }
218:         } else if (1.0 == ldb->beta) {
219:           /*  Compute summations for scalar scaling */
220:           yy_sum = 0; /*  No safeguard required */
221:           ys_sum = 0; /*  No safeguard required */
222:           ss_sum = 0; /*  No safeguard required */
223:           start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
224:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
225:             PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew));

227:             PetscCall(VecDotNorm2(lmvm->S[i], ldb->V, &ytDs, &ytDyr));

229:             yy_sum += ytDyr;
230:             ys_sum += PetscRealPart(ytDs);
231:             ss_sum += ldb->sts[i];
232:           }
233:         } else {
234:           PetscCall(VecCopy(ldb->invDnew, ldb->U));
235:           PetscCall(VecPow(ldb->U, ldb->beta - 1));

237:           /*  Compute summations for scalar scaling */
238:           yy_sum = 0; /*  No safeguard required */
239:           ys_sum = 0; /*  No safeguard required */
240:           ss_sum = 0; /*  No safeguard required */
241:           start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
242:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
243:             PetscCall(VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]));
244:             PetscCall(VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]));

246:             PetscCall(VecDotNorm2(ldb->W, ldb->V, &ytDs, &ytDyr));
247:             PetscCall(VecDot(ldb->W, ldb->W, &stDs));

249:             yy_sum += ytDyr;
250:             ys_sum += PetscRealPart(ytDs);
251:             ss_sum += PetscRealPart(stDs);
252:           }
253:         }

255:         if (0.0 == ldb->alpha) {
256:           /*  Safeguard ys_sum  */
257:           ys_sum = PetscMax(ldb->tol, ys_sum);

259:           sigma = ss_sum / ys_sum;
260:         } else if (1.0 == ldb->alpha) {
261:           /* yy_sum is never 0; if it were, we'd be at the minimum */
262:           sigma = ys_sum / yy_sum;
263:         } else {
264:           denom = 2.0 * ldb->alpha * yy_sum;

266:           /*  Safeguard denom */
267:           denom = PetscMax(ldb->tol, denom);

269:           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;
270:         }
271:       } else {
272:         sigma = 1.0;
273:       }
274:       /*  If Q has small values, then Q^(r_beta - 1)
275:        can have very large values.  Hence, ys_sum
276:        and ss_sum can be infinity.  In this case,
277:        sigma can either be not-a-number or infinity. */

279:       if (PetscIsInfOrNanScalar(sigma)) {
280:         /*  sigma is not-a-number; skip rescaling */
281:       } else if (0.0 == sigma) {
282:         /*  sigma is zero; this is a bad case; skip rescaling */
283:       } else {
284:         /*  sigma is positive */
285:         PetscCall(VecScale(ldb->invDnew, sigma));
286:       }

288:       /* Combine the old diagonal and the new diagonal using a convex limiter */
289:       if (1.0 == ldb->rho) {
290:         PetscCall(VecCopy(ldb->invDnew, ldb->invD));
291:       } else if (ldb->rho) PetscCall(VecAXPBY(ldb->invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew));
292:     } else {
293:       PetscCall(MatLMVMReset(B, PETSC_FALSE));
294:     }
295:     /* End DiagBrdn update */
296:   }
297:   /* Save the solution and function to be used in the next update */
298:   PetscCall(VecCopy(X, lmvm->Xprev));
299:   PetscCall(VecCopy(F, lmvm->Fprev));
300:   lmvm->prev_set = PETSC_TRUE;
301:   PetscFunctionReturn(PETSC_SUCCESS);
302: }

304: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
305: {
306:   Mat_LMVM     *bdata = (Mat_LMVM *)B->data;
307:   Mat_DiagBrdn *bctx  = (Mat_DiagBrdn *)bdata->ctx;
308:   Mat_LMVM     *mdata = (Mat_LMVM *)M->data;
309:   Mat_DiagBrdn *mctx  = (Mat_DiagBrdn *)mdata->ctx;
310:   PetscInt      i;

312:   PetscFunctionBegin;
313:   mctx->theta      = bctx->theta;
314:   mctx->alpha      = bctx->alpha;
315:   mctx->beta       = bctx->beta;
316:   mctx->rho        = bctx->rho;
317:   mctx->delta      = bctx->delta;
318:   mctx->delta_min  = bctx->delta_min;
319:   mctx->delta_max  = bctx->delta_max;
320:   mctx->tol        = bctx->tol;
321:   mctx->sigma      = bctx->sigma;
322:   mctx->sigma_hist = bctx->sigma_hist;
323:   mctx->forward    = bctx->forward;
324:   PetscCall(VecCopy(bctx->invD, mctx->invD));
325:   for (i = 0; i <= bdata->k; ++i) {
326:     mctx->yty[i] = bctx->yty[i];
327:     mctx->yts[i] = bctx->yts[i];
328:     mctx->sts[i] = bctx->sts[i];
329:   }
330:   PetscFunctionReturn(PETSC_SUCCESS);
331: }

333: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
334: {
335:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
336:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;
337:   PetscBool     isascii;

339:   PetscFunctionBegin;
340:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
341:   if (isascii) {
342:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", ldb->sigma_hist));
343:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho));
344:     PetscCall(PetscViewerASCIIPrintf(pv, "Convex factor: theta=%g\n", (double)ldb->theta));
345:   }
346:   PetscCall(MatView_LMVM(B, pv));
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: static PetscErrorCode MatSetFromOptions_DiagBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
351: {
352:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
353:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

355:   PetscFunctionBegin;
356:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
357:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
358:   PetscCall(PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL));
359:   PetscCall(PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL));
360:   PetscCall(PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL));
361:   PetscCall(PetscOptionsReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL));
362:   PetscCall(PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL));
363:   PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL));
364:   PetscCall(PetscOptionsInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", ldb->sigma_hist, &ldb->sigma_hist, NULL));
365:   PetscOptionsHeadEnd();
366:   PetscCheck(!(ldb->theta < 0.0) && !(ldb->theta > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
367:   PetscCheck(!(ldb->alpha < 0.0) && !(ldb->alpha > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
368:   PetscCheck(!(ldb->rho < 0.0) && !(ldb->rho > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex update limiter in the J0 scaling cannot be outside the range of [0, 1]");
369:   PetscCheck(ldb->sigma_hist >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
370:   PetscFunctionReturn(PETSC_SUCCESS);
371: }

373: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
374: {
375:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
376:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

378:   PetscFunctionBegin;
379:   PetscCall(VecSet(ldb->invD, ldb->delta));
380:   if (destructive && ldb->allocated) {
381:     PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
382:     PetscCall(VecDestroy(&ldb->invDnew));
383:     PetscCall(VecDestroy(&ldb->invD));
384:     PetscCall(VecDestroy(&ldb->BFGS));
385:     PetscCall(VecDestroy(&ldb->DFP));
386:     PetscCall(VecDestroy(&ldb->U));
387:     PetscCall(VecDestroy(&ldb->V));
388:     PetscCall(VecDestroy(&ldb->W));
389:     ldb->allocated = PETSC_FALSE;
390:   }
391:   PetscCall(MatReset_LMVM(B, destructive));
392:   PetscFunctionReturn(PETSC_SUCCESS);
393: }

395: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
396: {
397:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
398:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

400:   PetscFunctionBegin;
401:   PetscCall(MatAllocate_LMVM(B, X, F));
402:   if (!ldb->allocated) {
403:     PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
404:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
405:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
406:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
407:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
408:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
409:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
410:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
411:     ldb->allocated = PETSC_TRUE;
412:   }
413:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
417: {
418:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
419:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

421:   PetscFunctionBegin;
422:   if (ldb->allocated) {
423:     PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
424:     PetscCall(VecDestroy(&ldb->invDnew));
425:     PetscCall(VecDestroy(&ldb->invD));
426:     PetscCall(VecDestroy(&ldb->BFGS));
427:     PetscCall(VecDestroy(&ldb->DFP));
428:     PetscCall(VecDestroy(&ldb->U));
429:     PetscCall(VecDestroy(&ldb->V));
430:     PetscCall(VecDestroy(&ldb->W));
431:     ldb->allocated = PETSC_FALSE;
432:   }
433:   PetscCall(PetscFree(lmvm->ctx));
434:   PetscCall(MatDestroy_LMVM(B));
435:   PetscFunctionReturn(PETSC_SUCCESS);
436: }

438: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
439: {
440:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
441:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

443:   PetscFunctionBegin;
444:   PetscCall(MatSetUp_LMVM(B));
445:   if (!ldb->allocated) {
446:     PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
447:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
448:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
449:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
450:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
451:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
452:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
453:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
454:     ldb->allocated = PETSC_TRUE;
455:   }
456:   PetscFunctionReturn(PETSC_SUCCESS);
457: }

459: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
460: {
461:   Mat_LMVM     *lmvm;
462:   Mat_DiagBrdn *ldb;

464:   PetscFunctionBegin;
465:   PetscCall(MatCreate_LMVM(B));
466:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN));
467:   B->ops->setup          = MatSetUp_DiagBrdn;
468:   B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
469:   B->ops->destroy        = MatDestroy_DiagBrdn;
470:   B->ops->solve          = MatSolve_DiagBrdn;
471:   B->ops->view           = MatView_DiagBrdn;

473:   lmvm                = (Mat_LMVM *)B->data;
474:   lmvm->square        = PETSC_TRUE;
475:   lmvm->m             = 1;
476:   lmvm->ops->allocate = MatAllocate_DiagBrdn;
477:   lmvm->ops->reset    = MatReset_DiagBrdn;
478:   lmvm->ops->mult     = MatMult_DiagBrdn;
479:   lmvm->ops->update   = MatUpdate_DiagBrdn;
480:   lmvm->ops->copy     = MatCopy_DiagBrdn;

482:   PetscCall(PetscNew(&ldb));
483:   lmvm->ctx       = (void *)ldb;
484:   ldb->theta      = 0.0;
485:   ldb->alpha      = 1.0;
486:   ldb->rho        = 1.0;
487:   ldb->forward    = PETSC_TRUE;
488:   ldb->beta       = 0.5;
489:   ldb->sigma      = 1.0;
490:   ldb->delta      = 1.0;
491:   ldb->delta_min  = 1e-7;
492:   ldb->delta_max  = 100.0;
493:   ldb->tol        = 1e-8;
494:   ldb->sigma_hist = 1;
495:   ldb->allocated  = PETSC_FALSE;
496:   PetscFunctionReturn(PETSC_SUCCESS);
497: }

499: /*@
500:   MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
501:   for approximating Hessians.

503:   Collective

505:   Input Parameters:
506: + comm - MPI communicator
507: . n    - number of local rows for storage vectors
508: - N    - global size of the storage vectors

510:   Output Parameter:
511: . B - the matrix

513:   Options Database Keys:
514: + -mat_lmvm_theta      - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
515: . -mat_lmvm_rho        - (developer) update limiter for the J0 scaling
516: . -mat_lmvm_alpha      - (developer) coefficient factor for the quadratic subproblem in J0 scaling
517: . -mat_lmvm_beta       - (developer) exponential factor for the diagonal J0 scaling
518: . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
519: . -mat_lmvm_tol        - (developer) tolerance for bounding the denominator of the rescaling away from 0.
520: - -mat_lmvm_forward    - (developer) whether or not to use the forward or backward Broyden update to the diagonal

522:   Level: intermediate

524:   Notes:
525:   It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
526:   paradigm instead of this routine directly.

528:   It consists of a convex combination of DFP and BFGS
529:   diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
530:   To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
531:   We also ensure positive definiteness by taking the `VecAbs()` of the final vector.

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

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

544: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDIAGBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
545:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMSymBrdn()`
546: @*/
547: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
548: {
549:   PetscFunctionBegin;
550:   PetscCall(MatCreate(comm, B));
551:   PetscCall(MatSetSizes(*B, n, n, N, N));
552:   PetscCall(MatSetType(*B, MATLMVMDIAGBROYDEN));
553:   PetscCall(MatSetUp(*B));
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }