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:   PetscFunctionBegin;
 11:   PetscCall(VecPointwiseMult(dX, ldb->invD, F));
 12:   PetscFunctionReturn(PETSC_SUCCESS);
 13: }

 15: /*------------------------------------------------------------*/

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

 22:   PetscFunctionBegin;
 23:   PetscCall(VecPointwiseDivide(Z, X, ldb->invD));
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 25: }

 27: /*------------------------------------------------------------*/

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

 38:   PetscFunctionBegin;
 39:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
 40:   if (lmvm->prev_set) {
 41:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
 42:     PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
 43:     PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));

 45:     /* Test if the updates can be accepted */
 46:     PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ststmp));
 47:     if (ststmp < lmvm->eps) curvtol = 0.0;
 48:     else curvtol = lmvm->eps * ststmp;

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

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

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

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

 84:         if (1.0 != ldb->theta) {
 85:           /*  BFGS portion of the update */
 86:           /*  U = (inv(D)*s)*(inv(D)*s) */
 87:           PetscCall(VecPointwiseMult(ldb->U, ldb->W, ldb->W));

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

 97:           /*  Assemble */
 98:           PetscCall(VecAXPBY(ldb->DFP, stDs / ldb->yts[lmvm->k], 0.0, ldb->V));
 99:           PetscCall(VecAXPY(ldb->DFP, -2.0, ldb->U));
100:         }

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

111:         PetscCall(VecAXPY(ldb->invDnew, 1.0 / ldb->yts[lmvm->k], ldb->V));
112:         /*  Obtain inverse and ensure positive definite */
113:         PetscCall(VecReciprocal(ldb->invDnew));
114:         PetscCall(VecAbs(ldb->invDnew));

116:       } else {
117:         /* Inverse Hessian update instead. */
118:         PetscCall(VecCopy(ldb->invD, ldb->invDnew));

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

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

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

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

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

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

145:           /*  Assemble */
146:           PetscCall(VecAXPBY(ldb->DFP, -1.0 / ytDy, 0.0, ldb->U));
147:         }

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

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

171:             ss_sum = PetscRealPart(stDs);
172:             yy_sum = PetscRealPart(ytDy);
173:             ys_sum = ldb->yts[0];
174:           } else {
175:             PetscCall(VecCopy(ldb->invDnew, ldb->U));
176:             PetscCall(VecReciprocal(ldb->U));

178:             /*  Compute summations for scalar scaling */
179:             yy_sum = 0; /*  No safeguard required */
180:             ys_sum = 0; /*  No safeguard required */
181:             ss_sum = 0; /*  No safeguard required */
182:             start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
183:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
184:               PetscCall(VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U));
185:               PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));

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

190:               ss_sum += PetscRealPart(stDs);
191:               ys_sum += ldb->yts[i];
192:               yy_sum += PetscRealPart(ytDy);
193:             }
194:           }
195:         } else if (0.0 == ldb->beta) {
196:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
197:             /*  Compute summations for scalar scaling */
198:             PetscCall(VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew));

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

202:             ys_sum = PetscRealPart(ytDs);
203:             ss_sum = stDsr;
204:             yy_sum = ldb->yty[0];
205:           } else {
206:             PetscCall(VecCopy(ldb->invDnew, ldb->U));
207:             PetscCall(VecReciprocal(ldb->U));

209:             /*  Compute summations for scalar scaling */
210:             yy_sum = 0; /*  No safeguard required */
211:             ys_sum = 0; /*  No safeguard required */
212:             ss_sum = 0; /*  No safeguard required */
213:             start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
214:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
215:               PetscCall(VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U));

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

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

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

235:             yy_sum += ytDyr;
236:             ys_sum += PetscRealPart(ytDs);
237:             ss_sum += ldb->sts[i];
238:           }
239:         } else {
240:           PetscCall(VecCopy(ldb->invDnew, ldb->U));
241:           PetscCall(VecPow(ldb->U, ldb->beta - 1));

243:           /*  Compute summations for scalar scaling */
244:           yy_sum = 0; /*  No safeguard required */
245:           ys_sum = 0; /*  No safeguard required */
246:           ss_sum = 0; /*  No safeguard required */
247:           start  = PetscMax(0, lmvm->k - ldb->sigma_hist + 1);
248:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
249:             PetscCall(VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]));
250:             PetscCall(VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]));

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

255:             yy_sum += ytDyr;
256:             ys_sum += PetscRealPart(ytDs);
257:             ss_sum += PetscRealPart(stDs);
258:           }
259:         }

261:         if (0.0 == ldb->alpha) {
262:           /*  Safeguard ys_sum  */
263:           ys_sum = PetscMax(ldb->tol, ys_sum);

265:           sigma = ss_sum / ys_sum;
266:         } else if (1.0 == ldb->alpha) {
267:           /* yy_sum is never 0; if it were, we'd be at the minimum */
268:           sigma = ys_sum / yy_sum;
269:         } else {
270:           denom = 2.0 * ldb->alpha * yy_sum;

272:           /*  Safeguard denom */
273:           denom = PetscMax(ldb->tol, denom);

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

285:       if (PetscIsInfOrNanScalar(sigma)) {
286:         /*  sigma is not-a-number; skip rescaling */
287:       } else if (0.0 == sigma) {
288:         /*  sigma is zero; this is a bad case; skip rescaling */
289:       } else {
290:         /*  sigma is positive */
291:         PetscCall(VecScale(ldb->invDnew, sigma));
292:       }

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

310: /*------------------------------------------------------------*/

312: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
313: {
314:   Mat_LMVM     *bdata = (Mat_LMVM *)B->data;
315:   Mat_DiagBrdn *bctx  = (Mat_DiagBrdn *)bdata->ctx;
316:   Mat_LMVM     *mdata = (Mat_LMVM *)M->data;
317:   Mat_DiagBrdn *mctx  = (Mat_DiagBrdn *)mdata->ctx;
318:   PetscInt      i;

320:   PetscFunctionBegin;
321:   mctx->theta      = bctx->theta;
322:   mctx->alpha      = bctx->alpha;
323:   mctx->beta       = bctx->beta;
324:   mctx->rho        = bctx->rho;
325:   mctx->delta      = bctx->delta;
326:   mctx->delta_min  = bctx->delta_min;
327:   mctx->delta_max  = bctx->delta_max;
328:   mctx->tol        = bctx->tol;
329:   mctx->sigma      = bctx->sigma;
330:   mctx->sigma_hist = bctx->sigma_hist;
331:   mctx->forward    = bctx->forward;
332:   PetscCall(VecCopy(bctx->invD, mctx->invD));
333:   for (i = 0; i <= bdata->k; ++i) {
334:     mctx->yty[i] = bctx->yty[i];
335:     mctx->yts[i] = bctx->yts[i];
336:     mctx->sts[i] = bctx->sts[i];
337:   }
338:   PetscFunctionReturn(PETSC_SUCCESS);
339: }

341: /*------------------------------------------------------------*/

343: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
344: {
345:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
346:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;
347:   PetscBool     isascii;

349:   PetscFunctionBegin;
350:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
351:   if (isascii) {
352:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", ldb->sigma_hist));
353:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho));
354:     PetscCall(PetscViewerASCIIPrintf(pv, "Convex factor: theta=%g\n", (double)ldb->theta));
355:   }
356:   PetscCall(MatView_LMVM(B, pv));
357:   PetscFunctionReturn(PETSC_SUCCESS);
358: }

360: /*------------------------------------------------------------*/

362: static PetscErrorCode MatSetFromOptions_DiagBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
363: {
364:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
365:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

367:   PetscFunctionBegin;
368:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
369:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
370:   PetscCall(PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL));
371:   PetscCall(PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL));
372:   PetscCall(PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL));
373:   PetscCall(PetscOptionsReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL));
374:   PetscCall(PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL));
375:   PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL));
376:   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));
377:   PetscOptionsHeadEnd();
378:   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]");
379:   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]");
380:   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]");
381:   PetscCheck(ldb->sigma_hist >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
382:   PetscFunctionReturn(PETSC_SUCCESS);
383: }

385: /*------------------------------------------------------------*/

387: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
388: {
389:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
390:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

392:   PetscFunctionBegin;
393:   PetscCall(VecSet(ldb->invD, ldb->delta));
394:   if (destructive && ldb->allocated) {
395:     PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
396:     PetscCall(VecDestroy(&ldb->invDnew));
397:     PetscCall(VecDestroy(&ldb->invD));
398:     PetscCall(VecDestroy(&ldb->BFGS));
399:     PetscCall(VecDestroy(&ldb->DFP));
400:     PetscCall(VecDestroy(&ldb->U));
401:     PetscCall(VecDestroy(&ldb->V));
402:     PetscCall(VecDestroy(&ldb->W));
403:     ldb->allocated = PETSC_FALSE;
404:   }
405:   PetscCall(MatReset_LMVM(B, destructive));
406:   PetscFunctionReturn(PETSC_SUCCESS);
407: }

409: /*------------------------------------------------------------*/

411: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
412: {
413:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
414:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

416:   PetscFunctionBegin;
417:   PetscCall(MatAllocate_LMVM(B, X, F));
418:   if (!ldb->allocated) {
419:     PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
420:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
421:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
422:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
423:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
424:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
425:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
426:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
427:     ldb->allocated = PETSC_TRUE;
428:   }
429:   PetscFunctionReturn(PETSC_SUCCESS);
430: }

432: /*------------------------------------------------------------*/

434: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
435: {
436:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
437:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

439:   PetscFunctionBegin;
440:   if (ldb->allocated) {
441:     PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
442:     PetscCall(VecDestroy(&ldb->invDnew));
443:     PetscCall(VecDestroy(&ldb->invD));
444:     PetscCall(VecDestroy(&ldb->BFGS));
445:     PetscCall(VecDestroy(&ldb->DFP));
446:     PetscCall(VecDestroy(&ldb->U));
447:     PetscCall(VecDestroy(&ldb->V));
448:     PetscCall(VecDestroy(&ldb->W));
449:     ldb->allocated = PETSC_FALSE;
450:   }
451:   PetscCall(PetscFree(lmvm->ctx));
452:   PetscCall(MatDestroy_LMVM(B));
453:   PetscFunctionReturn(PETSC_SUCCESS);
454: }

456: /*------------------------------------------------------------*/

458: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
459: {
460:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
461:   Mat_DiagBrdn *ldb  = (Mat_DiagBrdn *)lmvm->ctx;

463:   PetscFunctionBegin;
464:   PetscCall(MatSetUp_LMVM(B));
465:   if (!ldb->allocated) {
466:     PetscCall(PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts));
467:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
468:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invD));
469:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
470:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
471:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
472:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
473:     PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
474:     ldb->allocated = PETSC_TRUE;
475:   }
476:   PetscFunctionReturn(PETSC_SUCCESS);
477: }

479: /*------------------------------------------------------------*/

481: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
482: {
483:   Mat_LMVM     *lmvm;
484:   Mat_DiagBrdn *ldb;

486:   PetscFunctionBegin;
487:   PetscCall(MatCreate_LMVM(B));
488:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN));
489:   B->ops->setup          = MatSetUp_DiagBrdn;
490:   B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
491:   B->ops->destroy        = MatDestroy_DiagBrdn;
492:   B->ops->solve          = MatSolve_DiagBrdn;
493:   B->ops->view           = MatView_DiagBrdn;

495:   lmvm                = (Mat_LMVM *)B->data;
496:   lmvm->square        = PETSC_TRUE;
497:   lmvm->m             = 1;
498:   lmvm->ops->allocate = MatAllocate_DiagBrdn;
499:   lmvm->ops->reset    = MatReset_DiagBrdn;
500:   lmvm->ops->mult     = MatMult_DiagBrdn;
501:   lmvm->ops->update   = MatUpdate_DiagBrdn;
502:   lmvm->ops->copy     = MatCopy_DiagBrdn;

504:   PetscCall(PetscNew(&ldb));
505:   lmvm->ctx       = (void *)ldb;
506:   ldb->theta      = 0.0;
507:   ldb->alpha      = 1.0;
508:   ldb->rho        = 1.0;
509:   ldb->forward    = PETSC_TRUE;
510:   ldb->beta       = 0.5;
511:   ldb->sigma      = 1.0;
512:   ldb->delta      = 1.0;
513:   ldb->delta_min  = 1e-7;
514:   ldb->delta_max  = 100.0;
515:   ldb->tol        = 1e-8;
516:   ldb->sigma_hist = 1;
517:   ldb->allocated  = PETSC_FALSE;
518:   PetscFunctionReturn(PETSC_SUCCESS);
519: }

521: /*------------------------------------------------------------*/

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

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

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

541:    Collective

543:    Input Parameters:
544: +  comm - MPI communicator
545: .  n - number of local rows for storage vectors
546: -  N - global size of the storage vectors

548:    Output Parameter:
549: .  B - the matrix

551:    Options Database Keys:
552: +   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
553: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
554: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
555: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
556: .   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
557: .   -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
558: -   -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal

560:    Level: intermediate

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

566: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDIAGBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
567:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMSymBrdn()`
568: @*/
569: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
570: {
571:   PetscFunctionBegin;
572:   PetscCall(MatCreate(comm, B));
573:   PetscCall(MatSetSizes(*B, n, n, N, N));
574:   PetscCall(MatSetType(*B, MATLMVMDIAGBROYDEN));
575:   PetscCall(MatSetUp(*B));
576:   PetscFunctionReturn(PETSC_SUCCESS);
577: }