Actual source code: symbrdn.c
petsc-3.14.6 2021-03-30
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;
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: }
178: VecCheckSameSize(X, 2, Z, 3);
179: VecCheckMatCompatible(B, X, 2, Z, 3);
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: }
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: }
299: /* Update the scaling */
300: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
301: MatLMVMUpdate(lsb->D, X, F);
302: }
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;
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;
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;
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).
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;
653: PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);
654: PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);
655: PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn);
656: PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn);
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.
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()
724: paradigm instead of this routine directly.
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(),
739: MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn()
740: @*/
741: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
742: {
743: PetscErrorCode ierr;
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;
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;
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;
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: }