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

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

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

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

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

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

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

 49:   VecCheckSameSize(F, 2, dX, 3);
 50:   VecCheckMatCompatible(B, dX, 3, F, 2);

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

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

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

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

134:   Z <- J0 * X

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

139:     rho = 1.0 / (Y[i]^T S[i])
140:     alpha = rho * (Y[i]^T F)
141:     zeta = 1.0 / (S[i]^T P[i])
142:     gamma = zeta * (S[i]^T dX)

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

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

167:   VecCheckSameSize(X, 2, Z, 3);
168:   VecCheckMatCompatible(B, X, 2, Z, 3);

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

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

211: static PetscErrorCode MatUpdate_LMVMSymBrdn(Mat B, Vec X, Vec F)
212: {
213:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
214:   Mat_SymBrdn  *lsb  = (Mat_SymBrdn *)lmvm->ctx;
215:   Mat_LMVM     *dbase;
216:   Mat_DiagBrdn *dctx;
217:   PetscInt      old_k, i;
218:   PetscReal     curvtol, ststmp;
219:   PetscScalar   curvature, ytytmp;

221:   PetscFunctionBegin;
222:   if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
223:   if (lmvm->prev_set) {
224:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
225:     PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
226:     PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));

228:     /* Test if the updates can be accepted */
229:     PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ststmp));
230:     if (ststmp < lmvm->eps) curvtol = 0.0;
231:     else curvtol = lmvm->eps * ststmp;

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

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

280:   if (lsb->watchdog > lsb->max_seq_rejects) {
281:     PetscCall(MatLMVMReset(B, PETSC_FALSE));
282:     if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMReset(lsb->D, PETSC_FALSE));
283:   }

285:   /* Save the solution and function to be used in the next update */
286:   PetscCall(VecCopy(X, lmvm->Xprev));
287:   PetscCall(VecCopy(F, lmvm->Fprev));
288:   lmvm->prev_set = PETSC_TRUE;
289:   PetscFunctionReturn(PETSC_SUCCESS);
290: }

292: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
293: {
294:   Mat_LMVM    *bdata = (Mat_LMVM *)B->data;
295:   Mat_SymBrdn *blsb  = (Mat_SymBrdn *)bdata->ctx;
296:   Mat_LMVM    *mdata = (Mat_LMVM *)M->data;
297:   Mat_SymBrdn *mlsb  = (Mat_SymBrdn *)mdata->ctx;
298:   PetscInt     i;

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

336: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
337: {
338:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
339:   Mat_SymBrdn  *lsb  = (Mat_SymBrdn *)lmvm->ctx;
340:   Mat_LMVM     *dbase;
341:   Mat_DiagBrdn *dctx;

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

385: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
386: {
387:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
388:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

390:   PetscFunctionBegin;
391:   PetscCall(MatAllocate_LMVM(B, X, F));
392:   if (!lsb->allocated) {
393:     PetscCall(VecDuplicate(X, &lsb->work));
394:     if (lmvm->m > 0) {
395:       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));
396:       PetscCall(PetscCalloc1(lmvm->m, &lsb->psi));
397:       PetscCall(VecDuplicateVecs(X, lmvm->m, &lsb->P));
398:       PetscCall(VecDuplicateVecs(X, lmvm->m, &lsb->Q));
399:     }
400:     switch (lsb->scale_type) {
401:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
402:       PetscCall(MatLMVMAllocate(lsb->D, X, F));
403:       break;
404:     default:
405:       break;
406:     }
407:     lsb->allocated = PETSC_TRUE;
408:   }
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
413: {
414:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
415:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

417:   PetscFunctionBegin;
418:   if (lsb->allocated) {
419:     PetscCall(VecDestroy(&lsb->work));
420:     PetscCall(PetscFree6(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts, lsb->workscalar));
421:     PetscCall(PetscFree(lsb->psi));
422:     PetscCall(VecDestroyVecs(lmvm->m, &lsb->P));
423:     PetscCall(VecDestroyVecs(lmvm->m, &lsb->Q));
424:     lsb->allocated = PETSC_FALSE;
425:   }
426:   PetscCall(MatDestroy(&lsb->D));
427:   PetscCall(PetscFree(lmvm->ctx));
428:   PetscCall(MatDestroy_LMVM(B));
429:   PetscFunctionReturn(PETSC_SUCCESS);
430: }

432: static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
433: {
434:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
435:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
436:   PetscInt     n, N;

438:   PetscFunctionBegin;
439:   PetscCall(MatSetUp_LMVM(B));
440:   if (!lsb->allocated) {
441:     PetscCall(VecDuplicate(lmvm->Xprev, &lsb->work));
442:     if (lmvm->m > 0) {
443:       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));
444:       PetscCall(PetscCalloc1(lmvm->m, &lsb->psi));
445:       PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->P));
446:       PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->Q));
447:     }
448:     switch (lsb->scale_type) {
449:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
450:       PetscCall(MatGetLocalSize(B, &n, &n));
451:       PetscCall(MatGetSize(B, &N, &N));
452:       PetscCall(MatSetSizes(lsb->D, n, n, N, N));
453:       PetscCall(MatSetUp(lsb->D));
454:       break;
455:     default:
456:       break;
457:     }
458:     lsb->allocated = PETSC_TRUE;
459:   }
460:   PetscFunctionReturn(PETSC_SUCCESS);
461: }

463: PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
464: {
465:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
466:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
467:   PetscBool    isascii;

469:   PetscFunctionBegin;
470:   PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
471:   if (isascii) {
472:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale type: %s\n", MatLMVMSymBroydenScaleTypes[lsb->scale_type]));
473:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale history: %" PetscInt_FMT "\n", lsb->sigma_hist));
474:     PetscCall(PetscViewerASCIIPrintf(pv, "Scale params: alpha=%g, beta=%g, rho=%g\n", (double)lsb->alpha, (double)lsb->beta, (double)lsb->rho));
475:     PetscCall(PetscViewerASCIIPrintf(pv, "Convex factors: phi=%g, theta=%g\n", (double)lsb->phi, (double)lsb->theta));
476:   }
477:   PetscCall(MatView_LMVM(B, pv));
478:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatView(lsb->D, pv));
479:   PetscFunctionReturn(PETSC_SUCCESS);
480: }

482: PetscErrorCode MatSetFromOptions_LMVMSymBrdn(Mat B, PetscOptionItems *PetscOptionsObject)
483: {
484:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
485:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

487:   PetscFunctionBegin;
488:   PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
489:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
490:   PetscCall(PetscOptionsReal("-mat_lmvm_phi", "(developer) convex ratio between BFGS and DFP components of the update", "", lsb->phi, &lsb->phi, NULL));
491:   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]");
492:   PetscCall(MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject));
493:   PetscOptionsHeadEnd();
494:   PetscFunctionReturn(PETSC_SUCCESS);
495: }

497: PetscErrorCode MatSetFromOptions_LMVMSymBrdn_Private(Mat B, PetscOptionItems *PetscOptionsObject)
498: {
499:   Mat_LMVM                  *lmvm = (Mat_LMVM *)B->data;
500:   Mat_SymBrdn               *lsb  = (Mat_SymBrdn *)lmvm->ctx;
501:   Mat_LMVM                  *dbase;
502:   Mat_DiagBrdn              *dctx;
503:   MatLMVMSymBroydenScaleType stype = lsb->scale_type;
504:   PetscBool                  flg;

506:   PetscFunctionBegin;
507:   PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", lsb->beta, &lsb->beta, NULL));
508:   PetscCall(PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", lsb->theta, &lsb->theta, NULL));
509:   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]");
510:   PetscCall(PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", lsb->rho, &lsb->rho, NULL));
511:   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]");
512:   PetscCall(PetscOptionsReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", lsb->alpha, &lsb->alpha, NULL));
513:   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]");
514:   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));
515:   PetscCall(PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBrdnScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg));
516:   if (flg) PetscCall(MatLMVMSymBroydenSetScaleType(B, stype));
517:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
518:     const char *prefix;

520:     PetscCall(MatGetOptionsPrefix(B, &prefix));
521:     PetscCall(MatSetOptionsPrefix(lsb->D, prefix));
522:     PetscCall(MatAppendOptionsPrefix(lsb->D, "J0_"));
523:     PetscCall(MatSetFromOptions(lsb->D));
524:     dbase            = (Mat_LMVM *)lsb->D->data;
525:     dctx             = (Mat_DiagBrdn *)dbase->ctx;
526:     dctx->delta_min  = lsb->delta_min;
527:     dctx->delta_max  = lsb->delta_max;
528:     dctx->theta      = lsb->theta;
529:     dctx->rho        = lsb->rho;
530:     dctx->alpha      = lsb->alpha;
531:     dctx->beta       = lsb->beta;
532:     dctx->sigma_hist = lsb->sigma_hist;
533:   }
534:   PetscFunctionReturn(PETSC_SUCCESS);
535: }

537: PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
538: {
539:   Mat_LMVM    *lmvm;
540:   Mat_SymBrdn *lsb;

542:   PetscFunctionBegin;
543:   PetscCall(MatCreate_LMVM(B));
544:   PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN));
545:   PetscCall(MatSetOption(B, MAT_SPD, PETSC_TRUE));
546:   PetscCall(MatSetOption(B, MAT_SPD_ETERNAL, PETSC_TRUE));
547:   B->ops->view           = MatView_LMVMSymBrdn;
548:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
549:   B->ops->setup          = MatSetUp_LMVMSymBrdn;
550:   B->ops->destroy        = MatDestroy_LMVMSymBrdn;
551:   B->ops->solve          = MatSolve_LMVMSymBrdn;

553:   lmvm                = (Mat_LMVM *)B->data;
554:   lmvm->square        = PETSC_TRUE;
555:   lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
556:   lmvm->ops->reset    = MatReset_LMVMSymBrdn;
557:   lmvm->ops->update   = MatUpdate_LMVMSymBrdn;
558:   lmvm->ops->mult     = MatMult_LMVMSymBrdn;
559:   lmvm->ops->copy     = MatCopy_LMVMSymBrdn;

561:   PetscCall(PetscNew(&lsb));
562:   lmvm->ctx      = (void *)lsb;
563:   lsb->allocated = PETSC_FALSE;
564:   lsb->needP = lsb->needQ = PETSC_TRUE;
565:   lsb->phi                = 0.125;
566:   lsb->theta              = 0.125;
567:   lsb->alpha              = 1.0;
568:   lsb->rho                = 1.0;
569:   lsb->beta               = 0.5;
570:   lsb->sigma              = 1.0;
571:   lsb->delta              = 1.0;
572:   lsb->delta_min          = 1e-7;
573:   lsb->delta_max          = 100.0;
574:   lsb->sigma_hist         = 1;
575:   lsb->scale_type         = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
576:   lsb->watchdog           = 0;
577:   lsb->max_seq_rejects    = lmvm->m / 2;

579:   PetscCall(MatCreate(PetscObjectComm((PetscObject)B), &lsb->D));
580:   PetscCall(MatSetType(lsb->D, MATLMVMDIAGBROYDEN));
581:   PetscFunctionReturn(PETSC_SUCCESS);
582: }

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

588:   Input Parameters:
589: + B     - `MATLMVM` matrix
590: - delta - initial value for diagonal scaling

592:   Level: intermediate

594: .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`
595: @*/
596: PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
597: {
598:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
599:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
600:   PetscBool    is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn;

602:   PetscFunctionBegin;
603:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs));
604:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp));
605:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn));
606:   PetscCall(PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn));
607:   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");
608:   lsb->delta = PetscAbsReal(PetscRealPart(delta));
609:   lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
610:   lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
611:   PetscFunctionReturn(PETSC_SUCCESS);
612: }

614: /*@
615:   MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.

617:   Input Parameters:
618: + B     - the `MATLMVM` matrix
619: - stype - scale type, see `MatLMVMSymBroydenScaleType`

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

624:   Level: intermediate

626:   MatLMVMSymBrdnScaleTypes\:
627: +   `MAT_LMVM_SYMBROYDEN_SCALE_NONE` - initial Hessian is the identity matrix
628: .   `MAT_LMVM_SYMBROYDEN_SCALE_SCALAR` - use the Shanno scalar as the initial Hessian
629: -   `MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL` - use a diagonalized BFGS update as the initial Hessian

631: .seealso: [](ch_ksp), `MATLMVMSYMBROYDEN`, `MatCreateLMVMSymBroyden()`, `MatLMVMSymBroydenScaleType`
632: @*/
633: PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
634: {
635:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
636:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

638:   PetscFunctionBegin;
640:   lsb->scale_type = stype;
641:   PetscFunctionReturn(PETSC_SUCCESS);
642: }

644: /*@
645:   MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
646:   for approximating Jacobians.

648:   Collective

650:   Input Parameters:
651: + comm - MPI communicator, set to `PETSC_COMM_SELF`
652: . n    - number of local rows for storage vectors
653: - N    - global size of the storage vectors

655:   Output Parameter:
656: . B - the matrix

658:   Options Database Keys:
659: + -mat_lmvm_phi        - (developer) convex ratio between BFGS and DFP components of the update
660: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
661: . -mat_lmvm_theta      - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
662: . -mat_lmvm_rho        - (developer) update limiter for the J0 scaling
663: . -mat_lmvm_alpha      - (developer) coefficient factor for the quadratic subproblem in J0 scaling
664: . -mat_lmvm_beta       - (developer) exponential factor for the diagonal J0 scaling
665: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

667:   Level: intermediate

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

673:   L-SymBrdn is a convex combination of L-DFP and
674:   L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor
675:   phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed
676:   to be symmetric positive-definite.

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

683: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMSYMBROYDEN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
684:           `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`
685: @*/
686: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
687: {
688:   PetscFunctionBegin;
689:   PetscCall(MatCreate(comm, B));
690:   PetscCall(MatSetSizes(*B, n, n, N, N));
691:   PetscCall(MatSetType(*B, MATLMVMSYMBROYDEN));
692:   PetscCall(MatSetUp(*B));
693:   PetscFunctionReturn(PETSC_SUCCESS);
694: }

696: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
697: {
698:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
699:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

701:   PetscFunctionBegin;
702:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
703:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
704:     PetscCall(MatLMVMApplyJ0Fwd(B, X, Z));
705:   } else {
706:     switch (lsb->scale_type) {
707:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
708:       PetscCall(VecAXPBY(Z, 1.0 / lsb->sigma, 0.0, X));
709:       break;
710:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
711:       PetscCall(MatMult(lsb->D, X, Z));
712:       break;
713:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
714:     default:
715:       PetscCall(VecCopy(X, Z));
716:       break;
717:     }
718:   }
719:   PetscFunctionReturn(PETSC_SUCCESS);
720: }

722: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
723: {
724:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
725:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

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

748: PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
749: {
750:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
751:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;
752:   PetscInt     i, start;
753:   PetscReal    a, b, c, sig1, sig2, signew;

755:   PetscFunctionBegin;
756:   if (lsb->sigma_hist == 0) {
757:     signew = 1.0;
758:   } else {
759:     start  = PetscMax(0, lmvm->k - lsb->sigma_hist + 1);
760:     signew = 0.0;
761:     if (lsb->alpha == 1.0) {
762:       for (i = start; i <= lmvm->k; ++i) signew += lsb->yts[i] / lsb->yty[i];
763:     } else if (lsb->alpha == 0.5) {
764:       for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yty[i];
765:       signew = PetscSqrtReal(signew);
766:     } else if (lsb->alpha == 0.0) {
767:       for (i = start; i <= lmvm->k; ++i) signew += lsb->sts[i] / lsb->yts[i];
768:     } else {
769:       /* compute coefficients of the quadratic */
770:       a = b = c = 0.0;
771:       for (i = start; i <= lmvm->k; ++i) {
772:         a += lsb->yty[i];
773:         b += lsb->yts[i];
774:         c += lsb->sts[i];
775:       }
776:       a *= lsb->alpha;
777:       b *= -(2.0 * lsb->alpha - 1.0);
778:       c *= lsb->alpha - 1.0;
779:       /* use quadratic formula to find roots */
780:       sig1 = (-b + PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
781:       sig2 = (-b - PetscSqrtReal(b * b - 4.0 * a * c)) / (2.0 * a);
782:       /* accept the positive root as the scalar */
783:       if (sig1 > 0.0) {
784:         signew = sig1;
785:       } else if (sig2 > 0.0) {
786:         signew = sig2;
787:       } else {
788:         SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
789:       }
790:     }
791:   }
792:   lsb->sigma = lsb->rho * signew + (1.0 - lsb->rho) * lsb->sigma;
793:   PetscFunctionReturn(PETSC_SUCCESS);
794: }