Actual source code: diagbrdn.c
1: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
3: /*------------------------------------------------------------*/
5: static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
6: {
7: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
8: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
10: VecCheckSameSize(F, 2, dX, 3);
11: VecCheckMatCompatible(B, dX, 3, F, 2);
12: VecPointwiseMult(dX, ldb->invD, F);
13: return 0;
14: }
16: /*------------------------------------------------------------*/
18: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
19: {
20: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
21: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
23: VecCheckSameSize(X, 2, Z, 3);
24: VecCheckMatCompatible(B, X, 2, Z, 3);
25: VecPointwiseDivide(Z, X, ldb->invD);
26: return 0;
27: }
29: /*------------------------------------------------------------*/
31: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
32: {
33: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
34: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
35: PetscInt old_k, i, start;
36: PetscScalar yty, ststmp, curvature, ytDy, stDs, ytDs;
37: PetscReal curvtol, sigma, yy_sum, ss_sum, ys_sum, denom;
39: if (!lmvm->m) return 0;
40: if (lmvm->prev_set) {
41: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
42: VecAYPX(lmvm->Xprev, -1.0, X);
43: VecAYPX(lmvm->Fprev, -1.0, F);
44: /* Compute tolerance for accepting the update */
45: VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
46: VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
47: VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
48: VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
49: if (PetscRealPart(ststmp) < lmvm->eps) {
50: curvtol = 0.0;
51: } else {
52: curvtol = lmvm->eps * PetscRealPart(ststmp);
53: }
54: /* Test the curvature for the update */
55: if (PetscRealPart(curvature) > curvtol) {
56: /* Update is good so we accept it */
57: old_k = lmvm->k;
58: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
59: /* If we hit the memory limit, shift the yty and yts arrays */
60: if (old_k == lmvm->k) {
61: for (i = 0; i <= lmvm->k-1; ++i) {
62: ldb->yty[i] = ldb->yty[i+1];
63: ldb->yts[i] = ldb->yts[i+1];
64: ldb->sts[i] = ldb->sts[i+1];
65: }
66: }
67: /* Accept dot products into the history */
68: VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
69: ldb->yty[lmvm->k] = PetscRealPart(yty);
70: ldb->yts[lmvm->k] = PetscRealPart(curvature);
71: ldb->sts[lmvm->k] = PetscRealPart(ststmp);
72: if (ldb->forward) {
73: /* We are doing diagonal scaling of the forward Hessian B */
74: /* BFGS = DFP = inv(D); */
75: VecCopy(ldb->invD, ldb->invDnew);
76: VecReciprocal(ldb->invDnew);
78: /* V = y*y */
79: VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]);
81: /* W = inv(D)*s */
82: VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]);
83: VecDot(ldb->W, lmvm->S[lmvm->k], &stDs);
85: /* Safeguard stDs */
86: stDs = PetscMax(PetscRealPart(stDs), ldb->tol);
88: if (1.0 != ldb->theta) {
89: /* BFGS portion of the update */
90: /* U = (inv(D)*s)*(inv(D)*s) */
91: VecPointwiseMult(ldb->U, ldb->W, ldb->W);
93: /* Assemble */
94: VecAXPBY(ldb->BFGS, -1.0/stDs, 0.0, ldb->U);
95: }
96: if (0.0 != ldb->theta) {
97: /* DFP portion of the update */
98: /* U = inv(D)*s*y */
99: VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]);
101: /* Assemble */
102: VecAXPBY(ldb->DFP, stDs/ldb->yts[lmvm->k], 0.0, ldb->V);
103: VecAXPY(ldb->DFP, -2.0, ldb->U);
104: }
106: if (0.0 == ldb->theta) {
107: VecAXPY(ldb->invDnew, 1.0, ldb->BFGS);
108: } else if (1.0 == ldb->theta) {
109: VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->DFP);
110: } else {
111: /* Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
112: VecAXPBYPCZ(ldb->invDnew, 1.0-ldb->theta, (ldb->theta)/ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP);
113: }
115: VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
116: /* Obtain inverse and ensure positive definite */
117: VecReciprocal(ldb->invDnew);
118: VecAbs(ldb->invDnew);
120: } else {
121: /* Inverse Hessian update instead. */
122: VecCopy(ldb->invD, ldb->invDnew);
124: /* V = s*s */
125: VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]);
127: /* W = D*y */
128: VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]);
129: VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy);
131: /* Safeguard ytDy */
132: ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);
134: if (1.0 != ldb->theta) {
135: /* BFGS portion of the update */
136: /* U = s*Dy */
137: VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]);
139: /* Assemble */
140: VecAXPBY(ldb->BFGS, ytDy/ldb->yts[lmvm->k], 0.0, ldb->V);
141: VecAXPY(ldb->BFGS, -2.0, ldb->U);
142: }
143: if (0.0 != ldb->theta) {
144: /* DFP portion of the update */
146: /* U = (inv(D)*y)*(inv(D)*y) */
147: VecPointwiseMult(ldb->U, ldb->W, ldb->W);
149: /* Assemble */
150: VecAXPBY(ldb->DFP, -1.0/ytDy, 0.0, ldb->U);
151: }
153: if (0.0 == ldb->theta) {
154: VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->BFGS);
155: } else if (1.0 == ldb->theta) {
156: VecAXPY(ldb->invDnew, 1.0, ldb->DFP);
157: } else {
158: /* Broyden update U=(1-theta)*P + theta*Q */
159: VecAXPBYPCZ(ldb->invDnew, (1.0-ldb->theta)/ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP);
160: }
161: VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
162: /* Ensure positive definite */
163: VecAbs(ldb->invDnew);
164: }
165: if (ldb->sigma_hist > 0) {
166: /* Start with re-scaling on the newly computed diagonal */
167: if (0.5 == ldb->beta) {
168: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
169: VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew);
170: VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);
172: VecDotBegin(ldb->V, lmvm->Y[0], &ytDy);
173: VecDotBegin(ldb->W, lmvm->S[0], &stDs);
174: VecDotEnd(ldb->V, lmvm->Y[0], &ytDy);
175: VecDotEnd(ldb->W, lmvm->S[0], &stDs);
177: ss_sum = PetscRealPart(stDs);
178: yy_sum = PetscRealPart(ytDy);
179: ys_sum = ldb->yts[0];
180: } else {
181: VecCopy(ldb->invDnew, ldb->U);
182: VecReciprocal(ldb->U);
184: /* Compute summations for scalar scaling */
185: yy_sum = 0; /* No safeguard required */
186: ys_sum = 0; /* No safeguard required */
187: ss_sum = 0; /* No safeguard required */
188: start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
189: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
190: VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U);
191: VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);
193: VecDotBegin(ldb->W, lmvm->S[i], &stDs);
194: VecDotBegin(ldb->V, lmvm->Y[i], &ytDy);
195: VecDotEnd(ldb->W, lmvm->S[i], &stDs);
196: VecDotEnd(ldb->V, lmvm->Y[i], &ytDy);
198: ss_sum += PetscRealPart(stDs);
199: ys_sum += ldb->yts[i];
200: yy_sum += PetscRealPart(ytDy);
201: }
202: }
203: } else if (0.0 == ldb->beta) {
204: if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
205: /* Compute summations for scalar scaling */
206: VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);
208: VecDotBegin(ldb->W, lmvm->Y[0], &ytDs);
209: VecDotBegin(ldb->W, ldb->W, &stDs);
210: VecDotEnd(ldb->W, lmvm->Y[0], &ytDs);
211: VecDotEnd(ldb->W, ldb->W, &stDs);
213: ys_sum = PetscRealPart(ytDs);
214: ss_sum = PetscRealPart(stDs);
215: yy_sum = ldb->yty[0];
216: } else {
217: VecCopy(ldb->invDnew, ldb->U);
218: VecReciprocal(ldb->U);
220: /* Compute summations for scalar scaling */
221: yy_sum = 0; /* No safeguard required */
222: ys_sum = 0; /* No safeguard required */
223: ss_sum = 0; /* No safeguard required */
224: start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
225: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
226: VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);
228: VecDotBegin(ldb->W, lmvm->Y[i], &ytDs);
229: VecDotBegin(ldb->W, ldb->W, &stDs);
230: VecDotEnd(ldb->W, lmvm->Y[i], &ytDs);
231: VecDotEnd(ldb->W, ldb->W, &stDs);
233: ss_sum += PetscRealPart(stDs);
234: ys_sum += PetscRealPart(ytDs);
235: yy_sum += ldb->yty[i];
236: }
237: }
238: } else if (1.0 == ldb->beta) {
239: /* Compute summations for scalar scaling */
240: yy_sum = 0; /* No safeguard required */
241: ys_sum = 0; /* No safeguard required */
242: ss_sum = 0; /* No safeguard required */
243: start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
244: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
245: VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew);
247: VecDotBegin(ldb->V, lmvm->S[i], &ytDs);
248: VecDotBegin(ldb->V, ldb->V, &ytDy);
249: VecDotEnd(ldb->V, lmvm->S[i], &ytDs);
250: VecDotEnd(ldb->V, ldb->V, &ytDy);
252: yy_sum += PetscRealPart(ytDy);
253: ys_sum += PetscRealPart(ytDs);
254: ss_sum += ldb->sts[i];
255: }
256: } else {
257: VecCopy(ldb->invDnew, ldb->U);
258: VecPow(ldb->U, ldb->beta-1);
260: /* Compute summations for scalar scaling */
261: yy_sum = 0; /* No safeguard required */
262: ys_sum = 0; /* No safeguard required */
263: ss_sum = 0; /* No safeguard required */
264: start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
265: for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
266: VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]);
267: VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]);
269: VecDotBegin(ldb->V, ldb->W, &ytDs);
270: VecDotBegin(ldb->V, ldb->V, &ytDy);
271: VecDotBegin(ldb->W, ldb->W, &stDs);
272: VecDotEnd(ldb->V, ldb->W, &ytDs);
273: VecDotEnd(ldb->V, ldb->V, &ytDy);
274: VecDotEnd(ldb->W, ldb->W, &stDs);
276: yy_sum += PetscRealPart(ytDy);
277: ys_sum += PetscRealPart(ytDs);
278: ss_sum += PetscRealPart(stDs);
279: }
280: }
282: if (0.0 == ldb->alpha) {
283: /* Safeguard ys_sum */
284: ys_sum = PetscMax(ldb->tol, ys_sum);
286: sigma = ss_sum / ys_sum;
287: } else if (1.0 == ldb->alpha) {
288: /* yy_sum is never 0; if it were, we'd be at the minimum */
289: sigma = ys_sum / yy_sum;
290: } else {
291: denom = 2.0*ldb->alpha*yy_sum;
293: /* Safeguard denom */
294: denom = PetscMax(ldb->tol, denom);
296: sigma = ((2.0*ldb->alpha-1)*ys_sum + PetscSqrtReal((2.0*ldb->alpha-1)*(2.0*ldb->alpha-1)*ys_sum*ys_sum - 4.0*ldb->alpha*(ldb->alpha-1)*yy_sum*ss_sum)) / denom;
297: }
298: } else {
299: sigma = 1.0;
300: }
301: /* If Q has small values, then Q^(r_beta - 1)
302: can have very large values. Hence, ys_sum
303: and ss_sum can be infinity. In this case,
304: sigma can either be not-a-number or infinity. */
306: if (PetscIsInfOrNanScalar(sigma)) {
307: /* sigma is not-a-number; skip rescaling */
308: } else if (0.0 == sigma) {
309: /* sigma is zero; this is a bad case; skip rescaling */
310: } else {
311: /* sigma is positive */
312: VecScale(ldb->invDnew, sigma);
313: }
315: /* Combine the old diagonal and the new diagonal using a convex limiter */
316: if (1.0 == ldb->rho) {
317: VecCopy(ldb->invDnew, ldb->invD);
318: } else if (ldb->rho) {
319: VecAXPBY(ldb->invD, 1.0-ldb->rho, ldb->rho, ldb->invDnew);
320: }
321: } else {
322: MatLMVMReset(B, PETSC_FALSE);
323: }
324: /* End DiagBrdn update */
326: }
327: /* Save the solution and function to be used in the next update */
328: VecCopy(X, lmvm->Xprev);
329: VecCopy(F, lmvm->Fprev);
330: lmvm->prev_set = PETSC_TRUE;
331: return 0;
332: }
334: /*------------------------------------------------------------*/
336: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
337: {
338: Mat_LMVM *bdata = (Mat_LMVM*)B->data;
339: Mat_DiagBrdn *bctx = (Mat_DiagBrdn*)bdata->ctx;
340: Mat_LMVM *mdata = (Mat_LMVM*)M->data;
341: Mat_DiagBrdn *mctx = (Mat_DiagBrdn*)mdata->ctx;
342: PetscInt i;
344: mctx->theta = bctx->theta;
345: mctx->alpha = bctx->alpha;
346: mctx->beta = bctx->beta;
347: mctx->rho = bctx->rho;
348: mctx->delta = bctx->delta;
349: mctx->delta_min = bctx->delta_min;
350: mctx->delta_max = bctx->delta_max;
351: mctx->tol = bctx->tol;
352: mctx->sigma = bctx->sigma;
353: mctx->sigma_hist = bctx->sigma_hist;
354: mctx->forward = bctx->forward;
355: VecCopy(bctx->invD, mctx->invD);
356: for (i=0; i<=bdata->k; ++i) {
357: mctx->yty[i] = bctx->yty[i];
358: mctx->yts[i] = bctx->yts[i];
359: mctx->sts[i] = bctx->sts[i];
360: }
361: return 0;
362: }
364: /*------------------------------------------------------------*/
366: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
367: {
368: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
369: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
370: PetscBool isascii;
372: PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
373: if (isascii) {
374: PetscViewerASCIIPrintf(pv,"Scale history: %d\n",ldb->sigma_hist);
375: PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)ldb->alpha, (double)ldb->beta, (double)ldb->rho);
376: PetscViewerASCIIPrintf(pv,"Convex factor: theta=%g\n", (double)ldb->theta);
377: }
378: MatView_LMVM(B, pv);
379: return 0;
380: }
382: /*------------------------------------------------------------*/
384: static PetscErrorCode MatSetFromOptions_DiagBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
385: {
386: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
387: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
389: MatSetFromOptions_LMVM(PetscOptionsObject, B);
390: PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
391: PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",ldb->theta,&ldb->theta,NULL);
392: PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",ldb->rho,&ldb->rho,NULL);
393: PetscOptionsReal("-mat_lmvm_tol","(developer) tolerance for bounding rescaling denominator","",ldb->tol,&ldb->tol,NULL);
394: PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",ldb->alpha,&ldb->alpha,NULL);
395: PetscOptionsBool("-mat_lmvm_forward","Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.","",ldb->forward,&ldb->forward,NULL);
396: PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",ldb->beta,&ldb->beta,NULL);
397: PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",ldb->sigma_hist,&ldb->sigma_hist,NULL);
398: PetscOptionsTail();
403: return 0;
404: }
406: /*------------------------------------------------------------*/
408: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
409: {
410: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
411: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
413: VecSet(ldb->invD, ldb->delta);
414: if (destructive && ldb->allocated) {
415: PetscFree3(ldb->yty, ldb->yts, ldb->sts);
416: VecDestroy(&ldb->invDnew);
417: VecDestroy(&ldb->invD);
418: VecDestroy(&ldb->BFGS);
419: VecDestroy(&ldb->DFP);
420: VecDestroy(&ldb->U);
421: VecDestroy(&ldb->V);
422: VecDestroy(&ldb->W);
423: ldb->allocated = PETSC_FALSE;
424: }
425: MatReset_LMVM(B, destructive);
426: return 0;
427: }
429: /*------------------------------------------------------------*/
431: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
432: {
433: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
434: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
436: MatAllocate_LMVM(B, X, F);
437: if (!ldb->allocated) {
438: PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
439: VecDuplicate(lmvm->Xprev, &ldb->invDnew);
440: VecDuplicate(lmvm->Xprev, &ldb->invD);
441: VecDuplicate(lmvm->Xprev, &ldb->BFGS);
442: VecDuplicate(lmvm->Xprev, &ldb->DFP);
443: VecDuplicate(lmvm->Xprev, &ldb->U);
444: VecDuplicate(lmvm->Xprev, &ldb->V);
445: VecDuplicate(lmvm->Xprev, &ldb->W);
446: ldb->allocated = PETSC_TRUE;
447: }
448: return 0;
449: }
451: /*------------------------------------------------------------*/
453: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
454: {
455: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
456: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
458: if (ldb->allocated) {
459: PetscFree3(ldb->yty, ldb->yts, ldb->sts);
460: VecDestroy(&ldb->invDnew);
461: VecDestroy(&ldb->invD);
462: VecDestroy(&ldb->BFGS);
463: VecDestroy(&ldb->DFP);
464: VecDestroy(&ldb->U);
465: VecDestroy(&ldb->V);
466: VecDestroy(&ldb->W);
467: ldb->allocated = PETSC_FALSE;
468: }
469: PetscFree(lmvm->ctx);
470: MatDestroy_LMVM(B);
471: return 0;
472: }
474: /*------------------------------------------------------------*/
476: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
477: {
478: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
479: Mat_DiagBrdn *ldb = (Mat_DiagBrdn*)lmvm->ctx;
481: MatSetUp_LMVM(B);
482: if (!ldb->allocated) {
483: PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
484: VecDuplicate(lmvm->Xprev, &ldb->invDnew);
485: VecDuplicate(lmvm->Xprev, &ldb->invD);
486: VecDuplicate(lmvm->Xprev, &ldb->BFGS);
487: VecDuplicate(lmvm->Xprev, &ldb->DFP);
488: VecDuplicate(lmvm->Xprev, &ldb->U);
489: VecDuplicate(lmvm->Xprev, &ldb->V);
490: VecDuplicate(lmvm->Xprev, &ldb->W);
491: ldb->allocated = PETSC_TRUE;
492: }
493: return 0;
494: }
496: /*------------------------------------------------------------*/
498: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
499: {
500: Mat_LMVM *lmvm;
501: Mat_DiagBrdn *ldb;
503: MatCreate_LMVM(B);
504: PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBROYDEN);
505: B->ops->setup = MatSetUp_DiagBrdn;
506: B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
507: B->ops->destroy = MatDestroy_DiagBrdn;
508: B->ops->solve = MatSolve_DiagBrdn;
509: B->ops->view = MatView_DiagBrdn;
511: lmvm = (Mat_LMVM*)B->data;
512: lmvm->square = PETSC_TRUE;
513: lmvm->m = 1;
514: lmvm->ops->allocate = MatAllocate_DiagBrdn;
515: lmvm->ops->reset = MatReset_DiagBrdn;
516: lmvm->ops->mult = MatMult_DiagBrdn;
517: lmvm->ops->update = MatUpdate_DiagBrdn;
518: lmvm->ops->copy = MatCopy_DiagBrdn;
520: PetscNewLog(B, &ldb);
521: lmvm->ctx = (void*)ldb;
522: ldb->theta = 0.0;
523: ldb->alpha = 1.0;
524: ldb->rho = 1.0;
525: ldb->forward = PETSC_TRUE;
526: ldb->beta = 0.5;
527: ldb->sigma = 1.0;
528: ldb->delta = 1.0;
529: ldb->delta_min = 1e-7;
530: ldb->delta_max = 100.0;
531: ldb->tol = 1e-8;
532: ldb->sigma_hist = 1;
533: ldb->allocated = PETSC_FALSE;
534: return 0;
535: }
537: /*------------------------------------------------------------*/
539: /*@
540: MatCreateLMVMDiagBroyden - DiagBrdn creates a symmetric Broyden-type diagonal matrix used
541: for approximating Hessians. It consists of a convex combination of DFP and BFGS
542: diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP.
543: To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1].
544: We also ensure positive definiteness by taking the VecAbs() of the final vector.
546: There are two ways of approximating the diagonal: using the forward (B) update
547: schemes for BFGS and DFP and then taking the inverse, or directly working with
548: the inverse (H) update schemes for the BFGS and DFP updates, derived using the
549: Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a
550: parameter below.
552: In order to use the DiagBrdn matrix with other vector types, i.e. doing MatMults
553: and MatSolves, the matrix must first be created using MatCreate() and MatSetType(),
554: followed by MatLMVMAllocate(). Then it will be available for updating
555: (via MatLMVMUpdate) in one's favored solver implementation.
556: This allows for MPI compatibility.
558: Collective
560: Input Parameters:
561: + comm - MPI communicator, set to PETSC_COMM_SELF
562: . n - number of local rows for storage vectors
563: - N - global size of the storage vectors
565: Output Parameter:
566: . B - the matrix
568: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
569: paradigm instead of this routine directly.
571: Options Database Keys:
572: + -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
573: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
574: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
575: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
576: . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
577: . -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
578: - -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal
580: Level: intermediate
582: .seealso: MatCreate(), MATLMVM, MATLMVMDIAGBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
583: MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMSymBrdn()
584: @*/
585: PetscErrorCode MatCreateLMVMDiagBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
586: {
587: MatCreate(comm, B);
588: MatSetSizes(*B, n, n, N, N);
589: MatSetType(*B, MATLMVMDIAGBROYDEN);
590: MatSetUp(*B);
591: return 0;
592: }