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:   PetscErrorCode    ierr;
 37:   PetscInt          i, j;
 38:   PetscReal         numer;
 39:   PetscScalar       sjtpi, yjtsi, wtsi, yjtqi, sjtyi, wtyi, ytx, stf, wtf, stp, ytq;

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

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

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

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

128:   return(0);
129: }

131: /*------------------------------------------------------------*/

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

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

143:   Z <- J0 * X

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

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

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

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

177:   VecCheckSameSize(X, 2, Z, 3);
178:   VecCheckMatCompatible(B, X, 2, Z, 3);

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

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

223: /*------------------------------------------------------------*/

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

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

298:   /* Update the scaling */
299:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
300:     MatLMVMUpdate(lsb->D, X, F);
301:   }

303:   if (lsb->watchdog > lsb->max_seq_rejects) {
304:     MatLMVMReset(B, PETSC_FALSE);
305:     if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
306:       MatLMVMReset(lsb->D, PETSC_FALSE);
307:     }
308:   }

310:   /* Save the solution and function to be used in the next update */
311:   VecCopy(X, lmvm->Xprev);
312:   VecCopy(F, lmvm->Fprev);
313:   lmvm->prev_set = PETSC_TRUE;
314:   return(0);
315: }

317: /*------------------------------------------------------------*/

319: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
320: {
321:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
322:   Mat_SymBrdn       *blsb = (Mat_SymBrdn*)bdata->ctx;
323:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
324:   Mat_SymBrdn       *mlsb = (Mat_SymBrdn*)mdata->ctx;
325:   PetscErrorCode    ierr;
326:   PetscInt          i;

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

364: /*------------------------------------------------------------*/

366: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
367: {
368:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
369:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
370:   Mat_LMVM          *dbase;
371:   Mat_DiagBrdn      *dctx;
372:   PetscErrorCode    ierr;

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

416: /*------------------------------------------------------------*/

418: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
419: {
420:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
421:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
422:   PetscErrorCode    ierr;

425:   MatAllocate_LMVM(B, X, F);
426:   if (!lsb->allocated) {
427:     VecDuplicate(X, &lsb->work);
428:     if (lmvm->m > 0) {
429:       PetscMalloc5(lmvm->m,&lsb->stp,lmvm->m,&lsb->ytq,lmvm->m,&lsb->yts,lmvm->m,&lsb->yty,lmvm->m,&lsb->sts);
430:       PetscCalloc1(lmvm->m,&lsb->psi);
431:       VecDuplicateVecs(X, lmvm->m, &lsb->P);
432:       VecDuplicateVecs(X, lmvm->m, &lsb->Q);
433:     }
434:     switch (lsb->scale_type) {
435:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
436:       MatLMVMAllocate(lsb->D, X, F);
437:       break;
438:     default:
439:       break;
440:     }
441:     lsb->allocated = PETSC_TRUE;
442:   }
443:   return(0);
444: }

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

448: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
449: {
450:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
451:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
452:   PetscErrorCode    ierr;

455:   if (lsb->allocated) {
456:     VecDestroy(&lsb->work);
457:     PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);
458:     PetscFree(lsb->psi);
459:     VecDestroyVecs(lmvm->m, &lsb->P);
460:     VecDestroyVecs(lmvm->m, &lsb->Q);
461:     lsb->allocated = PETSC_FALSE;
462:   }
463:   MatDestroy(&lsb->D);
464:   PetscFree(lmvm->ctx);
465:   MatDestroy_LMVM(B);
466:   return(0);
467: }

469: /*------------------------------------------------------------*/

471: static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
472: {
473:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
474:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
475:   PetscErrorCode    ierr;
476:   PetscInt          n, N;

479:   MatSetUp_LMVM(B);
480:   if (!lsb->allocated) {
481:     VecDuplicate(lmvm->Xprev, &lsb->work);
482:     if (lmvm->m > 0) {
483:       PetscMalloc5(lmvm->m,&lsb->stp,lmvm->m,&lsb->ytq,lmvm->m,&lsb->yts,lmvm->m,&lsb->yty,lmvm->m,&lsb->sts);
484:       PetscCalloc1(lmvm->m,&lsb->psi);
485:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->P);
486:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->Q);
487:     }
488:     switch (lsb->scale_type) {
489:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
490:       MatGetLocalSize(B, &n, &n);
491:       MatGetSize(B, &N, &N);
492:       MatSetSizes(lsb->D, n, n, N, N);
493:       MatSetUp(lsb->D);
494:       break;
495:     default:
496:       break;
497:     }
498:     lsb->allocated = PETSC_TRUE;
499:   }
500:   return(0);
501: }

503: /*------------------------------------------------------------*/

505: PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
506: {
507:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
508:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
509:   PetscErrorCode    ierr;
510:   PetscBool         isascii;

513:   PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
514:   if (isascii) {
515:     PetscViewerASCIIPrintf(pv,"Scale type: %s\n",MatLMVMSymBroydenScaleTypes[lsb->scale_type]);
516:     PetscViewerASCIIPrintf(pv,"Scale history: %d\n",lsb->sigma_hist);
517:     PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)lsb->alpha, (double)lsb->beta, (double)lsb->rho);
518:     PetscViewerASCIIPrintf(pv,"Convex factors: phi=%g, theta=%g\n",(double)lsb->phi, (double)lsb->theta);
519:   }
520:   MatView_LMVM(B, pv);
521:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
522:     MatView(lsb->D, pv);
523:   }
524:   return(0);
525: }

527: /*------------------------------------------------------------*/

529: PetscErrorCode MatSetFromOptions_LMVMSymBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
530: {
531:   Mat_LMVM                     *lmvm = (Mat_LMVM*)B->data;
532:   Mat_SymBrdn                  *lsb = (Mat_SymBrdn*)lmvm->ctx;
533:   PetscErrorCode               ierr;

536:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
537:   PetscOptionsHead(PetscOptionsObject,"Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
538:   PetscOptionsReal("-mat_lmvm_phi","(developer) convex ratio between BFGS and DFP components of the update","",lsb->phi,&lsb->phi,NULL);
539:   if ((lsb->phi < 0.0) || (lsb->phi > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the update formula cannot be outside the range of [0, 1]");
540:   MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
541:   PetscOptionsTail();
542:   return(0);
543: }

545: PetscErrorCode MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionItems *PetscOptionsObject, Mat B)
546: {
547:   Mat_LMVM                     *lmvm = (Mat_LMVM*)B->data;
548:   Mat_SymBrdn                  *lsb = (Mat_SymBrdn*)lmvm->ctx;
549:   Mat_LMVM                     *dbase;
550:   Mat_DiagBrdn                 *dctx;
551:   MatLMVMSymBroydenScaleType   stype = lsb->scale_type;
552:   PetscBool                    flg;
553:   PetscErrorCode               ierr;

556:   PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",lsb->beta,&lsb->beta,NULL);
557:   PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",lsb->theta,&lsb->theta,NULL);
558:   if ((lsb->theta < 0.0) || (lsb->theta > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
559:   PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",lsb->rho,&lsb->rho,NULL);
560:   if ((lsb->rho < 0.0) || (lsb->rho > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "update limiter in the J0 scaling cannot be outside the range of [0, 1]");
561:   PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",lsb->alpha,&lsb->alpha,NULL);
562:   if ((lsb->alpha < 0.0) || (lsb->alpha > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
563:   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);
564:   PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0","MatLMVMSymBrdnScaleType",MatLMVMSymBroydenScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);
565:   if (flg) {MatLMVMSymBroydenSetScaleType(B, stype);}
566:   if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
567:     MatSetFromOptions(lsb->D);
568:     dbase = (Mat_LMVM*)lsb->D->data;
569:     dctx = (Mat_DiagBrdn*)dbase->ctx;
570:     dctx->delta_min  = lsb->delta_min;
571:     dctx->delta_max  = lsb->delta_max;
572:     dctx->theta      = lsb->theta;
573:     dctx->rho        = lsb->rho;
574:     dctx->alpha      = lsb->alpha;
575:     dctx->beta       = lsb->beta;
576:     dctx->sigma_hist = lsb->sigma_hist;
577:   }
578:   return(0);
579: }

581: /*------------------------------------------------------------*/

583: PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
584: {
585:   Mat_LMVM          *lmvm;
586:   Mat_SymBrdn       *lsb;
587:   PetscErrorCode    ierr;

590:   MatCreate_LMVM(B);
591:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN);
592:   MatSetOption(B, MAT_SPD, PETSC_TRUE);
593:   B->ops->view = MatView_LMVMSymBrdn;
594:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
595:   B->ops->setup = MatSetUp_LMVMSymBrdn;
596:   B->ops->destroy = MatDestroy_LMVMSymBrdn;
597:   B->ops->solve = MatSolve_LMVMSymBrdn;

599:   lmvm = (Mat_LMVM*)B->data;
600:   lmvm->square = PETSC_TRUE;
601:   lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
602:   lmvm->ops->reset = MatReset_LMVMSymBrdn;
603:   lmvm->ops->update = MatUpdate_LMVMSymBrdn;
604:   lmvm->ops->mult = MatMult_LMVMSymBrdn;
605:   lmvm->ops->copy = MatCopy_LMVMSymBrdn;

607:   PetscNewLog(B, &lsb);
608:   lmvm->ctx = (void*)lsb;
609:   lsb->allocated       = PETSC_FALSE;
610:   lsb->needP           = lsb->needQ = PETSC_TRUE;
611:   lsb->phi             = 0.125;
612:   lsb->theta           = 0.125;
613:   lsb->alpha           = 1.0;
614:   lsb->rho             = 1.0;
615:   lsb->beta            = 0.5;
616:   lsb->sigma           = 1.0;
617:   lsb->delta           = 1.0;
618:   lsb->delta_min       = 1e-7;
619:   lsb->delta_max       = 100.0;
620:   lsb->sigma_hist      = 1;
621:   lsb->scale_type      = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
622:   lsb->watchdog        = 0;
623:   lsb->max_seq_rejects = lmvm->m/2;

625:   MatCreate(PetscObjectComm((PetscObject)B), &lsb->D);
626:   MatSetType(lsb->D, MATLMVMDIAGBROYDEN);
627:   MatSetOptionsPrefix(lsb->D, "J0_");
628:   return(0);
629: }

631: /*------------------------------------------------------------*/

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

637:    Input Parameters:
638: +  B - LMVM matrix
639: -  delta - initial value for diagonal scaling

641:    Level: intermediate
642: @*/

644: PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
645: {
646:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
647:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
648:   PetscErrorCode    ierr;
649:   PetscBool         is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn;

652:   PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);
653:   PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);
654:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn);
655:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn);
656:   if (!is_bfgs && !is_dfp && !is_symbrdn && !is_symbadbrdn) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "diagonal scaling is only available for DFP, BFGS and SymBrdn matrices");
657:   lsb->delta = PetscAbsReal(PetscRealPart(delta));
658:   lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
659:   lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
660:   return(0);
661: }

663: /*------------------------------------------------------------*/

665: /*@
666:     MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.

668:     Input Parameters:
669: +   snes - the iterative context
670: -   rtype - restart type

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

675:     Level: intermediate

677:     MatLMVMSymBrdnScaleTypes:
678: +   MAT_LMVM_SYMBROYDEN_SCALE_NONE - initial Hessian is the identity matrix
679: .   MAT_LMVM_SYMBROYDEN_SCALE_SCALAR - use the Shanno scalar as the initial Hessian
680: -   MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL - use a diagonalized BFGS update as the initial Hessian

682: .seealso: MATLMVMSYMBROYDEN, MatCreateLMVMSymBroyden()
683: @*/
684: PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
685: {
686:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
687:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;

691:   lsb->scale_type = stype;
692:   return(0);
693: }

695: /*------------------------------------------------------------*/

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

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

712:    Collective

714:    Input Parameters:
715: +  comm - MPI communicator, set to PETSC_COMM_SELF
716: .  n - number of local rows for storage vectors
717: -  N - global size of the storage vectors

719:    Output Parameter:
720: .  B - the matrix

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

725:    Options Database Keys:
726: +   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
727: .   -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
728: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
729: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
730: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
731: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
732: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
733: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

735:    Level: intermediate

737: .seealso: MatCreate(), MATLMVM, MATLMVMSYMBROYDEN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
738:           MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn()
739: @*/
740: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
741: {
742:   PetscErrorCode    ierr;

745:   MatCreate(comm, B);
746:   MatSetSizes(*B, n, n, N, N);
747:   MatSetType(*B, MATLMVMSYMBROYDEN);
748:   MatSetUp(*B);
749:   return(0);
750: }

752: /*------------------------------------------------------------*/

754: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
755: {
756:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
757:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
758:   PetscErrorCode    ierr;

761:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
762:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
763:     MatLMVMApplyJ0Fwd(B, X, Z);
764:   } else {
765:     switch (lsb->scale_type) {
766:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
767:       VecCopy(X, Z);
768:       VecScale(Z, 1.0/lsb->sigma);
769:       break;
770:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
771:       MatMult(lsb->D, X, Z);
772:       break;
773:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
774:     default:
775:       VecCopy(X, Z);
776:       break;
777:     }
778:   }
779:   return(0);
780: }

782: /*------------------------------------------------------------*/

784: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
785: {
786:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
787:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
788:   PetscErrorCode    ierr;

791:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
792:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
793:     MatLMVMApplyJ0Inv(B, F, dX);
794:   } else {
795:     switch (lsb->scale_type) {
796:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
797:       VecCopy(F, dX);
798:       VecScale(dX, lsb->sigma);
799:       break;
800:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
801:       MatSolve(lsb->D, F, dX);
802:       break;
803:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
804:     default:
805:       VecCopy(F, dX);
806:       break;
807:     }
808:   }
809:   return(0);
810: }

812: /*------------------------------------------------------------*/

814: PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
815: {
816:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
817:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
818:   PetscInt          i, start;
819:   PetscReal         a, b, c, sig1, sig2, signew;

822:   if (lsb->sigma_hist == 0) {
823:     signew = 1.0;
824:   } else {
825:     start = PetscMax(0, lmvm->k-lsb->sigma_hist+1);
826:     signew = 0.0;
827:     if (lsb->alpha == 1.0) {
828:       for (i = start; i <= lmvm->k; ++i) {
829:         signew += lsb->yts[i]/lsb->yty[i];
830:       }
831:     } else if (lsb->alpha == 0.5) {
832:       for (i = start; i <= lmvm->k; ++i) {
833:         signew += lsb->sts[i]/lsb->yty[i];
834:       }
835:       signew = PetscSqrtReal(signew);
836:     } else if (lsb->alpha == 0.0) {
837:       for (i = start; i <= lmvm->k; ++i) {
838:         signew += lsb->sts[i]/lsb->yts[i];
839:       }
840:     } else {
841:       /* compute coefficients of the quadratic */
842:       a = b = c = 0.0;
843:       for (i = start; i <= lmvm->k; ++i) {
844:         a += lsb->yty[i];
845:         b += lsb->yts[i];
846:         c += lsb->sts[i];
847:       }
848:       a *= lsb->alpha;
849:       b *= -(2.0*lsb->alpha - 1.0);
850:       c *= lsb->alpha - 1.0;
851:       /* use quadratic formula to find roots */
852:       sig1 = (-b + PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
853:       sig2 = (-b - PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
854:       /* accept the positive root as the scalar */
855:       if (sig1 > 0.0) {
856:         signew = sig1;
857:       } else if (sig2 > 0.0) {
858:         signew = sig2;
859:       } else {
860:         SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
861:       }
862:     }
863:   }
864:   lsb->sigma = lsb->rho*signew + (1.0 - lsb->rho)*lsb->sigma;
865:   return(0);
866: }