Actual source code: symbrdn.c

petsc-master 2020-08-10
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: 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).
12:
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.
16:
17:   dX <- J0^{-1} * F
18:
19:   for i=0,1,2,...,k
20:     # Q[i] = (B_i)^T{-1} Y[i]
21:
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)
26:
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;
40:
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:   }
51:
52:   VecCheckSameSize(F, 2, dX, 3);
53:   VecCheckMatCompatible(B, dX, 3, F, 2);
54:
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 Section 1.5 Writing Application Codes with PETSc*/
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 Section 1.5 Writing Application Codes with PETSc*/
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:   }
111:
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).
137:
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.
142:
143:   Z <- J0 * X
144:
145:   for i=0,1,2,...,k
146:     # P[i] = (B_k) * S[i]
147:
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)
152:
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;
165:
166:
168:   /* Efficient shortcuts for pure BFGS and pure DFP configurations */
169:   if (lsb->phi == 0.0) {
170:     MatMult_LMVMBFGS(B, X, Z);
171:     return(0);
172:   }
173:   if (lsb->phi == 1.0) {
174:     MatMult_LMVMDFP(B, X, Z);
175:     return(0);
176:   }
177:
178:   VecCheckSameSize(X, 2, Z, 3);
179:   VecCheckMatCompatible(B, X, 2, Z, 3);
180:
181:   if (lsb->needP) {
182:     /* Start the loop for (P[k] = (B_k) * S[k]) */
183:     for (i = 0; i <= lmvm->k; ++i) {
184:       MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);
185:       for (j = 0; j <= i-1; ++j) {
186:         /* Compute the necessary dot products */
187:         VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);
188:         VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
189:         VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);
190:         VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
191:         /* Compute the pure BFGS component of the forward product */
192:         VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi)/lsb->stp[j], PetscRealPart(yjtsi)/lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);
193:         /* Tack on the convexly scaled extras to the forward product */
194:         if (lsb->phi > 0.0) {
195:           VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);
196:           VecDot(lsb->work, lmvm->S[i], &wtsi);
197:           VecAXPY(lsb->P[i], lsb->phi*lsb->stp[j]*PetscRealPart(wtsi), lsb->work);
198:         }
199:       }
200:       VecDot(lmvm->S[i], lsb->P[i], &stp);
201:       lsb->stp[i] = PetscRealPart(stp);
202:     }
203:     lsb->needP = PETSC_FALSE;
204:   }
205:
206:   /* Start the outer iterations for (B * X) */
207:   MatSymBrdnApplyJ0Fwd(B, X, Z);
208:   for (i = 0; i <= lmvm->k; ++i) {
209:     /* Compute the necessary dot products */
210:     VecDotBegin(lmvm->S[i], Z, &stz);
211:     VecDotBegin(lmvm->Y[i], X, &ytx);
212:     VecDotEnd(lmvm->S[i], Z, &stz);
213:     VecDotEnd(lmvm->Y[i], X, &ytx);
214:     /* Compute the pure BFGS component */
215:     VecAXPBYPCZ(Z, -PetscRealPart(stz)/lsb->stp[i], PetscRealPart(ytx)/lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]);
216:     /* Tack on the convexly scaled extras */
217:     VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]);
218:     VecDot(lsb->work, X, &wtx);
219:     VecAXPY(Z, lsb->phi*lsb->stp[i]*PetscRealPart(wtx), lsb->work);
220:   }
221:   return(0);
222: }

224: /*------------------------------------------------------------*/

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

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

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

318: /*------------------------------------------------------------*/

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

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

365: /*------------------------------------------------------------*/

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

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

417: /*------------------------------------------------------------*/

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

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

447: /*------------------------------------------------------------*/

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

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

470: /*------------------------------------------------------------*/

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

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

504: /*------------------------------------------------------------*/

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

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

528: /*------------------------------------------------------------*/

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

537:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
538:   PetscOptionsHead(PetscOptionsObject,"Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
539:   PetscOptionsReal("-mat_lmvm_phi","(developer) convex ratio between BFGS and DFP components of the update","",lsb->phi,&lsb->phi,NULL);
540:   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]");
541:   MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
542:   PetscOptionsTail();
543:   return(0);
544: }

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

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

582: /*------------------------------------------------------------*/

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

591:   MatCreate_LMVM(B);
592:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN);
593:   MatSetOption(B, MAT_SPD, PETSC_TRUE);
594:   B->ops->view = MatView_LMVMSymBrdn;
595:   B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
596:   B->ops->setup = MatSetUp_LMVMSymBrdn;
597:   B->ops->destroy = MatDestroy_LMVMSymBrdn;
598:   B->ops->solve = MatSolve_LMVMSymBrdn;
599:
600:   lmvm = (Mat_LMVM*)B->data;
601:   lmvm->square = PETSC_TRUE;
602:   lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
603:   lmvm->ops->reset = MatReset_LMVMSymBrdn;
604:   lmvm->ops->update = MatUpdate_LMVMSymBrdn;
605:   lmvm->ops->mult = MatMult_LMVMSymBrdn;
606:   lmvm->ops->copy = MatCopy_LMVMSymBrdn;
607:
608:   PetscNewLog(B, &lsb);
609:   lmvm->ctx = (void*)lsb;
610:   lsb->allocated       = PETSC_FALSE;
611:   lsb->needP           = lsb->needQ = PETSC_TRUE;
612:   lsb->phi             = 0.125;
613:   lsb->theta           = 0.125;
614:   lsb->alpha           = 1.0;
615:   lsb->rho             = 1.0;
616:   lsb->beta            = 0.5;
617:   lsb->sigma           = 1.0;
618:   lsb->delta           = 1.0;
619:   lsb->delta_min       = 1e-7;
620:   lsb->delta_max       = 100.0;
621:   lsb->sigma_hist      = 1;
622:   lsb->scale_type      = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
623:   lsb->watchdog        = 0;
624:   lsb->max_seq_rejects = lmvm->m/2;
625:
626:   MatCreate(PetscObjectComm((PetscObject)B), &lsb->D);
627:   MatSetType(lsb->D, MATLMVMDIAGBROYDEN);
628:   MatSetOptionsPrefix(lsb->D, "J0_");
629:   return(0);
630: }

632: /*------------------------------------------------------------*/

634: /*@
635:    MatLMVMSymBroydenSetDelta - Sets the starting value for the diagonal scaling vector computed
636:    in the SymBrdn approximations (also works for BFGS and DFP).
637:
638:    Input Parameters:
639: +  B - LMVM matrix
640: -  delta - initial value for diagonal scaling

642:    Level: intermediate
643: @*/

645: PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
646: {
647:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
648:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
649:   PetscErrorCode    ierr;
650:   PetscBool         is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn;
651:
653:   PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);
654:   PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);
655:   PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn);
657:   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");
658:   lsb->delta = PetscAbsReal(PetscRealPart(delta));
659:   lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
660:   lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
661:   return(0);
662: }

664: /*------------------------------------------------------------*/

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

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

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

676:     Level: intermediate

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

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

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

696: /*------------------------------------------------------------*/

698: /*@
699:    MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
700:    for approximating Jacobians. L-SymBrdn is a convex combination of L-DFP and
701:    L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor
702:    phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed
703:    to be symmetric positive-definite.
704:
705:    The provided local and global sizes must match the solution and function vectors
706:    used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have
707:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
708:    parallel. To use the L-SymBrdn matrix with other vector types, the matrix must be
709:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
710:    This ensures that the internal storage and work vectors are duplicated from the
711:    correct type of vector.

713:    Collective

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

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

723:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()

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

736:    Level: intermediate

738: .seealso: MatCreate(), MATLMVM, MATLMVMSYMBROYDEN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
740: @*/
741: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
742: {
743:   PetscErrorCode    ierr;
744:
746:   MatCreate(comm, B);
747:   MatSetSizes(*B, n, n, N, N);
748:   MatSetType(*B, MATLMVMSYMBROYDEN);
749:   MatSetUp(*B);
750:   return(0);
751: }

753: /*------------------------------------------------------------*/

755: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
756: {
757:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
758:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
759:   PetscErrorCode    ierr;
760:
762:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
763:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
764:     MatLMVMApplyJ0Fwd(B, X, Z);
765:   } else {
766:     switch (lsb->scale_type) {
767:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
768:       VecCopy(X, Z);
769:       VecScale(Z, 1.0/lsb->sigma);
770:       break;
771:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
772:       MatMult(lsb->D, X, Z);
773:       break;
774:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
775:     default:
776:       VecCopy(X, Z);
777:       break;
778:     }
779:   }
780:   return(0);
781: }

783: /*------------------------------------------------------------*/

785: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
786: {
787:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
788:   Mat_SymBrdn       *lsb = (Mat_SymBrdn*)lmvm->ctx;
789:   PetscErrorCode    ierr;
790:
792:   if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
793:     lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
794:     MatLMVMApplyJ0Inv(B, F, dX);
795:   } else {
796:     switch (lsb->scale_type) {
797:     case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
798:       VecCopy(F, dX);
799:       VecScale(dX, lsb->sigma);
800:       break;
801:     case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
802:       MatSolve(lsb->D, F, dX);
803:       break;
804:     case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
805:     default:
806:       VecCopy(F, dX);
807:       break;
808:     }
809:   }
810:   return(0);
811: }

813: /*------------------------------------------------------------*/

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