Actual source code: symbrdn.c
1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
4: const char *const MatLMVMSymBroydenScaleTypes[] = {"NONE","SCALAR","DIAGONAL","USER","MatLMVMSymBrdnScaleType","MAT_LMVM_SYMBROYDEN_SCALING_",NULL};
6: /*------------------------------------------------------------*/
8: /*
9: The solution method below is the matrix-free implementation of
10: Equation 8.6a in Dennis and More "Quasi-Newton Methods, Motivation
11: and Theory" (https://epubs.siam.org/doi/abs/10.1137/1019005).
13: Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
14: the matrix is updated with a new (S[i], Y[i]) pair. This allows
15: repeated calls of MatSolve without incurring redundant computation.
17: dX <- J0^{-1} * F
19: for i=0,1,2,...,k
20: # Q[i] = (B_i)^T{-1} Y[i]
22: rho = 1.0 / (Y[i]^T S[i])
23: alpha = rho * (S[i]^T F)
24: zeta = 1.0 / (Y[i]^T Q[i])
25: gamma = zeta * (Y[i]^T dX)
27: dX <- dX - (gamma * Q[i]) + (alpha * Y[i])
28: W <- (rho * S[i]) - (zeta * Q[i])
29: dX <- dX + (psi[i] * (Y[i]^T Q[i]) * (W^T F) * W)
30: end
31: */
32: static PetscErrorCode MatSolve_LMVMSymBrdn(Mat B, Vec F, Vec dX)
33: {
34: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
35: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
36: PetscErrorCode ierr;
37: PetscInt i, j;
38: PetscReal numer;
39: PetscScalar sjtpi, yjtsi, wtsi, yjtqi, sjtyi, wtyi, ytx, stf, wtf, stp, ytq;
42: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
43: if (lsb->phi == 0.0) {
44: MatSolve_LMVMBFGS(B, F, dX);
45: return(0);
46: }
47: if (lsb->phi == 1.0) {
48: MatSolve_LMVMDFP(B, F, dX);
49: return(0);
50: }
52: VecCheckSameSize(F, 2, dX, 3);
53: VecCheckMatCompatible(B, dX, 3, F, 2);
55: if (lsb->needP) {
56: /* Start the loop for (P[k] = (B_k) * S[k]) */
57: for (i = 0; i <= lmvm->k; ++i) {
58: MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);
59: for (j = 0; j <= i-1; ++j) {
60: /* Compute the necessary dot products */
61: VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);
62: VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
63: VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);
64: VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
65: /* Compute the pure BFGS component of the forward product */
66: VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi)/lsb->stp[j], PetscRealPart(yjtsi)/lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);
67: /* Tack on the convexly scaled extras to the forward product */
68: if (lsb->phi > 0.0) {
69: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);
70: VecDot(lsb->work, lmvm->S[i], &wtsi);
71: VecAXPY(lsb->P[i], lsb->phi*lsb->stp[j]*PetscRealPart(wtsi), lsb->work);
72: }
73: }
74: VecDot(lmvm->S[i], lsb->P[i], &stp);
75: lsb->stp[i] = PetscRealPart(stp);
76: }
77: lsb->needP = PETSC_FALSE;
78: }
79: if (lsb->needQ) {
80: /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
81: for (i = 0; i <= lmvm->k; ++i) {
82: MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]);
83: for (j = 0; j <= i-1; ++j) {
84: /* Compute the necessary dot products */
85: VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi);
86: VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
87: VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi);
88: VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
89: /* Compute the pure DFP component of the inverse application*/
90: VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi)/lsb->ytq[j], PetscRealPart(sjtyi)/lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]);
91: /* Tack on the convexly scaled extras to the inverse application*/
92: if (lsb->psi[j] > 0.0) {
93: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]);
94: VecDot(lsb->work, lmvm->Y[i], &wtyi);
95: VecAXPY(lsb->Q[i], lsb->psi[j]*lsb->ytq[j]*PetscRealPart(wtyi), lsb->work);
96: }
97: }
98: VecDot(lmvm->Y[i], lsb->Q[i], &ytq);
99: lsb->ytq[i] = PetscRealPart(ytq);
100: if (lsb->phi == 1.0) {
101: lsb->psi[i] = 0.0;
102: } else if (lsb->phi == 0.0) {
103: lsb->psi[i] = 1.0;
104: } else {
105: numer = (1.0 - lsb->phi)*lsb->yts[i]*lsb->yts[i];
106: lsb->psi[i] = numer / (numer + (lsb->phi*lsb->ytq[i]*lsb->stp[i]));
107: }
108: }
109: lsb->needQ = PETSC_FALSE;
110: }
112: /* Start the outer iterations for ((B^{-1}) * dX) */
113: MatSymBrdnApplyJ0Inv(B, F, dX);
114: for (i = 0; i <= lmvm->k; ++i) {
115: /* Compute the necessary dot products -- store yTs and yTp for inner iterations later */
116: VecDotBegin(lmvm->Y[i], dX, &ytx);
117: VecDotBegin(lmvm->S[i], F, &stf);
118: VecDotEnd(lmvm->Y[i], dX, &ytx);
119: VecDotEnd(lmvm->S[i], F, &stf);
120: /* Compute the pure DFP component */
121: VecAXPBYPCZ(dX, -PetscRealPart(ytx)/lsb->ytq[i], PetscRealPart(stf)/lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]);
122: /* Tack on the convexly scaled extras */
123: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]);
124: VecDot(lsb->work, F, &wtf);
125: VecAXPY(dX, lsb->psi[i]*lsb->ytq[i]*PetscRealPart(wtf), lsb->work);
126: }
128: return(0);
129: }
131: /*------------------------------------------------------------*/
133: /*
134: The forward-product below is the matrix-free implementation of
135: Equation 16 in Dennis and Wolkowicz "Sizing and Least Change Secant
136: Methods" (http://www.caam.rice.edu/caam/trs/90/TR90-05.pdf).
138: P[i] = (B_i)*S[i] terms are computed ahead of time whenever
139: the matrix is updated with a new (S[i], Y[i]) pair. This allows
140: repeated calls of MatMult inside KSP solvers without unnecessarily
141: recomputing P[i] terms in expensive nested-loops.
143: Z <- J0 * X
145: for i=0,1,2,...,k
146: # P[i] = (B_k) * S[i]
148: rho = 1.0 / (Y[i]^T S[i])
149: alpha = rho * (Y[i]^T F)
150: zeta = 1.0 / (S[i]^T P[i])
151: gamma = zeta * (S[i]^T dX)
153: dX <- dX - (gamma * P[i]) + (alpha * S[i])
154: W <- (rho * Y[i]) - (zeta * P[i])
155: dX <- dX + (phi * (S[i]^T P[i]) * (W^T F) * W)
156: end
157: */
158: static PetscErrorCode MatMult_LMVMSymBrdn(Mat B, Vec X, Vec Z)
159: {
160: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
161: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
162: PetscErrorCode ierr;
163: PetscInt i, j;
164: PetscScalar sjtpi, yjtsi, wtsi, stz, ytx, wtx, stp;
167: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
168: if (lsb->phi == 0.0) {
169: MatMult_LMVMBFGS(B, X, Z);
170: return(0);
171: }
172: if (lsb->phi == 1.0) {
173: MatMult_LMVMDFP(B, X, Z);
174: return(0);
175: }
177: VecCheckSameSize(X, 2, Z, 3);
178: VecCheckMatCompatible(B, X, 2, Z, 3);
180: if (lsb->needP) {
181: /* Start the loop for (P[k] = (B_k) * S[k]) */
182: for (i = 0; i <= lmvm->k; ++i) {
183: MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);
184: for (j = 0; j <= i-1; ++j) {
185: /* Compute the necessary dot products */
186: VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);
187: VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
188: VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);
189: VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
190: /* Compute the pure BFGS component of the forward product */
191: VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi)/lsb->stp[j], PetscRealPart(yjtsi)/lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);
192: /* Tack on the convexly scaled extras to the forward product */
193: if (lsb->phi > 0.0) {
194: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);
195: VecDot(lsb->work, lmvm->S[i], &wtsi);
196: VecAXPY(lsb->P[i], lsb->phi*lsb->stp[j]*PetscRealPart(wtsi), lsb->work);
197: }
198: }
199: VecDot(lmvm->S[i], lsb->P[i], &stp);
200: lsb->stp[i] = PetscRealPart(stp);
201: }
202: lsb->needP = PETSC_FALSE;
203: }
205: /* Start the outer iterations for (B * X) */
206: MatSymBrdnApplyJ0Fwd(B, X, Z);
207: for (i = 0; i <= lmvm->k; ++i) {
208: /* Compute the necessary dot products */
209: VecDotBegin(lmvm->S[i], Z, &stz);
210: VecDotBegin(lmvm->Y[i], X, &ytx);
211: VecDotEnd(lmvm->S[i], Z, &stz);
212: VecDotEnd(lmvm->Y[i], X, &ytx);
213: /* Compute the pure BFGS component */
214: VecAXPBYPCZ(Z, -PetscRealPart(stz)/lsb->stp[i], PetscRealPart(ytx)/lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]);
215: /* Tack on the convexly scaled extras */
216: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]);
217: VecDot(lsb->work, X, &wtx);
218: VecAXPY(Z, lsb->phi*lsb->stp[i]*PetscRealPart(wtx), lsb->work);
219: }
220: return(0);
221: }
223: /*------------------------------------------------------------*/
225: static PetscErrorCode MatUpdate_LMVMSymBrdn(Mat B, Vec X, Vec F)
226: {
227: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
228: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
229: Mat_LMVM *dbase;
230: Mat_DiagBrdn *dctx;
231: PetscErrorCode ierr;
232: PetscInt old_k, i;
233: PetscReal curvtol;
234: PetscScalar curvature, ytytmp, ststmp;
237: if (!lmvm->m) return(0);
238: if (lmvm->prev_set) {
239: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
240: VecAYPX(lmvm->Xprev, -1.0, X);
241: VecAYPX(lmvm->Fprev, -1.0, F);
242: /* Test if the updates can be accepted */
243: VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
244: VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
245: VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
246: VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
247: if (PetscRealPart(ststmp) < lmvm->eps) {
248: curvtol = 0.0;
249: } else {
250: curvtol = lmvm->eps * PetscRealPart(ststmp);
251: }
252: if (PetscRealPart(curvature) > curvtol) {
253: /* Update is good, accept it */
254: lsb->watchdog = 0;
255: lsb->needP = lsb->needQ = PETSC_TRUE;
256: old_k = lmvm->k;
257: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
258: /* If we hit the memory limit, shift the yts, yty and sts arrays */
259: if (old_k == lmvm->k) {
260: for (i = 0; i <= lmvm->k-1; ++i) {
261: lsb->yts[i] = lsb->yts[i+1];
262: lsb->yty[i] = lsb->yty[i+1];
263: lsb->sts[i] = lsb->sts[i+1];
264: }
265: }
266: /* Update history of useful scalars */
267: VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &ytytmp);
268: lsb->yts[lmvm->k] = PetscRealPart(curvature);
269: lsb->yty[lmvm->k] = PetscRealPart(ytytmp);
270: lsb->sts[lmvm->k] = PetscRealPart(ststmp);
271: /* Compute the scalar scale if necessary */
272: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) {
273: MatSymBrdnComputeJ0Scalar(B);
274: }
275: } else {
276: /* Update is bad, skip it */
277: ++lmvm->nrejects;
278: ++lsb->watchdog;
279: }
280: } else {
281: switch (lsb->scale_type) {
282: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
283: dbase = (Mat_LMVM*)lsb->D->data;
284: dctx = (Mat_DiagBrdn*)dbase->ctx;
285: VecSet(dctx->invD, lsb->delta);
286: break;
287: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
288: lsb->sigma = lsb->delta;
289: break;
290: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
291: lsb->sigma = 1.0;
292: break;
293: default:
294: break;
295: }
296: }
298: /* Update the scaling */
299: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
300: MatLMVMUpdate(lsb->D, X, F);
301: }
303: if (lsb->watchdog > lsb->max_seq_rejects) {
304: MatLMVMReset(B, PETSC_FALSE);
305: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
306: MatLMVMReset(lsb->D, PETSC_FALSE);
307: }
308: }
310: /* Save the solution and function to be used in the next update */
311: VecCopy(X, lmvm->Xprev);
312: VecCopy(F, lmvm->Fprev);
313: lmvm->prev_set = PETSC_TRUE;
314: return(0);
315: }
317: /*------------------------------------------------------------*/
319: static PetscErrorCode MatCopy_LMVMSymBrdn(Mat B, Mat M, MatStructure str)
320: {
321: Mat_LMVM *bdata = (Mat_LMVM*)B->data;
322: Mat_SymBrdn *blsb = (Mat_SymBrdn*)bdata->ctx;
323: Mat_LMVM *mdata = (Mat_LMVM*)M->data;
324: Mat_SymBrdn *mlsb = (Mat_SymBrdn*)mdata->ctx;
325: PetscErrorCode ierr;
326: PetscInt i;
329: mlsb->phi = blsb->phi;
330: mlsb->needP = blsb->needP;
331: mlsb->needQ = blsb->needQ;
332: for (i=0; i<=bdata->k; ++i) {
333: mlsb->stp[i] = blsb->stp[i];
334: mlsb->ytq[i] = blsb->ytq[i];
335: mlsb->yts[i] = blsb->yts[i];
336: mlsb->psi[i] = blsb->psi[i];
337: VecCopy(blsb->P[i], mlsb->P[i]);
338: VecCopy(blsb->Q[i], mlsb->Q[i]);
339: }
340: mlsb->scale_type = blsb->scale_type;
341: mlsb->alpha = blsb->alpha;
342: mlsb->beta = blsb->beta;
343: mlsb->rho = blsb->rho;
344: mlsb->delta = blsb->delta;
345: mlsb->sigma_hist = blsb->sigma_hist;
346: mlsb->watchdog = blsb->watchdog;
347: mlsb->max_seq_rejects = blsb->max_seq_rejects;
348: switch (blsb->scale_type) {
349: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
350: mlsb->sigma = blsb->sigma;
351: break;
352: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
353: MatCopy(blsb->D, mlsb->D, SAME_NONZERO_PATTERN);
354: break;
355: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
356: mlsb->sigma = 1.0;
357: break;
358: default:
359: break;
360: }
361: return(0);
362: }
364: /*------------------------------------------------------------*/
366: static PetscErrorCode MatReset_LMVMSymBrdn(Mat B, PetscBool destructive)
367: {
368: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
369: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
370: Mat_LMVM *dbase;
371: Mat_DiagBrdn *dctx;
372: PetscErrorCode ierr;
375: lsb->watchdog = 0;
376: lsb->needP = lsb->needQ = PETSC_TRUE;
377: if (lsb->allocated) {
378: if (destructive) {
379: VecDestroy(&lsb->work);
380: PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);
381: PetscFree(lsb->psi);
382: VecDestroyVecs(lmvm->m, &lsb->P);
383: VecDestroyVecs(lmvm->m, &lsb->Q);
384: switch (lsb->scale_type) {
385: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
386: MatLMVMReset(lsb->D, PETSC_TRUE);
387: break;
388: default:
389: break;
390: }
391: lsb->allocated = PETSC_FALSE;
392: } else {
393: PetscMemzero(lsb->psi, lmvm->m);
394: switch (lsb->scale_type) {
395: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
396: lsb->sigma = lsb->delta;
397: break;
398: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
399: MatLMVMReset(lsb->D, PETSC_FALSE);
400: dbase = (Mat_LMVM*)lsb->D->data;
401: dctx = (Mat_DiagBrdn*)dbase->ctx;
402: VecSet(dctx->invD, lsb->delta);
403: break;
404: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
405: lsb->sigma = 1.0;
406: break;
407: default:
408: break;
409: }
410: }
411: }
412: MatReset_LMVM(B, destructive);
413: return(0);
414: }
416: /*------------------------------------------------------------*/
418: static PetscErrorCode MatAllocate_LMVMSymBrdn(Mat B, Vec X, Vec F)
419: {
420: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
421: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
422: PetscErrorCode ierr;
425: MatAllocate_LMVM(B, X, F);
426: if (!lsb->allocated) {
427: VecDuplicate(X, &lsb->work);
428: if (lmvm->m > 0) {
429: PetscMalloc5(lmvm->m,&lsb->stp,lmvm->m,&lsb->ytq,lmvm->m,&lsb->yts,lmvm->m,&lsb->yty,lmvm->m,&lsb->sts);
430: PetscCalloc1(lmvm->m,&lsb->psi);
431: VecDuplicateVecs(X, lmvm->m, &lsb->P);
432: VecDuplicateVecs(X, lmvm->m, &lsb->Q);
433: }
434: switch (lsb->scale_type) {
435: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
436: MatLMVMAllocate(lsb->D, X, F);
437: break;
438: default:
439: break;
440: }
441: lsb->allocated = PETSC_TRUE;
442: }
443: return(0);
444: }
446: /*------------------------------------------------------------*/
448: static PetscErrorCode MatDestroy_LMVMSymBrdn(Mat B)
449: {
450: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
451: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
452: PetscErrorCode ierr;
455: if (lsb->allocated) {
456: VecDestroy(&lsb->work);
457: PetscFree5(lsb->stp, lsb->ytq, lsb->yts, lsb->yty, lsb->sts);
458: PetscFree(lsb->psi);
459: VecDestroyVecs(lmvm->m, &lsb->P);
460: VecDestroyVecs(lmvm->m, &lsb->Q);
461: lsb->allocated = PETSC_FALSE;
462: }
463: MatDestroy(&lsb->D);
464: PetscFree(lmvm->ctx);
465: MatDestroy_LMVM(B);
466: return(0);
467: }
469: /*------------------------------------------------------------*/
471: static PetscErrorCode MatSetUp_LMVMSymBrdn(Mat B)
472: {
473: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
474: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
475: PetscErrorCode ierr;
476: PetscInt n, N;
479: MatSetUp_LMVM(B);
480: if (!lsb->allocated) {
481: VecDuplicate(lmvm->Xprev, &lsb->work);
482: if (lmvm->m > 0) {
483: PetscMalloc5(lmvm->m,&lsb->stp,lmvm->m,&lsb->ytq,lmvm->m,&lsb->yts,lmvm->m,&lsb->yty,lmvm->m,&lsb->sts);
484: PetscCalloc1(lmvm->m,&lsb->psi);
485: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->P);
486: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lsb->Q);
487: }
488: switch (lsb->scale_type) {
489: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
490: MatGetLocalSize(B, &n, &n);
491: MatGetSize(B, &N, &N);
492: MatSetSizes(lsb->D, n, n, N, N);
493: MatSetUp(lsb->D);
494: break;
495: default:
496: break;
497: }
498: lsb->allocated = PETSC_TRUE;
499: }
500: return(0);
501: }
503: /*------------------------------------------------------------*/
505: PetscErrorCode MatView_LMVMSymBrdn(Mat B, PetscViewer pv)
506: {
507: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
508: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
509: PetscErrorCode ierr;
510: PetscBool isascii;
513: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
514: if (isascii) {
515: PetscViewerASCIIPrintf(pv,"Scale type: %s\n",MatLMVMSymBroydenScaleTypes[lsb->scale_type]);
516: PetscViewerASCIIPrintf(pv,"Scale history: %d\n",lsb->sigma_hist);
517: PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)lsb->alpha, (double)lsb->beta, (double)lsb->rho);
518: PetscViewerASCIIPrintf(pv,"Convex factors: phi=%g, theta=%g\n",(double)lsb->phi, (double)lsb->theta);
519: }
520: MatView_LMVM(B, pv);
521: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
522: MatView(lsb->D, pv);
523: }
524: return(0);
525: }
527: /*------------------------------------------------------------*/
529: PetscErrorCode MatSetFromOptions_LMVMSymBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
530: {
531: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
532: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
533: PetscErrorCode ierr;
536: MatSetFromOptions_LMVM(PetscOptionsObject, B);
537: PetscOptionsHead(PetscOptionsObject,"Restricted/Symmetric Broyden method for approximating SPD Jacobian actions (MATLMVMSYMBRDN)");
538: PetscOptionsReal("-mat_lmvm_phi","(developer) convex ratio between BFGS and DFP components of the update","",lsb->phi,&lsb->phi,NULL);
539: if ((lsb->phi < 0.0) || (lsb->phi > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the update formula cannot be outside the range of [0, 1]");
540: MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
541: PetscOptionsTail();
542: return(0);
543: }
545: PetscErrorCode MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionItems *PetscOptionsObject, Mat B)
546: {
547: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
548: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
549: Mat_LMVM *dbase;
550: Mat_DiagBrdn *dctx;
551: MatLMVMSymBroydenScaleType stype = lsb->scale_type;
552: PetscBool flg;
553: PetscErrorCode ierr;
556: PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",lsb->beta,&lsb->beta,NULL);
557: PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",lsb->theta,&lsb->theta,NULL);
558: if ((lsb->theta < 0.0) || (lsb->theta > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
559: PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",lsb->rho,&lsb->rho,NULL);
560: if ((lsb->rho < 0.0) || (lsb->rho > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "update limiter in the J0 scaling cannot be outside the range of [0, 1]");
561: PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",lsb->alpha,&lsb->alpha,NULL);
562: if ((lsb->alpha < 0.0) || (lsb->alpha > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
563: PetscOptionsBoundedInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",lsb->sigma_hist,&lsb->sigma_hist,NULL,1);
564: PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0","MatLMVMSymBrdnScaleType",MatLMVMSymBroydenScaleTypes,(PetscEnum)stype,(PetscEnum*)&stype,&flg);
565: if (flg) {MatLMVMSymBroydenSetScaleType(B, stype);}
566: if (lsb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
567: MatSetFromOptions(lsb->D);
568: dbase = (Mat_LMVM*)lsb->D->data;
569: dctx = (Mat_DiagBrdn*)dbase->ctx;
570: dctx->delta_min = lsb->delta_min;
571: dctx->delta_max = lsb->delta_max;
572: dctx->theta = lsb->theta;
573: dctx->rho = lsb->rho;
574: dctx->alpha = lsb->alpha;
575: dctx->beta = lsb->beta;
576: dctx->sigma_hist = lsb->sigma_hist;
577: }
578: return(0);
579: }
581: /*------------------------------------------------------------*/
583: PetscErrorCode MatCreate_LMVMSymBrdn(Mat B)
584: {
585: Mat_LMVM *lmvm;
586: Mat_SymBrdn *lsb;
587: PetscErrorCode ierr;
590: MatCreate_LMVM(B);
591: PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBROYDEN);
592: MatSetOption(B, MAT_SPD, PETSC_TRUE);
593: B->ops->view = MatView_LMVMSymBrdn;
594: B->ops->setfromoptions = MatSetFromOptions_LMVMSymBrdn;
595: B->ops->setup = MatSetUp_LMVMSymBrdn;
596: B->ops->destroy = MatDestroy_LMVMSymBrdn;
597: B->ops->solve = MatSolve_LMVMSymBrdn;
599: lmvm = (Mat_LMVM*)B->data;
600: lmvm->square = PETSC_TRUE;
601: lmvm->ops->allocate = MatAllocate_LMVMSymBrdn;
602: lmvm->ops->reset = MatReset_LMVMSymBrdn;
603: lmvm->ops->update = MatUpdate_LMVMSymBrdn;
604: lmvm->ops->mult = MatMult_LMVMSymBrdn;
605: lmvm->ops->copy = MatCopy_LMVMSymBrdn;
607: PetscNewLog(B, &lsb);
608: lmvm->ctx = (void*)lsb;
609: lsb->allocated = PETSC_FALSE;
610: lsb->needP = lsb->needQ = PETSC_TRUE;
611: lsb->phi = 0.125;
612: lsb->theta = 0.125;
613: lsb->alpha = 1.0;
614: lsb->rho = 1.0;
615: lsb->beta = 0.5;
616: lsb->sigma = 1.0;
617: lsb->delta = 1.0;
618: lsb->delta_min = 1e-7;
619: lsb->delta_max = 100.0;
620: lsb->sigma_hist = 1;
621: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL;
622: lsb->watchdog = 0;
623: lsb->max_seq_rejects = lmvm->m/2;
625: MatCreate(PetscObjectComm((PetscObject)B), &lsb->D);
626: MatSetType(lsb->D, MATLMVMDIAGBROYDEN);
627: MatSetOptionsPrefix(lsb->D, "J0_");
628: return(0);
629: }
631: /*------------------------------------------------------------*/
633: /*@
634: MatLMVMSymBroydenSetDelta - Sets the starting value for the diagonal scaling vector computed
635: in the SymBrdn approximations (also works for BFGS and DFP).
637: Input Parameters:
638: + B - LMVM matrix
639: - delta - initial value for diagonal scaling
641: Level: intermediate
642: @*/
644: PetscErrorCode MatLMVMSymBroydenSetDelta(Mat B, PetscScalar delta)
645: {
646: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
647: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
648: PetscErrorCode ierr;
649: PetscBool is_bfgs, is_dfp, is_symbrdn, is_symbadbrdn;
652: PetscObjectTypeCompare((PetscObject)B, MATLMVMBFGS, &is_bfgs);
653: PetscObjectTypeCompare((PetscObject)B, MATLMVMDFP, &is_dfp);
654: PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBROYDEN, &is_symbrdn);
655: PetscObjectTypeCompare((PetscObject)B, MATLMVMSYMBADBROYDEN, &is_symbadbrdn);
656: if (!is_bfgs && !is_dfp && !is_symbrdn && !is_symbadbrdn) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_INCOMP, "diagonal scaling is only available for DFP, BFGS and SymBrdn matrices");
657: lsb->delta = PetscAbsReal(PetscRealPart(delta));
658: lsb->delta = PetscMin(lsb->delta, lsb->delta_max);
659: lsb->delta = PetscMax(lsb->delta, lsb->delta_min);
660: return(0);
661: }
663: /*------------------------------------------------------------*/
665: /*@
666: MatLMVMSymBroydenSetScaleType - Sets the scale type for symmetric Broyden-type updates.
668: Input Parameters:
669: + snes - the iterative context
670: - rtype - restart type
672: Options Database:
673: . -mat_lmvm_scale_type <none,scalar,diagonal> - set the scaling type
675: Level: intermediate
677: MatLMVMSymBrdnScaleTypes:
678: + MAT_LMVM_SYMBROYDEN_SCALE_NONE - initial Hessian is the identity matrix
679: . MAT_LMVM_SYMBROYDEN_SCALE_SCALAR - use the Shanno scalar as the initial Hessian
680: - MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL - use a diagonalized BFGS update as the initial Hessian
682: .seealso: MATLMVMSYMBROYDEN, MatCreateLMVMSymBroyden()
683: @*/
684: PetscErrorCode MatLMVMSymBroydenSetScaleType(Mat B, MatLMVMSymBroydenScaleType stype)
685: {
686: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
687: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
691: lsb->scale_type = stype;
692: return(0);
693: }
695: /*------------------------------------------------------------*/
697: /*@
698: MatCreateLMVMSymBroyden - Creates a limited-memory Symmetric Broyden-type matrix used
699: for approximating Jacobians. L-SymBrdn is a convex combination of L-DFP and
700: L-BFGS such that SymBrdn = (1 - phi)*BFGS + phi*DFP. The combination factor
701: phi is restricted to the range [0, 1], where the L-SymBrdn matrix is guaranteed
702: to be symmetric positive-definite.
704: The provided local and global sizes must match the solution and function vectors
705: used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have
706: storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
707: parallel. To use the L-SymBrdn matrix with other vector types, the matrix must be
708: created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
709: This ensures that the internal storage and work vectors are duplicated from the
710: correct type of vector.
712: Collective
714: Input Parameters:
715: + comm - MPI communicator, set to PETSC_COMM_SELF
716: . n - number of local rows for storage vectors
717: - N - global size of the storage vectors
719: Output Parameter:
720: . B - the matrix
722: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
723: paradigm instead of this routine directly.
725: Options Database Keys:
726: + -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
727: . -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
728: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
729: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
730: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
731: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
732: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
733: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
735: Level: intermediate
737: .seealso: MatCreate(), MATLMVM, MATLMVMSYMBROYDEN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
738: MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn()
739: @*/
740: PetscErrorCode MatCreateLMVMSymBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
741: {
742: PetscErrorCode ierr;
745: MatCreate(comm, B);
746: MatSetSizes(*B, n, n, N, N);
747: MatSetType(*B, MATLMVMSYMBROYDEN);
748: MatSetUp(*B);
749: return(0);
750: }
752: /*------------------------------------------------------------*/
754: PetscErrorCode MatSymBrdnApplyJ0Fwd(Mat B, Vec X, Vec Z)
755: {
756: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
757: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
758: PetscErrorCode ierr;
761: if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
762: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
763: MatLMVMApplyJ0Fwd(B, X, Z);
764: } else {
765: switch (lsb->scale_type) {
766: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
767: VecCopy(X, Z);
768: VecScale(Z, 1.0/lsb->sigma);
769: break;
770: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
771: MatMult(lsb->D, X, Z);
772: break;
773: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
774: default:
775: VecCopy(X, Z);
776: break;
777: }
778: }
779: return(0);
780: }
782: /*------------------------------------------------------------*/
784: PetscErrorCode MatSymBrdnApplyJ0Inv(Mat B, Vec F, Vec dX)
785: {
786: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
787: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
788: PetscErrorCode ierr;
791: if (lmvm->J0 || lmvm->user_pc || lmvm->user_ksp || lmvm->user_scale) {
792: lsb->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_USER;
793: MatLMVMApplyJ0Inv(B, F, dX);
794: } else {
795: switch (lsb->scale_type) {
796: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
797: VecCopy(F, dX);
798: VecScale(dX, lsb->sigma);
799: break;
800: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
801: MatSolve(lsb->D, F, dX);
802: break;
803: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
804: default:
805: VecCopy(F, dX);
806: break;
807: }
808: }
809: return(0);
810: }
812: /*------------------------------------------------------------*/
814: PetscErrorCode MatSymBrdnComputeJ0Scalar(Mat B)
815: {
816: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
817: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
818: PetscInt i, start;
819: PetscReal a, b, c, sig1, sig2, signew;
822: if (lsb->sigma_hist == 0) {
823: signew = 1.0;
824: } else {
825: start = PetscMax(0, lmvm->k-lsb->sigma_hist+1);
826: signew = 0.0;
827: if (lsb->alpha == 1.0) {
828: for (i = start; i <= lmvm->k; ++i) {
829: signew += lsb->yts[i]/lsb->yty[i];
830: }
831: } else if (lsb->alpha == 0.5) {
832: for (i = start; i <= lmvm->k; ++i) {
833: signew += lsb->sts[i]/lsb->yty[i];
834: }
835: signew = PetscSqrtReal(signew);
836: } else if (lsb->alpha == 0.0) {
837: for (i = start; i <= lmvm->k; ++i) {
838: signew += lsb->sts[i]/lsb->yts[i];
839: }
840: } else {
841: /* compute coefficients of the quadratic */
842: a = b = c = 0.0;
843: for (i = start; i <= lmvm->k; ++i) {
844: a += lsb->yty[i];
845: b += lsb->yts[i];
846: c += lsb->sts[i];
847: }
848: a *= lsb->alpha;
849: b *= -(2.0*lsb->alpha - 1.0);
850: c *= lsb->alpha - 1.0;
851: /* use quadratic formula to find roots */
852: sig1 = (-b + PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
853: sig2 = (-b - PetscSqrtReal(b*b - 4.0*a*c))/(2.0*a);
854: /* accept the positive root as the scalar */
855: if (sig1 > 0.0) {
856: signew = sig1;
857: } else if (sig2 > 0.0) {
858: signew = sig2;
859: } else {
860: SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
861: }
862: }
863: }
864: lsb->sigma = lsb->rho*signew + (1.0 - lsb->rho)*lsb->sigma;
865: return(0);
866: }