Actual source code: symbrdn.c

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

  4: const char *const MatLMVMSymBroydenScaleTypes[] = {"NONE", "SCALAR", "DIAGONAL", "USER", "MatLMVMSymBrdnScaleType", "MAT_LMVM_SYMBROYDEN_SCALING_", NULL};

  6: /*------------------------------------------------------------*/

  8: /*
  9:   The solution method below is the matrix-free implementation of
 10:   Equation 8.6a in Dennis and More "Quasi-Newton Methods, Motivation
 11:   and Theory" (https://epubs.siam.org/doi/abs/10.1137/1019005).

 13:   Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
 14:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
 15:   repeated calls of MatSolve without incurring redundant computation.

 17:   dX <- J0^{-1} * F

 19:   for i=0,1,2,...,k
 20:     # Q[i] = (B_i)^T{-1} Y[i]

 22:     rho = 1.0 / (Y[i]^T S[i])
 23:     alpha = rho * (S[i]^T F)
 24:     zeta = 1.0 / (Y[i]^T Q[i])
 25:     gamma = zeta * (Y[i]^T dX)

 27:     dX <- dX - (gamma * Q[i]) + (alpha * Y[i])
 28:     W <- (rho * S[i]) - (zeta * Q[i])
 29:     dX <- dX + (psi[i] * (Y[i]^T Q[i]) * (W^T F) * W)
 30:   end
 31: */
 32: static PetscErrorCode MatSolve_LMVMSymBrdn(Mat B, Vec F, Vec dX)
 33: {
 34:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
 35:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
 36:   PetscInt     i, j;
 37:   PetscReal    numer;
 38:   PetscScalar  sjtpi, yjtsi, wtsi, yjtqi, sjtyi, wtyi, ytx, stf, wtf, stp, ytq;

 40:   PetscFunctionBegin;
 41:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
 42:   if (lsb->phi == 0.0) {
 43:     PetscCall(MatSolve_LMVMBFGS(B, F, dX));
 44:     PetscFunctionReturn(PETSC_SUCCESS);
 45:   }
 46:   if (lsb->phi == 1.0) {
 47:     PetscCall(MatSolve_LMVMDFP(B, F, dX));
 48:     PetscFunctionReturn(PETSC_SUCCESS);
 49:   }

 51:   VecCheckSameSize(F, 2, dX, 3);
 52:   VecCheckMatCompatible(B, dX, 3, F, 2);

 54:   if (lsb->needP) {
 55:     /* Start the loop for (P[k] = (B_k) * S[k]) */
 56:     for (i = 0; i <= lmvm->k; ++i) {
 57:       PetscCall(MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]));
 58:       /* Compute the necessary dot products */
 59:       PetscCall(VecMDot(lmvm->S[i], i, lmvm->Y, lsb->workscalar));
 60:       for (j = 0; j < i; ++j) {
 61:         yjtsi = lsb->workscalar[j];
 62:         PetscCall(VecDot(lmvm->S[j], lsb->P[i], &sjtpi));
 63:         /* Compute the pure BFGS component of the forward product */
 64:         PetscCall(VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi) / lsb->stp[j], PetscRealPart(yjtsi) / lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]));
 65:         /* Tack on the convexly scaled extras to the forward product */
 66:         if (lsb->phi > 0.0) {
 67:           PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]));
 68:           PetscCall(VecDot(lsb->work, lmvm->S[i], &wtsi));
 69:           PetscCall(VecAXPY(lsb->P[i], lsb->phi * lsb->stp[j] * PetscRealPart(wtsi), lsb->work));
 70:         }
 71:       }
 72:       PetscCall(VecDot(lmvm->S[i], lsb->P[i], &stp));
 73:       lsb->stp[i] = PetscRealPart(stp);
 74:     }
 75:     lsb->needP = PETSC_FALSE;
 76:   }
 77:   if (lsb->needQ) {
 78:     /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
 79:     for (i = 0; i <= lmvm->k; ++i) {
 80:       PetscCall(MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]));
 81:       /* Compute the necessary dot products */
 82:       PetscCall(VecMDot(lmvm->Y[i], i, lmvm->S, lsb->workscalar));
 83:       for (j = 0; j < i; ++j) {
 84:         sjtyi = lsb->workscalar[j];
 85:         PetscCall(VecDot(lmvm->Y[j], lsb->Q[i], &yjtqi));
 86:         /* Compute the pure DFP component of the inverse application*/
 87:         PetscCall(VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi) / lsb->ytq[j], PetscRealPart(sjtyi) / lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]));
 88:         /* Tack on the convexly scaled extras to the inverse application*/
 89:         if (lsb->psi[j] > 0.0) {
 90:           PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]));
 91:           PetscCall(VecDot(lsb->work, lmvm->Y[i], &wtyi));
 92:           PetscCall(VecAXPY(lsb->Q[i], lsb->psi[j] * lsb->ytq[j] * PetscRealPart(wtyi), lsb->work));
 93:         }
 94:       }
 95:       PetscCall(VecDot(lmvm->Y[i], lsb->Q[i], &ytq));
 96:       lsb->ytq[i] = PetscRealPart(ytq);
 97:       if (lsb->phi == 1.0) {
 98:         lsb->psi[i] = 0.0;
 99:       } else if (lsb->phi == 0.0) {
100:         lsb->psi[i] = 1.0;
101:       } else {
102:         numer       = (1.0 - lsb->phi) * lsb->yts[i] * lsb->yts[i];
103:         lsb->psi[i] = numer / (numer + (lsb->phi * lsb->ytq[i] * lsb->stp[i]));
104:       }
105:     }
106:     lsb->needQ = PETSC_FALSE;
107:   }

109:   /* Start the outer iterations for ((B^{-1}) * dX) */
110:   PetscCall(MatSymBrdnApplyJ0Inv(B, F, dX));
111:   /* Get all the dot products we need */
112:   PetscCall(VecMDot(F, lmvm->k + 1, lmvm->S, lsb->workscalar));
113:   for (i = 0; i <= lmvm->k; ++i) {
114:     stf = lsb->workscalar[i];
115:     PetscCall(VecDot(lmvm->Y[i], dX, &ytx));
116:     /* Compute the pure DFP component */
117:     PetscCall(VecAXPBYPCZ(dX, -PetscRealPart(ytx) / lsb->ytq[i], PetscRealPart(stf) / lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]));
118:     /* Tack on the convexly scaled extras */
119:     PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]));
120:     PetscCall(VecDot(lsb->work, F, &wtf));
121:     PetscCall(VecAXPY(dX, lsb->psi[i] * lsb->ytq[i] * PetscRealPart(wtf), lsb->work));
122:   }
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: /*------------------------------------------------------------*/

128: /*
129:   The forward-product below is the matrix-free implementation of
130:   Equation 16 in Dennis and Wolkowicz "Sizing and Least Change Secant
131:   Methods" (http://www.caam.rice.edu/caam/trs/90/TR90-05.pdf).

133:   P[i] = (B_i)*S[i] terms are computed ahead of time whenever
134:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
135:   repeated calls of MatMult inside KSP solvers without unnecessarily
136:   recomputing P[i] terms in expensive nested-loops.

138:   Z <- J0 * X

140:   for i=0,1,2,...,k
141:     # P[i] = (B_k) * S[i]

143:     rho = 1.0 / (Y[i]^T S[i])
144:     alpha = rho * (Y[i]^T F)
145:     zeta = 1.0 / (S[i]^T P[i])
146:     gamma = zeta * (S[i]^T dX)

148:     dX <- dX - (gamma * P[i]) + (alpha * S[i])
149:     W <- (rho * Y[i]) - (zeta * P[i])
150:     dX <- dX + (phi * (S[i]^T P[i]) * (W^T F) * W)
151:   end
152: */
153: static PetscErrorCode MatMult_LMVMSymBrdn(Mat B, Vec X, Vec Z)
154: {
155:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
156:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
157:   PetscInt     i, j;
158:   PetscScalar  sjtpi, yjtsi, wtsi, stz, ytx, wtx, stp;

160:   PetscFunctionBegin;
161:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
162:   if (lsb->phi == 0.0) {
163:     PetscCall(MatMult_LMVMBFGS(B, X, Z));
164:     PetscFunctionReturn(PETSC_SUCCESS);
165:   }
166:   if (lsb->phi == 1.0) {
167:     PetscCall(MatMult_LMVMDFP(B, X, Z));
168:     PetscFunctionReturn(PETSC_SUCCESS);
169:   }

171:   VecCheckSameSize(X, 2, Z, 3);
172:   VecCheckMatCompatible(B, X, 2, Z, 3);

174:   if (lsb->needP) {
175:     /* Start the loop for (P[k] = (B_k) * S[k]) */
176:     for (i = 0; i <= lmvm->k; ++i) {
177:       PetscCall(MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]));
178:       /* Compute the necessary dot products */
179:       PetscCall(VecMDot(lmvm->S[i], i, lmvm->Y, lsb->workscalar));
180:       for (j = 0; j < i; ++j) {
181:         yjtsi = lsb->workscalar[j];
182:         PetscCall(VecDot(lmvm->S[j], lsb->P[i], &sjtpi));
183:         /* Compute the pure BFGS component of the forward product */
184:         PetscCall(VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi) / lsb->stp[j], PetscRealPart(yjtsi) / lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]));
185:         /* Tack on the convexly scaled extras to the forward product */
186:         if (lsb->phi > 0.0) {
187:           PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[j], -1.0 / lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]));
188:           PetscCall(VecDot(lsb->work, lmvm->S[i], &wtsi));
189:           PetscCall(VecAXPY(lsb->P[i], lsb->phi * lsb->stp[j] * PetscRealPart(wtsi), lsb->work));
190:         }
191:       }
192:       PetscCall(VecDot(lmvm->S[i], lsb->P[i], &stp));
193:       lsb->stp[i] = PetscRealPart(stp);
194:     }
195:     lsb->needP = PETSC_FALSE;
196:   }

198:   /* Start the outer iterations for (B * X) */
199:   PetscCall(MatSymBrdnApplyJ0Fwd(B, X, Z));
200:   /* Get all the dot products we need */
201:   PetscCall(VecMDot(X, lmvm->k + 1, lmvm->Y, lsb->workscalar));
202:   for (i = 0; i <= lmvm->k; ++i) {
203:     ytx = lsb->workscalar[i];
204:     PetscCall(VecDot(lmvm->S[i], Z, &stz));
205:     /* Compute the pure BFGS component */
206:     PetscCall(VecAXPBYPCZ(Z, -PetscRealPart(stz) / lsb->stp[i], PetscRealPart(ytx) / lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]));
207:     /* Tack on the convexly scaled extras */
208:     PetscCall(VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]));
209:     PetscCall(VecDot(lsb->work, X, &wtx));
210:     PetscCall(VecAXPY(Z, lsb->phi * lsb->stp[i] * PetscRealPart(wtx), lsb->work));
211:   }
212:   PetscFunctionReturn(PETSC_SUCCESS);
213: }

215: /*------------------------------------------------------------*/

217: static PetscErrorCode MatUpdate_LMVMSymBrdn(Mat B, Vec X, Vec F)
218: {
219:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
220:   Mat_SymBrdn  *lsb  = (Mat_SymBrdn *)lmvm->ctx;
221:   Mat_LMVM     *dbase;
222:   Mat_DiagBrdn *dctx;
223:   PetscInt      old_k, i;
224:   PetscReal     curvtol, ststmp;
225:   PetscScalar   curvature, ytytmp;

227:   PetscFunctionBegin;
228:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
229:   if (lmvm->prev_set) {
230:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
231:     PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
232:     PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));

234:     /* Test if the updates can be accepted */
235:     PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ststmp));
236:     if (ststmp < lmvm->eps) curvtol = 0.0;
237:     else curvtol = lmvm->eps * ststmp;

239:     if (PetscRealPart(curvature) > curvtol) {
240:       /* Update is good, accept it */
241:       lsb->watchdog = 0;
242:       lsb->needP = lsb->needQ = PETSC_TRUE;
243:       old_k                   = lmvm->k;
244:       PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
245:       /* If we hit the memory limit, shift the yts, yty and sts arrays */
246:       if (old_k == lmvm->k) {
247:         for (i = 0; i <= lmvm->k - 1; ++i) {
248:           lsb->yts[i] = lsb->yts[i + 1];
249:           lsb->yty[i] = lsb->yty[i + 1];
250:           lsb->sts[i] = lsb->sts[i + 1];
251:         }
252:       }
253:       /* Update history of useful scalars */
254:       PetscCall(VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp));
255:       lsb->yts[lmvm->k] = PetscRealPart(curvature);
256:       lsb->yty[lmvm->k] = PetscRealPart(ytytmp);
257:       lsb->sts[lmvm->k] = ststmp;
258:       /* Compute the scalar scale if necessary */
259:       if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(MatSymBrdnComputeJ0Scalar(B));
260:     } else {
261:       /* Update is bad, skip it */
262:       ++lmvm->nrejects;
263:       ++lsb->watchdog;
264:     }
265:   } else {
266:     switch (lsb->scale_type) {
267:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
268:       dbase = (Mat_LMVM *)lsb->D->data;
269:       dctx  = (Mat_DiagBrdn *)dbase->ctx;
270:       PetscCall(VecSet(dctx->invD, lsb->delta));
271:       break;
272:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
273:       lsb->sigma = lsb->delta;
274:       break;
275:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
276:       lsb->sigma = 1.0;
277:       break;
278:     default:
279:       break;
280:     }
281:   }

283:   /* Update the scaling */
284:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMUpdate(lsb->D, X, F));

286:   if (lsb->watchdog > lsb->max_seq_rejects) {
287:     PetscCall(MatLMVMReset(B, PETSC_FALSE));
288:     if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMReset(lsb->D, PETSC_FALSE));
289:   }

291:   /* Save the solution and function to be used in the next update */
292:   PetscCall(VecCopy(X, lmvm->Xprev));
293:   PetscCall(VecCopy(F, lmvm->Fprev));
294:   lmvm->prev_set = PETSC_TRUE;
295:   PetscFunctionReturn(PETSC_SUCCESS);
296: }

298: /*------------------------------------------------------------*/

300: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
301: {
302:   Mat_LMVM    *bdata = (Mat_LMVM *)B->data;
303:   Mat_SymBrdn *blsb  = (Mat_SymBrdn *)bdata->ctx;
304:   Mat_LMVM    *mdata = (Mat_LMVM *)M->data;
305:   Mat_SymBrdn *mlsb  = (Mat_SymBrdn *)mdata->ctx;
306:   PetscInt     i;

308:   PetscFunctionBegin;
309:   mlsb->phi   = blsb->phi;
310:   mlsb->needP = blsb->needP;
311:   mlsb->needQ = blsb->needQ;
312:   for (i = 0; i <= bdata->k; ++i) {
313:     mlsb->stp[i] = blsb->stp[i];
314:     mlsb->ytq[i] = blsb->ytq[i];
315:     mlsb->yts[i] = blsb->yts[i];
316:     mlsb->psi[i] = blsb->psi[i];
317:     PetscCall(VecCopy(blsb->P[i], mlsb->P[i]));
318:     PetscCall(VecCopy(blsb->Q[i], mlsb->Q[i]));
319:   }
320:   mlsb->scale_type      = blsb->scale_type;
321:   mlsb->alpha           = blsb->alpha;
322:   mlsb->beta            = blsb->beta;
323:   mlsb->rho             = blsb->rho;
324:   mlsb->delta           = blsb->delta;
325:   mlsb->sigma_hist      = blsb->sigma_hist;
326:   mlsb->watchdog        = blsb->watchdog;
327:   mlsb->max_seq_rejects = blsb->max_seq_rejects;
328:   switch (blsb->scale_type) {
329:   case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
330:     mlsb->sigma = blsb->sigma;
331:     break;
332:   case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
333:     PetscCall(MatCopy(blsb->D, mlsb->D, SAME_NONZERO_PATTERN));
334:     break;
335:   case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
336:     mlsb->sigma = 1.0;
337:     break;
338:   default:
339:     break;
340:   }
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }

344: /*------------------------------------------------------------*/

346: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
347: {
348:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
349:   Mat_SymBrdn  *lsb  = (Mat_SymBrdn *)lmvm->ctx;
350:   Mat_LMVM     *dbase;
351:   Mat_DiagBrdn *dctx;

353:   PetscFunctionBegin;
354:   lsb->watchdog = 0;
355:   lsb->needP = lsb->needQ = PETSC_TRUE;
356:   if (lsb->allocated) {
357:     if (destructive) {
358:       PetscCall(VecDestroy(&lsb->work));
359:       PetscCall(PetscFree6(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts, lsb->workscalar));
360:       PetscCall(PetscFree(lsb->psi));
361:       PetscCall(VecDestroyVecs(lmvm->m, &lsb->P));
362:       PetscCall(VecDestroyVecs(lmvm->m, &lsb->Q));
363:       switch (lsb->scale_type) {
364:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
365:         PetscCall(MatLMVMReset(lsb->D, PETSC_TRUE));
366:         break;
367:       default:
368:         break;
369:       }
370:       lsb->allocated = PETSC_FALSE;
371:     } else {
372:       PetscCall(PetscMemzero(lsb->psi, lmvm->m));
373:       switch (lsb->scale_type) {
374:       case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
375:         lsb->sigma = lsb->delta;
376:         break;
377:       case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
378:         PetscCall(MatLMVMReset(lsb->D, PETSC_FALSE));
379:         dbase = (Mat_LMVM *)lsb->D->data;
380:         dctx  = (Mat_DiagBrdn *)dbase->ctx;
381:         PetscCall(VecSet(dctx->invD, lsb->delta));
382:         break;
383:       case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
384:         lsb->sigma = 1.0;
385:         break;
386:       default:
387:         break;
388:       }
389:     }
390:   }
391:   PetscCall(MatReset_LMVM(B, destructive));
392:   PetscFunctionReturn(PETSC_SUCCESS);
393: }

395: /*------------------------------------------------------------*/

397: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
398: {
399:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
400:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

402:   PetscFunctionBegin;
403:   PetscCall(MatAllocate_LMVM(B, X, F));
404:   if (!lsb->allocated) {
405:     PetscCall(VecDuplicate(X, &lsb->work));
406:     if (lmvm->m > 0) {
407:       PetscCall(PetscMalloc6(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts, lmvm->m, &lsb->workscalar));
408:       PetscCall(PetscCalloc1(lmvm->m, &lsb->psi));
409:       PetscCall(VecDuplicateVecs(X, lmvm->m, &lsb->P));
410:       PetscCall(VecDuplicateVecs(X, lmvm->m, &lsb->Q));
411:     }
412:     switch (lsb->scale_type) {
413:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
414:       PetscCall(MatLMVMAllocate(lsb->D, X, F));
415:       break;
416:     default:
417:       break;
418:     }
419:     lsb->allocated = PETSC_TRUE;
420:   }
421:   PetscFunctionReturn(PETSC_SUCCESS);
422: }

424: /*------------------------------------------------------------*/

426: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
427: {
428:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
429:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

431:   PetscFunctionBegin;
432:   if (lsb->allocated) {
433:     PetscCall(VecDestroy(&lsb->work));
434:     PetscCall(PetscFree6(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts, lsb->workscalar));
435:     PetscCall(PetscFree(lsb->psi));
436:     PetscCall(VecDestroyVecs(lmvm->m, &lsb->P));
437:     PetscCall(VecDestroyVecs(lmvm->m, &lsb->Q));
438:     lsb->allocated = PETSC_FALSE;
439:   }
440:   PetscCall(MatDestroy(&lsb->D));
441:   PetscCall(PetscFree(lmvm->ctx));
442:   PetscCall(MatDestroy_LMVM(B));
443:   PetscFunctionReturn(PETSC_SUCCESS);
444: }

446: /*------------------------------------------------------------*/

448: static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
449: {
450:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
451:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
452:   PetscInt     n, N;

454:   PetscFunctionBegin;
455:   PetscCall(MatSetUp_LMVM(B));
456:   if (!lsb->allocated) {
457:     PetscCall(VecDuplicate(lmvm->Xprev, &lsb->work));
458:     if (lmvm->m > 0) {
459:       PetscCall(PetscMalloc6(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts, lmvm->m, &lsb->workscalar));
460:       PetscCall(PetscCalloc1(lmvm->m, &lsb->psi));
461:       PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->P));
462:       PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->Q));
463:     }
464:     switch (lsb->scale_type) {
465:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
466:       PetscCall(MatGetLocalSize(B, &n, &n));
467:       PetscCall(MatGetSize(B, &N, &N));
468:       PetscCall(MatSetSizes(lsb->D, n, n, N, N));
469:       PetscCall(MatSetUp(lsb->D));
470:       break;
471:     default:
472:       break;
473:     }
474:     lsb->allocated = PETSC_TRUE;
475:   }
476:   PetscFunctionReturn(PETSC_SUCCESS);
477: }

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

481: PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
482: {
483:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
484:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
485:   PetscBool    isascii;

487:   PetscFunctionBegin;
488:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
489:   if (isascii) {
490:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale type: %s\n", MatLMVMSymBroydenScaleTypes[lsb->scale_type]));
491:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", lsb->sigma_hist));
492:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)lsb->alpha, (double)lsb->beta, (double)lsb->rho));
493:     PetscCall(PetscViewerASCIIPrintf(pv, "Convex factors: phi=%g, theta=%g\n", (double)lsb->phi, (double)lsb->theta));
494:   }
495:   PetscCall(MatView_LMVM(B, pv));
496:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatView(lsb->D, pv));
497:   PetscFunctionReturn(PETSC_SUCCESS);
498: }

500: /*------------------------------------------------------------*/

502: PetscErrorCode MatSetFromOptions_LMVMSymBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
503: {
504:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
505:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

507:   PetscFunctionBegin;
508:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
509:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
510:   PetscCall(PetscOptionsReal("-mat_lmvm_phi", "(developer) convex ratio between BFGS and DFP components of the update", "", lsb->phi, &lsb->phi, NULL));
511:   PetscCheck(!(lsb->phi < 0.0) && !(lsb->phi > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the update formula cannot be outside the range of [0, 1]");
512:   PetscCall(MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject));
513:   PetscOptionsHeadEnd();
514:   PetscFunctionReturn(PETSC_SUCCESS);
515: }

517: PetscErrorCode MatSetFromOptions_LMVMSymBrdn_Private(Mat B, PetscOptionItems *PetscOptionsObject)
518: {
519:   Mat_LMVM                  *lmvm = (Mat_LMVM *)B->data;
520:   Mat_SymBrdn               *lsb  = (Mat_SymBrdn *)lmvm->ctx;
521:   Mat_LMVM                  *dbase;
522:   Mat_DiagBrdn              *dctx;
523:   MatLMVMSymBroydenScaleType stype = lsb->scale_type;
524:   PetscBool                  flg;

526:   PetscFunctionBegin;
527:   PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", lsb->beta, &lsb->beta, NULL));
528:   PetscCall(PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", lsb->theta, &lsb->theta, NULL));
529:   PetscCheck(!(lsb->theta < 0.0) && !(lsb->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]");
530:   PetscCall(PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", lsb->rho, &lsb->rho, NULL));
531:   PetscCheck(!(lsb->rho < 0.0) && !(lsb->rho > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "update limiter in the J0 scaling cannot be outside the range of [0, 1]");
532:   PetscCall(PetscOptionsReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", lsb->alpha, &lsb->alpha, NULL));
533:   PetscCheck(!(lsb->alpha < 0.0) && !(lsb->alpha > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
534:   PetscCall(PetscOptionsBoundedInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", lsb->sigma_hist, &lsb->sigma_hist, NULL, 1));
535:   PetscCall(PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBrdnScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg));
536:   if (flg) PetscCall(MatLMVMSymBroydenSetScaleType(B, stype));
537:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
538:     const char *prefix;

540:     PetscCall(MatGetOptionsPrefix(B, &prefix));
541:     PetscCall(MatSetOptionsPrefix(lsb->D, prefix));
542:     PetscCall(MatAppendOptionsPrefix(lsb->D, "J0_"));
543:     PetscCall(MatSetFromOptions(lsb->D));
544:     dbase            = (Mat_LMVM *)lsb->D->data;
545:     dctx             = (Mat_DiagBrdn *)dbase->ctx;
546:     dctx->delta_min  = lsb->delta_min;
547:     dctx->delta_max  = lsb->delta_max;
548:     dctx->theta      = lsb->theta;
549:     dctx->rho        = lsb->rho;
550:     dctx->alpha      = lsb->alpha;
551:     dctx->beta       = lsb->beta;
552:     dctx->sigma_hist = lsb->sigma_hist;
553:   }
554:   PetscFunctionReturn(PETSC_SUCCESS);
555: }

557: /*------------------------------------------------------------*/

559: PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
560: {
561:   Mat_LMVM    *lmvm;
562:   Mat_SymBrdn *lsb;

564:   PetscFunctionBegin;
565:   PetscCall(MatCreate_LMVM(B));
566:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN));
567:   PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
568:   PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
569:   B->ops->view           = MatView_LMVMSymBrdn;
570:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
571:   B->ops->setup          = MatSetUp_LMVMSymBrdn;
572:   B->ops->destroy        = MatDestroy_LMVMSymBrdn;
573:   B->ops->solve          = MatSolve_LMVMSymBrdn;

575:   lmvm                = (Mat_LMVM *)B->data;
576:   lmvm->square        = PETSC_TRUE;
577:   lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
578:   lmvm->ops->reset    = MatReset_LMVMSymBrdn;
579:   lmvm->ops->update   = MatUpdate_LMVMSymBrdn;
580:   lmvm->ops->mult     = MatMult_LMVMSymBrdn;
581:   lmvm->ops->copy     = MatCopy_LMVMSymBrdn;

583:   PetscCall(PetscNew(&lsb));
584:   lmvm->ctx      = (void *)lsb;
585:   lsb->allocated = PETSC_FALSE;
586:   lsb->needP = lsb->needQ = PETSC_TRUE;
587:   lsb->phi                = 0.125;
588:   lsb->theta              = 0.125;
589:   lsb->alpha              = 1.0;
590:   lsb->rho                = 1.0;
591:   lsb->beta               = 0.5;
592:   lsb->sigma              = 1.0;
593:   lsb->delta              = 1.0;
594:   lsb->delta_min          = 1e-7;
595:   lsb->delta_max          = 100.0;
596:   lsb->sigma_hist         = 1;
597:   lsb->scale_type         = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
598:   lsb->watchdog           = 0;
599:   lsb->max_seq_rejects    = lmvm->m / 2;

601:   PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &lsb->D));
602:   PetscCall(MatSetType(lsb->D, MATLMVMDIAGBROYDEN));
603:   PetscFunctionReturn(PETSC_SUCCESS);
604: }

606: /*------------------------------------------------------------*/

608: /*@
609:    MatLMVMSymBroydenSetDelta - Sets the starting value for the diagonal scaling vector computed
610:    in the SymBrdn approximations (also works for BFGS and DFP).

612:    Input Parameters:
613: +  B - LMVM matrix
614: -  delta - initial value for diagonal scaling

616:    Level: intermediate

618: @*/

620: PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
621: {
622:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
623:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
624:   PetscBool    is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn;

626:   PetscFunctionBegin;
627:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs));
628:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp));
629:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn));
630:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn));
631:   PetscCheck(is_bfgs || is_dfp || is_symbrdn || is_symbadbrdn, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "diagonal scaling is only available for DFP, BFGS and SymBrdn matrices");
632:   lsb->delta = PetscAbsReal(PetscRealPart(delta));
633:   lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
634:   lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
635:   PetscFunctionReturn(PETSC_SUCCESS);
636: }

638: /*------------------------------------------------------------*/

640: /*@
641:     MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.

643:     Input Parameters:
644: +   snes - the iterative context
645: -   rtype - restart type

647:     Options Database Key:
648: .   -mat_lmvm_scale_type <none,scalar,diagonal> - set the scaling type

650:     Level: intermediate

652:     MatLMVMSymBrdnScaleTypes:
653: +   MAT_LMVM_SYMBROYDEN_SCALE_NONE - initial Hessian is the identity matrix
654: .   MAT_LMVM_SYMBROYDEN_SCALE_SCALAR - use the Shanno scalar as the initial Hessian
655: -   MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL - use a diagonalized BFGS update as the initial Hessian

657: .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`, `MatCreateLMVMSymBroyden()`
658: @*/
659: PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
660: {
661:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
662:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

664:   PetscFunctionBegin;
666:   lsb->scale_type = stype;
667:   PetscFunctionReturn(PETSC_SUCCESS);
668: }

670: /*------------------------------------------------------------*/

672: /*@
673:    MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
674:    for approximating Jacobians. L-SymBrdn is a convex combination of L-DFP and
675:    L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor
676:    phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed
677:    to be symmetric positive-definite.

679:    To use the L-SymBrdn matrix with other vector types, the matrix must be
680:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
681:    This ensures that the internal storage and work vectors are duplicated from the
682:    correct type of vector.

684:    Collective

686:    Input Parameters:
687: +  comm - MPI communicator, set to PETSC_COMM_SELF
688: .  n - number of local rows for storage vectors
689: -  N - global size of the storage vectors

691:    Output Parameter:
692: .  B - the matrix

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

697:    Options Database Keys:
698: +   -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
699: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
700: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
701: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
702: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
703: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
704: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

706:    Level: intermediate

708: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
709:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`
710: @*/
711: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
712: {
713:   PetscFunctionBegin;
714:   PetscCall(MatCreate(comm, B));
715:   PetscCall(MatSetSizes(*B, n, n, N, N));
716:   PetscCall(MatSetType(*B, MATLMVMSYMBROYDEN));
717:   PetscCall(MatSetUp(*B));
718:   PetscFunctionReturn(PETSC_SUCCESS);
719: }

721: /*------------------------------------------------------------*/

723: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
724: {
725:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
726:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

728:   PetscFunctionBegin;
729:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
730:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
731:     PetscCall(MatLMVMApplyJ0Fwd(B, X, Z));
732:   } else {
733:     switch (lsb->scale_type) {
734:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
735:       PetscCall(VecAXPBY(Z, 1.0 / lsb->sigma, 0.0, X));
736:       break;
737:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
738:       PetscCall(MatMult(lsb->D, X, Z));
739:       break;
740:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
741:     default:
742:       PetscCall(VecCopy(X, Z));
743:       break;
744:     }
745:   }
746:   PetscFunctionReturn(PETSC_SUCCESS);
747: }

749: /*------------------------------------------------------------*/

751: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
752: {
753:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
754:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

756:   PetscFunctionBegin;
757:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
758:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
759:     PetscCall(MatLMVMApplyJ0Inv(B, F, dX));
760:   } else {
761:     switch (lsb->scale_type) {
762:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
763:       PetscCall(VecAXPBY(dX, lsb->sigma, 0.0, F));
764:       break;
765:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
766:       PetscCall(MatSolve(lsb->D, F, dX));
767:       break;
768:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
769:     default:
770:       PetscCall(VecCopy(F, dX));
771:       break;
772:     }
773:   }
774:   PetscFunctionReturn(PETSC_SUCCESS);
775: }

777: /*------------------------------------------------------------*/

779: PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
780: {
781:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
782:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
783:   PetscInt     i, start;
784:   PetscReal    a, b, c, sig1, sig2, signew;

786:   PetscFunctionBegin;
787:   if (lsb->sigma_hist == 0) {
788:     signew = 1.0;
789:   } else {
790:     start  = PetscMax(0, lmvm->k - lsb->sigma_hist + 1);
791:     signew = 0.0;
792:     if (lsb->alpha == 1.0) {
793:       for (i = start; i <= lmvm->k; ++i) signew += lsb->yts[i] / lsb->yty[i];
794:     } else if (lsb->alpha == 0.5) {
795:       for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yty[i];
796:       signew = PetscSqrtReal(signew);
797:     } else if (lsb->alpha == 0.0) {
798:       for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yts[i];
799:     } else {
800:       /* compute coefficients of the quadratic */
801:       a = b = c = 0.0;
802:       for (i = start; i <= lmvm->k; ++i) {
803:         a += lsb->yty[i];
804:         b += lsb->yts[i];
805:         c += lsb->sts[i];
806:       }
807:       a *= lsb->alpha;
808:       b *= -(2.0 * lsb->alpha - 1.0);
809:       c *= lsb->alpha - 1.0;
810:       /* use quadratic formula to find roots */
811:       sig1 = (-b + PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
812:       sig2 = (-b - PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
813:       /* accept the positive root as the scalar */
814:       if (sig1 > 0.0) {
815:         signew = sig1;
816:       } else if (sig2 > 0.0) {
817:         signew = sig2;
818:       } else {
819:         SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
820:       }
821:     }
822:   }
823:   lsb->sigma = lsb->rho * signew + (1.0 - lsb->rho) * lsb->sigma;
824:   PetscFunctionReturn(PETSC_SUCCESS);
825: }