Actual source code: symbrdn.c

petsc-3.12.2 2019-11-22
Report Typos and Errors
  1:  #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
  2:  #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>

  4: /*------------------------------------------------------------*/

  6: /*
  7:   The solution method below is the matrix-free implementation of 
  8:   Equation 8.6a in Dennis and More "Quasi-Newton Methods, Motivation 
  9:   and Theory" (https://epubs.siam.org/doi/abs/10.1137/1019005).
 10:   
 11:   Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever 
 12:   the matrix is updated with a new (S[i], Y[i]) pair. This allows 
 13:   repeated calls of MatSolve without incurring redundant computation.
 14:   
 15:   dX <- J0^{-1} * F
 16:   
 17:   for i=0,1,2,...,k
 18:     # Q[i] = (B_i)^T{-1} Y[i]
 19:     
 20:     rho = 1.0 / (Y[i]^T S[i])
 21:     alpha = rho * (S[i]^T F)
 22:     zeta = 1.0 / (Y[i]^T Q[i])
 23:     gamma = zeta * (Y[i]^T dX)
 24:     
 25:     dX <- dX - (gamma * Q[i]) + (alpha * Y[i])
 26:     W <- (rho * S[i]) - (zeta * Q[i])
 27:     dX <- dX + (psi[i] * (Y[i]^T Q[i]) * (W^T F) * W)
 28:   end
 29: */
 30: static PetscErrorCode MatSolve_LMVMSymBrdn(Mat B, Vec F, Vec dX)
 31: {
 32:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
 33:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
 34:   PetscErrorCode    ierr;
 35:   PetscInt          i, j;
 36:   PetscReal         numer;
 37:   PetscScalar       sjtpi, yjtsi, wtsi, yjtqi, sjtyi, wtyi, ytx, stf, wtf, stp, ytq;
 38: 
 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:   }
 49: 
 50:   VecCheckSameSize(F, 2, dX, 3);
 51:   VecCheckMatCompatible(B, dX, 3, F, 2);
 52: 
 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:   }
109: 
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).
135:   
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.
140:   
141:   Z <- J0 * X
142:   
143:   for i=0,1,2,...,k
144:     # P[i] = (B_k) * S[i]
145:     
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)
150:     
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:   PetscErrorCode    ierr;
161:   PetscInt          i, j;
162:   PetscScalar         sjtpi, yjtsi, wtsi, stz, ytx, wtx, stp;
163: 
164: 
166:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
167:   if (lsb->phi == 0.0) {
168:     MatMult_LMVMBFGS(B, X, Z);
169:     return(0);
170:   }
171:   if (lsb->phi == 1.0) {
172:     MatMult_LMVMDFP(B, X, Z);
173:     return(0);
174:   }
175: 
176:   VecCheckSameSize(X, 2, Z, 3);
177:   VecCheckMatCompatible(B, X, 2, Z, 3);
178: 
179:   if (lsb->needP) {
180:     /* Start the loop for (P[k] = (B_k) * S[k]) */
181:     for (i = 0; i <= lmvm->k; ++i) {
182:       MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);
183:       for (j = 0; j <= i-1; ++j) {
184:         /* Compute the necessary dot products */
185:         VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);
186:         VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
187:         VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);
188:         VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
189:         /* Compute the pure BFGS component of the forward product */
190:         VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi)/lsb->stp[j], PetscRealPart(yjtsi)/lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);
191:         /* Tack on the convexly scaled extras to the forward product */
192:         if (lsb->phi > 0.0) {
193:           VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);
194:           VecDot(lsb->work, lmvm->S[i], &wtsi);
195:           VecAXPY(lsb->P[i], lsb->phi*lsb->stp[j]*PetscRealPart(wtsi), lsb->work);
196:         }
197:       }
198:       VecDot(lmvm->S[i], lsb->P[i], &stp);
199:       lsb->stp[i] = PetscRealPart(stp);
200:     }
201:     lsb->needP = PETSC_FALSE;
202:   }
203: 
204:   /* Start the outer iterations for (B * X) */
205:   MatSymBrdnApplyJ0Fwd(B, X, Z);
206:   for (i = 0; i <= lmvm->k; ++i) {
207:     /* Compute the necessary dot products */
208:     VecDotBegin(lmvm->S[i], Z, &stz);
209:     VecDotBegin(lmvm->Y[i], X, &ytx);
210:     VecDotEnd(lmvm->S[i], Z, &stz);
211:     VecDotEnd(lmvm->Y[i], X, &ytx);
212:     /* Compute the pure BFGS component */
213:     VecAXPBYPCZ(Z, -PetscRealPart(stz)/lsb->stp[i], PetscRealPart(ytx)/lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]);
214:     /* Tack on the convexly scaled extras */
215:     VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]);
216:     VecDot(lsb->work, X, &wtx);
217:     VecAXPY(Z, lsb->phi*lsb->stp[i]*PetscRealPart(wtx), lsb->work);
218:   }
219:   return(0);
220: }

222: /*------------------------------------------------------------*/

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

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

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

316: /*------------------------------------------------------------*/

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

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

363: /*------------------------------------------------------------*/

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

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

415: /*------------------------------------------------------------*/

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

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

445: /*------------------------------------------------------------*/

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

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

468: /*------------------------------------------------------------*/

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

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

502: /*------------------------------------------------------------*/

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

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

526: /*------------------------------------------------------------*/

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

537:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
538:   PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
539:   PetscOptionsEList("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "", Scale_Table, SYMBRDN_SCALE_SIZE, Scale_Table[lsb->scale_type], &lsb->scale_type,NULL);
540:   PetscOptionsReal("-mat_lmvm_phi","(developer) convex ratio between BFGS and DFP components of the update","",lsb->phi,&lsb->phi,NULL);
541:   PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",lsb->theta,&lsb->theta,NULL);
542:   PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",lsb->rho,&lsb->rho,NULL);
543:   PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",lsb->alpha,&lsb->alpha,NULL);
544:   PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",lsb->beta,&lsb->beta,NULL);
545:   PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",lsb->sigma_hist,&lsb->sigma_hist,NULL);
546:   PetscOptionsTail();
547:   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]");
548:   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]");
549:   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]");
550:   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]");
551:   if (lsb->sigma_hist < 0) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
552:   if (lsb->scale_type == SYMBRDN_SCALE_DIAG) {
553:     MatSetFromOptions(lsb->D);
554:     dbase = (Mat_LMVM*)lsb->D->data;
555:     dctx = (Mat_DiagBrdn*)dbase->ctx;
556:     dctx->delta_min  = lsb->delta_min;
557:     dctx->delta_max  = lsb->delta_max;
558:     dctx->theta      = lsb->theta;
559:     dctx->rho        = lsb->rho;
560:     dctx->alpha      = lsb->alpha;
561:     dctx->beta       = lsb->beta;
562:     dctx->sigma_hist = lsb->sigma_hist;
563:     dctx->forward    = PETSC_TRUE;
564:   }
565:   return(0);
566: }

568: /*------------------------------------------------------------*/

570: PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
571: {
572:   Mat_LMVM          *lmvm;
573:   Mat_SymBrdn       *lsb;
574:   PetscErrorCode    ierr;

577:   MatCreate_LMVM(B);
578:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBRDN);
579:   MatSetOption(B, MAT_SPD, PETSC_TRUE);
580:   B->ops->view = MatView_LMVMSymBrdn;
581:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
582:   B->ops->setup = MatSetUp_LMVMSymBrdn;
583:   B->ops->destroy = MatDestroy_LMVMSymBrdn;
584:   B->ops->solve = MatSolve_LMVMSymBrdn;
585: 
586:   lmvm = (Mat_LMVM*)B->data;
587:   lmvm->square = PETSC_TRUE;
588:   lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
589:   lmvm->ops->reset = MatReset_LMVMSymBrdn;
590:   lmvm->ops->update = MatUpdate_LMVMSymBrdn;
591:   lmvm->ops->mult = MatMult_LMVMSymBrdn;
592:   lmvm->ops->copy = MatCopy_LMVMSymBrdn;
593: 
594:   PetscNewLog(B, &lsb);
595:   lmvm->ctx = (void*)lsb;
596:   lsb->allocated       = PETSC_FALSE;
597:   lsb->needP           = lsb->needQ = PETSC_TRUE;
598:   lsb->phi             = 0.125;
599:   lsb->theta           = 0.125;
600:   lsb->alpha           = 1.0;
601:   lsb->rho             = 1.0;
602:   lsb->beta            = 0.5;
603:   lsb->sigma           = 1.0;
604:   lsb->delta           = 1.0;
605:   lsb->delta_min       = 1e-7;
606:   lsb->delta_max       = 100.0;
607:   lsb->sigma_hist      = 1;
608:   lsb->scale_type      = SYMBRDN_SCALE_DIAG;
609:   lsb->watchdog        = 0;
610:   lsb->max_seq_rejects = lmvm->m/2;
611: 
612:   MatCreate(PetscObjectComm((PetscObject)B), &lsb->D);
613:   MatSetType(lsb->D, MATLMVMDIAGBRDN);
614:   MatSetOptionsPrefix(lsb->D, "J0_");
615:   return(0);
616: }

618: /*------------------------------------------------------------*/

620: /*@
621:    MatSymBrdnSetDelta - Sets the starting value for the diagonal scaling vector computed 
622:    in the SymBrdn approximations (also works for BFGS and DFP).
623:    
624:    Input Parameters:
625: +  B - LMVM matrix
626: -  delta - initial value for diagonal scaling

628:    Level: intermediate
629: @*/

631: PetscErrorCode MatSymBrdnSetDelta(Mat B, PetscScalar delta)
632: {
633:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
634:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
635:   PetscErrorCode    ierr;
636:   PetscBool         is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn;
637: 
639:   PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);
640:   PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);
641:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBRDN, &is_symbrdn);
642:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBRDN, &is_symbadbrdn);
643:   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");
644:   lsb->delta = PetscAbsReal(PetscRealPart(delta));
645:   lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
646:   lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
647:   return(0);
648: }

650: /*------------------------------------------------------------*/

652: /*@
653:    MatCreateLMVMSymBrdn - Creates a limited-memory Symmetric Broyden-type matrix used 
654:    for approximating Jacobians. L-SymBrdn is a convex combination of L-DFP and 
655:    L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor 
656:    phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed 
657:    to be symmetric positive-definite.
658:    
659:    The provided local and global sizes must match the solution and function vectors 
660:    used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have 
661:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in 
662:    parallel. To use the L-SymBrdn matrix with other vector types, the matrix must be 
663:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate(). 
664:    This ensures that the internal storage and work vectors are duplicated from the 
665:    correct type of vector.

667:    Collective

669:    Input Parameters:
670: +  comm - MPI communicator, set to PETSC_COMM_SELF
671: .  n - number of local rows for storage vectors
672: -  N - global size of the storage vectors

674:    Output Parameter:
675: .  B - the matrix

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

680:    Options Database Keys:
681: +   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
682: .   -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
683: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
684: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
685: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
686: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
687: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
688: -   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

690:    Level: intermediate

692: .seealso: MatCreate(), MATLMVM, MATLMVMSYMBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(), 
693:           MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn()
694: @*/
695: PetscErrorCode MatCreateLMVMSymBrdn(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
696: {
697:   PetscErrorCode    ierr;
698: 
700:   MatCreate(comm, B);
701:   MatSetSizes(*B, n, n, N, N);
702:   MatSetType(*B, MATLMVMSYMBRDN);
703:   MatSetUp(*B);
704:   return(0);
705: }

707: /*------------------------------------------------------------*/

709: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
710: {
711:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
712:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
713:   PetscErrorCode    ierr;
714: 
716:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
717:     MatLMVMApplyJ0Fwd(B, X, Z);
718:   } else {
719:     switch (lsb->scale_type) {
720:     case SYMBRDN_SCALE_SCALAR:
721:       VecCopy(X, Z);
722:       VecScale(Z, 1.0/lsb->sigma);
723:       break;
724:     case SYMBRDN_SCALE_DIAG:
725:       MatMult(lsb->D, X, Z);
726:       break;
727:     case SYMBRDN_SCALE_NONE:
728:     default:
729:       VecCopy(X, Z);
730:       break;
731:     }
732:   }
733:   return(0);
734: }

736: /*------------------------------------------------------------*/

738: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
739: {
740:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
741:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
742:   PetscErrorCode    ierr;
743: 
745:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
746:     MatLMVMApplyJ0Inv(B, F, dX);
747:   } else {
748:     switch (lsb->scale_type) {
749:     case SYMBRDN_SCALE_SCALAR:
750:       VecCopy(F, dX);
751:       VecScale(dX, lsb->sigma);
752:       break;
753:     case SYMBRDN_SCALE_DIAG:
754:       MatSolve(lsb->D, F, dX);
755:       break;
756:     case SYMBRDN_SCALE_NONE:
757:     default:
758:       VecCopy(F, dX);
759:       break;
760:     }
761:   }
762:   return(0);
763: }

765: /*------------------------------------------------------------*/

767: PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
768: {
769:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
770:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
771:   PetscInt          i, start;
772:   PetscReal         a, b, c, sig1, sig2, signew;
773: 
775:   if (lsb->sigma_hist == 0) {
776:     signew = 1.0;
777:   } else {
778:     start = PetscMax(0, lmvm->k-lsb->sigma_hist+1);
779:     signew = 0.0;
780:     if (lsb->alpha == 1.0) {
781:       for (i = start; i <= lmvm->k; ++i) {
782:         signew += lsb->yts[i]/lsb->yty[i];
783:       }
784:     } else if (lsb->alpha == 0.5) {
785:       for (i = start; i <= lmvm->k; ++i) {
786:         signew += lsb->sts[i]/lsb->yty[i];
787:       }
788:       signew = PetscSqrtReal(signew);
789:     } else if (lsb->alpha == 0.0) {
790:       for (i = start; i <= lmvm->k; ++i) {
791:         signew += lsb->sts[i]/lsb->yts[i];
792:       }
793:     } else {
794:       /* compute coefficients of the quadratic */
795:       a = b = c = 0.0;
796:       for (i = start; i <= lmvm->k; ++i) {
797:         a += lsb->yty[i];
798:         b += lsb->yts[i];
799:         c += lsb->sts[i];
800:       }
801:       a *= lsb->alpha;
802:       b *= -(2.0*lsb->alpha - 1.0);
803:       c *= lsb->alpha - 1.0;
804:       /* use quadratic formula to find roots */
805:       sig1 = (-b + PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
806:       sig2 = (-b - PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
807:       /* accept the positive root as the scalar */
808:       if (sig1 > 0.0) {
809:         signew = sig1;
810:       } else if (sig2 > 0.0) {
811:         signew = sig2;
812:       } else {
813:         SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
814:       }
815:     }
816:   }
817:   lsb->sigma = lsb->rho*signew + (1.0 - lsb->rho)*lsb->sigma;
818:   return(0);
819: }