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:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
 41:   if (lsb->phi == 0.0) {
 42:     MatSolve_LMVMBFGS(B, F, dX);
 43:     return 0;
 44:   }
 45:   if (lsb->phi == 1.0) {
 46:     MatSolve_LMVMDFP(B, F, dX);
 47:     return 0;
 48:   }

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

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

110:   /* Start the outer iterations for ((B^{-1}) * dX) */
111:   MatSymBrdnApplyJ0Inv(B, F, dX);
112:   for (i = 0; i <= lmvm->k; ++i) {
113:     /* Compute the necessary dot products -- store yTs and yTp for inner iterations later */
114:     VecDotBegin(lmvm->Y[i], dX, &ytx);
115:     VecDotBegin(lmvm->S[i], F, &stf);
116:     VecDotEnd(lmvm->Y[i], dX, &ytx);
117:     VecDotEnd(lmvm->S[i], F, &stf);
118:     /* Compute the pure DFP component */
119:     VecAXPBYPCZ(dX, -PetscRealPart(ytx) / lsb->ytq[i], PetscRealPart(stf) / lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]);
120:     /* Tack on the convexly scaled extras */
121:     VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]);
122:     VecDot(lsb->work, F, &wtf);
123:     VecAXPY(dX, lsb->psi[i] * lsb->ytq[i] * PetscRealPart(wtf), lsb->work);
124:   }

126:   return 0;
127: }

129: /*------------------------------------------------------------*/

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

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

141:   Z <- J0 * X

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

146:     rho = 1.0 / (Y[i]^T S[i])
147:     alpha = rho * (Y[i]^T F)
148:     zeta = 1.0 / (S[i]^T P[i])
149:     gamma = zeta * (S[i]^T dX)

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

163:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
164:   if (lsb->phi == 0.0) {
165:     MatMult_LMVMBFGS(B, X, Z);
166:     return 0;
167:   }
168:   if (lsb->phi == 1.0) {
169:     MatMult_LMVMDFP(B, X, Z);
170:     return 0;
171:   }

173:   VecCheckSameSize(X, 2, Z, 3);
174:   VecCheckMatCompatible(B, X, 2, Z, 3);

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

201:   /* Start the outer iterations for (B * X) */
202:   MatSymBrdnApplyJ0Fwd(B, X, Z);
203:   for (i = 0; i <= lmvm->k; ++i) {
204:     /* Compute the necessary dot products */
205:     VecDotBegin(lmvm->S[i], Z, &stz);
206:     VecDotBegin(lmvm->Y[i], X, &ytx);
207:     VecDotEnd(lmvm->S[i], Z, &stz);
208:     VecDotEnd(lmvm->Y[i], X, &ytx);
209:     /* Compute the pure BFGS component */
210:     VecAXPBYPCZ(Z, -PetscRealPart(stz) / lsb->stp[i], PetscRealPart(ytx) / lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]);
211:     /* Tack on the convexly scaled extras */
212:     VecAXPBYPCZ(lsb->work, 1.0 / lsb->yts[i], -1.0 / lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]);
213:     VecDot(lsb->work, X, &wtx);
214:     VecAXPY(Z, lsb->phi * lsb->stp[i] * PetscRealPart(wtx), lsb->work);
215:   }
216:   return 0;
217: }

219: /*------------------------------------------------------------*/

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

231:   if (!lmvm->m) return 0;
232:   if (lmvm->prev_set) {
233:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
234:     VecAYPX(lmvm->Xprev, -1.0, X);
235:     VecAYPX(lmvm->Fprev, -1.0, F);
236:     /* Test if the updates can be accepted */
237:     VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
238:     VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
239:     VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
240:     VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
241:     if (PetscRealPart(ststmp) < lmvm->eps) {
242:       curvtol = 0.0;
243:     } else {
244:       curvtol = lmvm->eps * PetscRealPart(ststmp);
245:     }
246:     if (PetscRealPart(curvature) > curvtol) {
247:       /* Update is good, accept it */
248:       lsb->watchdog = 0;
249:       lsb->needP = lsb->needQ = PETSC_TRUE;
250:       old_k                   = lmvm->k;
251:       MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
252:       /* If we hit the memory limit, shift the yts, yty and sts arrays */
253:       if (old_k == lmvm->k) {
254:         for (i = 0; i <= lmvm->k - 1; ++i) {
255:           lsb->yts[i] = lsb->yts[i + 1];
256:           lsb->yty[i] = lsb->yty[i + 1];
257:           lsb->sts[i] = lsb->sts[i + 1];
258:         }
259:       }
260:       /* Update history of useful scalars */
261:       VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);
262:       lsb->yts[lmvm->k] = PetscRealPart(curvature);
263:       lsb->yty[lmvm->k] = PetscRealPart(ytytmp);
264:       lsb->sts[lmvm->k] = PetscRealPart(ststmp);
265:       /* Compute the scalar scale if necessary */
266:       if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) MatSymBrdnComputeJ0Scalar(B);
267:     } else {
268:       /* Update is bad, skip it */
269:       ++lmvm->nrejects;
270:       ++lsb->watchdog;
271:     }
272:   } else {
273:     switch (lsb->scale_type) {
274:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
275:       dbase = (Mat_LMVM *)lsb->D->data;
276:       dctx  = (Mat_DiagBrdn *)dbase->ctx;
277:       VecSet(dctx->invD, lsb->delta);
278:       break;
279:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
280:       lsb->sigma = lsb->delta;
281:       break;
282:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
283:       lsb->sigma = 1.0;
284:       break;
285:     default:
286:       break;
287:     }
288:   }

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

293:   if (lsb->watchdog > lsb->max_seq_rejects) {
294:     MatLMVMReset(B, PETSC_FALSE);
295:     if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) MatLMVMReset(lsb->D, PETSC_FALSE);
296:   }

298:   /* Save the solution and function to be used in the next update */
299:   VecCopy(X, lmvm->Xprev);
300:   VecCopy(F, lmvm->Fprev);
301:   lmvm->prev_set = PETSC_TRUE;
302:   return 0;
303: }

305: /*------------------------------------------------------------*/

307: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
308: {
309:   Mat_LMVM    *bdata = (Mat_LMVM *)B->data;
310:   Mat_SymBrdn *blsb  = (Mat_SymBrdn *)bdata->ctx;
311:   Mat_LMVM    *mdata = (Mat_LMVM *)M->data;
312:   Mat_SymBrdn *mlsb  = (Mat_SymBrdn *)mdata->ctx;
313:   PetscInt     i;

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

350: /*------------------------------------------------------------*/

352: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
353: {
354:   Mat_LMVM     *lmvm = (Mat_LMVM *)B->data;
355:   Mat_SymBrdn  *lsb  = (Mat_SymBrdn *)lmvm->ctx;
356:   Mat_LMVM     *dbase;
357:   Mat_DiagBrdn *dctx;

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

400: /*------------------------------------------------------------*/

402: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
403: {
404:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
405:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

407:   MatAllocate_LMVM(B, X, F);
408:   if (!lsb->allocated) {
409:     VecDuplicate(X, &lsb->work);
410:     if (lmvm->m > 0) {
411:       PetscMalloc5(lmvm->m, &lsb->stp, lmvm->m, &lsb->ytq, lmvm->m, &lsb->yts, lmvm->m, &lsb->yty, lmvm->m, &lsb->sts);
412:       PetscCalloc1(lmvm->m, &lsb->psi);
413:       VecDuplicateVecs(X, lmvm->m, &lsb->P);
414:       VecDuplicateVecs(X, lmvm->m, &lsb->Q);
415:     }
416:     switch (lsb->scale_type) {
417:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
418:       MatLMVMAllocate(lsb->D, X, F);
419:       break;
420:     default:
421:       break;
422:     }
423:     lsb->allocated = PETSC_TRUE;
424:   }
425:   return 0;
426: }

428: /*------------------------------------------------------------*/

430: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
431: {
432:   Mat_LMVM    *lmvm = (Mat_LMVM *)B->data;
433:   Mat_SymBrdn *lsb  = (Mat_SymBrdn *)lmvm->ctx;

435:   if (lsb->allocated) {
436:     VecDestroy(&lsb->work);
437:     PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);
438:     PetscFree(lsb->psi);
439:     VecDestroyVecs(lmvm->m, &lsb->P);
440:     VecDestroyVecs(lmvm->m, &lsb->Q);
441:     lsb->allocated = PETSC_FALSE;
442:   }
443:   MatDestroy(&lsb->D);
444:   PetscFree(lmvm->ctx);
445:   MatDestroy_LMVM(B);
446:   return 0;
447: }

449: /*------------------------------------------------------------*/

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

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

481: /*------------------------------------------------------------*/

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

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

501: /*------------------------------------------------------------*/

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

508:   MatSetFromOptions_LMVM(B, PetscOptionsObject);
509:   PetscOptionsHeadBegin(PetscOptionsObject, "Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
510:   PetscOptionsReal("-mat_lmvm_phi", "(developer) convex ratio between BFGS and DFP components of the update", "", lsb->phi, &lsb->phi, NULL);
512:   MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject);
513:   PetscOptionsHeadEnd();
514:   return 0;
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:   PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", lsb->beta, &lsb->beta, NULL);
527:   PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", lsb->theta, &lsb->theta, NULL);
529:   PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", lsb->rho, &lsb->rho, NULL);
531:   PetscOptionsReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", lsb->alpha, &lsb->alpha, NULL);
533:   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);
534:   PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBrdnScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg);
535:   if (flg) MatLMVMSymBroydenSetScaleType(B, stype);
536:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
537:     const char *prefix;

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

556: /*------------------------------------------------------------*/

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

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

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

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

599:   MatCreate(PetscObjectComm((PetscObject)B), &lsb->D);
600:   MatSetType(lsb->D, MATLMVMDIAGBROYDEN);
601:   return 0;
602: }

604: /*------------------------------------------------------------*/

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

610:    Input Parameters:
611: +  B - LMVM matrix
612: -  delta - initial value for diagonal scaling

614:    Level: intermediate

616: @*/

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

624:   PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);
625:   PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);
626:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn);
627:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn);
629:   lsb->delta = PetscAbsReal(PetscRealPart(delta));
630:   lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
631:   lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
632:   return 0;
633: }

635: /*------------------------------------------------------------*/

637: /*@
638:     MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.

640:     Input Parameters:
641: +   snes - the iterative context
642: -   rtype - restart type

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

647:     Level: intermediate

649:     MatLMVMSymBrdnScaleTypes:
650: +   MAT_LMVM_SYMBROYDEN_SCALE_NONE - initial Hessian is the identity matrix
651: .   MAT_LMVM_SYMBROYDEN_SCALE_SCALAR - use the Shanno scalar as the initial Hessian
652: -   MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL - use a diagonalized BFGS update as the initial Hessian

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

662:   lsb->scale_type = stype;
663:   return 0;
664: }

666: /*------------------------------------------------------------*/

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

675:    The provided local and global sizes must match the solution and function vectors
676:    used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have
677:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
678:    parallel. 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:    Collective

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

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

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

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

705:    Level: intermediate

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

719: /*------------------------------------------------------------*/

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

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

747: /*------------------------------------------------------------*/

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

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

775: /*------------------------------------------------------------*/

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

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