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) {
267:         MatSymBrdnComputeJ0Scalar(B);
268:       }
269:     } else {
270:       /* Update is bad, skip it */
271:       ++lmvm->nrejects;
272:       ++lsb->watchdog;
273:     }
274:   } else {
275:     switch (lsb->scale_type) {
276:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
277:       dbase = (Mat_LMVM*)lsb->D->data;
278:       dctx = (Mat_DiagBrdn*)dbase->ctx;
279:       VecSet(dctx->invD, lsb->delta);
280:       break;
281:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
282:       lsb->sigma = lsb->delta;
283:       break;
284:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
285:       lsb->sigma = 1.0;
286:       break;
287:     default:
288:       break;
289:     }
290:   }

292:   /* Update the scaling */
293:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
294:     MatLMVMUpdate(lsb->D, X, F);
295:   }

297:   if (lsb->watchdog > lsb->max_seq_rejects) {
298:     MatLMVMReset(B, PETSC_FALSE);
299:     if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
300:       MatLMVMReset(lsb->D, PETSC_FALSE);
301:     }
302:   }

304:   /* Save the solution and function to be used in the next update */
305:   VecCopy(X, lmvm->Xprev);
306:   VecCopy(F, lmvm->Fprev);
307:   lmvm->prev_set = PETSC_TRUE;
308:   return 0;
309: }

311: /*------------------------------------------------------------*/

313: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
314: {
315:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
316:   Mat_SymBrdn       *blsb = (Mat_SymBrdn*)bdata->ctx;
317:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
318:   Mat_SymBrdn       *mlsb = (Mat_SymBrdn*)mdata->ctx;
319:   PetscInt          i;

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

356: /*------------------------------------------------------------*/

358: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
359: {
360:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
361:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
362:   Mat_LMVM          *dbase;
363:   Mat_DiagBrdn      *dctx;

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

406: /*------------------------------------------------------------*/

408: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
409: {
410:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
411:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;

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

434: /*------------------------------------------------------------*/

436: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
437: {
438:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
439:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;

441:   if (lsb->allocated) {
442:     VecDestroy(&lsb->work);
443:     PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);
444:     PetscFree(lsb->psi);
445:     VecDestroyVecs(lmvm->m, &lsb->P);
446:     VecDestroyVecs(lmvm->m, &lsb->Q);
447:     lsb->allocated = PETSC_FALSE;
448:   }
449:   MatDestroy(&lsb->D);
450:   PetscFree(lmvm->ctx);
451:   MatDestroy_LMVM(B);
452:   return 0;
453: }

455: /*------------------------------------------------------------*/

457: static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
458: {
459:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
460:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
461:   PetscInt          n, N;

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

487: /*------------------------------------------------------------*/

489: PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
490: {
491:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
492:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
493:   PetscBool         isascii;

495:   PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
496:   if (isascii) {
497:     PetscViewerASCIIPrintf(pv,"Scale type: %s\n",MatLMVMSymBroydenScaleTypes[lsb->scale_type]);
498:     PetscViewerASCIIPrintf(pv,"Scale history: %d\n",lsb->sigma_hist);
499:     PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)lsb->alpha, (double)lsb->beta, (double)lsb->rho);
500:     PetscViewerASCIIPrintf(pv,"Convex factors: phi=%g, theta=%g\n",(double)lsb->phi, (double)lsb->theta);
501:   }
502:   MatView_LMVM(B, pv);
503:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
504:     MatView(lsb->D, pv);
505:   }
506:   return 0;
507: }

509: /*------------------------------------------------------------*/

511: PetscErrorCode MatSetFromOptions_LMVMSymBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
512: {
513:   Mat_LMVM                     *lmvm = (Mat_LMVM*)B->data;
514:   Mat_SymBrdn                  *lsb = (Mat_SymBrdn*)lmvm->ctx;

516:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
517:   PetscOptionsHead(PetscOptionsObject,"Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
518:   PetscOptionsReal("-mat_lmvm_phi","(developer) convex ratio between BFGS and DFP components of the update","",lsb->phi,&lsb->phi,NULL);
520:   MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
521:   PetscOptionsTail();
522:   return 0;
523: }

525: PetscErrorCode MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionItems *PetscOptionsObject, Mat B)
526: {
527:   Mat_LMVM                     *lmvm = (Mat_LMVM*)B->data;
528:   Mat_SymBrdn                  *lsb = (Mat_SymBrdn*)lmvm->ctx;
529:   Mat_LMVM                     *dbase;
530:   Mat_DiagBrdn                 *dctx;
531:   MatLMVMSymBroydenScaleType   stype = lsb->scale_type;
532:   PetscBool                    flg;

534:   PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",lsb->beta,&lsb->beta,NULL);
535:   PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",lsb->theta,&lsb->theta,NULL);
537:   PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",lsb->rho,&lsb->rho,NULL);
539:   PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",lsb->alpha,&lsb->alpha,NULL);
541:   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);
542:   PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0","MatLMVMSymBrdnScaleType",MatLMVMSymBroydenScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);
543:   if (flg) MatLMVMSymBroydenSetScaleType(B, stype);
544:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
545:     MatSetFromOptions(lsb->D);
546:     dbase = (Mat_LMVM*)lsb->D->data;
547:     dctx = (Mat_DiagBrdn*)dbase->ctx;
548:     dctx->delta_min  = lsb->delta_min;
549:     dctx->delta_max  = lsb->delta_max;
550:     dctx->theta      = lsb->theta;
551:     dctx->rho        = lsb->rho;
552:     dctx->alpha      = lsb->alpha;
553:     dctx->beta       = lsb->beta;
554:     dctx->sigma_hist = lsb->sigma_hist;
555:   }
556:   return 0;
557: }

559: /*------------------------------------------------------------*/

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

566:   MatCreate_LMVM(B);
567:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN);
568:   MatSetOption(B, MAT_SPD, 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:   PetscNewLog(B, &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:   MatCreate(PetscObjectComm((PetscObject)B), &lsb->D);
602:   MatSetType(lsb->D, MATLMVMDIAGBROYDEN);
603:   MatSetOptionsPrefix(lsb->D, "J0_");
604:   return 0;
605: }

607: /*------------------------------------------------------------*/

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

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

617:    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:   PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);
627:   PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);
628:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn);
629:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn);
631:   lsb->delta = PetscAbsReal(PetscRealPart(delta));
632:   lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
633:   lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
634:   return 0;
635: }

637: /*------------------------------------------------------------*/

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

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

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

649:     Level: intermediate

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

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

664:   lsb->scale_type = stype;
665:   return 0;
666: }

668: /*------------------------------------------------------------*/

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

677:    The provided local and global sizes must match the solution and function vectors
678:    used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have
679:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
680:    parallel. To use the L-SymBrdn matrix with other vector types, the matrix must be
681:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
682:    This ensures that the internal storage and work vectors are duplicated from the
683:    correct type of vector.

685:    Collective

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

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

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

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

708:    Level: intermediate

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

722: /*------------------------------------------------------------*/

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

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

750: /*------------------------------------------------------------*/

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

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

778: /*------------------------------------------------------------*/

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

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) {
794:         signew += lsb->yts[i]/lsb->yty[i];
795:       }
796:     } else if (lsb->alpha == 0.5) {
797:       for (i = start; i <= lmvm->k; ++i) {
798:         signew += lsb->sts[i]/lsb->yty[i];
799:       }
800:       signew = PetscSqrtReal(signew);
801:     } else if (lsb->alpha == 0.0) {
802:       for (i = start; i <= lmvm->k; ++i) {
803:         signew += lsb->sts[i]/lsb->yts[i];
804:       }
805:     } else {
806:       /* compute coefficients of the quadratic */
807:       a = b = c = 0.0;
808:       for (i = start; i <= lmvm->k; ++i) {
809:         a += lsb->yty[i];
810:         b += lsb->yts[i];
811:         c += lsb->sts[i];
812:       }
813:       a *= lsb->alpha;
814:       b *= -(2.0*lsb->alpha - 1.0);
815:       c *= lsb->alpha - 1.0;
816:       /* use quadratic formula to find roots */
817:       sig1 = (-b + PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
818:       sig2 = (-b - PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
819:       /* accept the positive root as the scalar */
820:       if (sig1 > 0.0) {
821:         signew = sig1;
822:       } else if (sig2 > 0.0) {
823:         signew = sig2;
824:       } else {
825:         SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
826:       }
827:     }
828:   }
829:   lsb->sigma = lsb->rho*signew + (1.0 - lsb->rho)*lsb->sigma;
830:   return 0;
831: }