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