Actual source code: symbrdnrescale.c
1: #include <petscdevice.h>
2: #include "symbrdnrescale.h"
4: PetscLogEvent SBRDN_Rescale;
6: const char *const MatLMVMSymBroydenScaleTypes[] = {"none", "scalar", "diagonal", "user", "decide", "MatLMVMSymBroydenScaleType", "MAT_LMVM_SYMBROYDEN_SCALING_", NULL};
8: static PetscErrorCode SymBroydenRescaleUpdateScalar(Mat B, SymBroydenRescale ldb)
9: {
10: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
11: PetscReal a, b, c, signew;
12: PetscReal sigma_inv, sigma;
13: PetscInt oldest, next;
15: PetscFunctionBegin;
16: next = ldb->k;
17: oldest = PetscMax(0, ldb->k - ldb->sigma_hist);
18: PetscCall(MatNorm(lmvm->J0, NORM_INFINITY, &sigma_inv));
19: sigma = 1.0 / sigma_inv;
20: if (ldb->sigma_hist == 0) {
21: signew = 1.0;
22: } else {
23: signew = 0.0;
24: if (ldb->alpha == 1.0) {
25: for (PetscInt i = 0; i < next - oldest; ++i) signew += ldb->yts[i] / ldb->yty[i];
26: } else if (ldb->alpha == 0.5) {
27: for (PetscInt i = 0; i < next - oldest; ++i) signew += ldb->sts[i] / ldb->yty[i];
28: signew = PetscSqrtReal(signew);
29: } else if (ldb->alpha == 0.0) {
30: for (PetscInt i = 0; i < next - oldest; ++i) signew += ldb->sts[i] / ldb->yts[i];
31: } else {
32: /* compute coefficients of the quadratic */
33: a = b = c = 0.0;
34: for (PetscInt i = 0; i < next - oldest; ++i) {
35: a += ldb->yty[i];
36: b += ldb->yts[i];
37: c += ldb->sts[i];
38: }
39: a *= ldb->alpha;
40: b *= -(2.0 * ldb->alpha - 1.0);
41: c *= ldb->alpha - 1.0;
42: /* use quadratic formula to find roots */
43: PetscReal sqrtdisc = PetscSqrtReal(b * b - 4 * a * c);
44: if (b >= 0.0) {
45: if (a >= 0.0) {
46: signew = (2 * c) / (-b - sqrtdisc);
47: } else {
48: signew = (-b - sqrtdisc) / (2 * a);
49: }
50: } else {
51: if (a >= 0.0) {
52: signew = (-b + sqrtdisc) / (2 * a);
53: } else {
54: signew = (2 * c) / (-b + sqrtdisc);
55: }
56: }
57: PetscCheck(signew > 0.0, PetscObjectComm((PetscObject)B), PETSC_ERR_CONV_FAILED, "Cannot find positive scalar");
58: }
59: }
60: sigma = ldb->rho * signew + (1.0 - ldb->rho) * sigma;
61: PetscCall(MatLMVMSetJ0Scale(B, 1.0 / sigma));
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }
65: static PetscErrorCode DiagonalUpdate(SymBroydenRescale ldb, Vec D, Vec s, Vec y, Vec V, Vec W, Vec BFGS, Vec DFP, PetscReal theta, PetscReal yts)
66: {
67: PetscFunctionBegin;
68: /* V = |y o y| */
69: PetscCall(VecPointwiseMult(V, y, y));
70: if (PetscDefined(USE_COMPLEX)) PetscCall(VecAbs(V));
72: /* W = D o s */
73: PetscReal stDs;
74: PetscCall(VecPointwiseMult(W, D, s));
75: PetscCall(VecDotRealPart(W, s, &stDs));
77: PetscCall(VecAXPY(D, 1.0 / yts, ldb->V));
79: /* Safeguard stDs */
80: stDs = PetscMax(stDs, ldb->tol);
82: if (theta != 1.0) {
83: /* BFGS portion of the update */
85: /* U = |(D o s) o (D o s)| */
86: PetscCall(VecPointwiseMult(BFGS, W, W));
87: if (PetscDefined(USE_COMPLEX)) PetscCall(VecAbs(BFGS));
89: /* Assemble */
90: PetscCall(VecScale(BFGS, -1.0 / stDs));
91: }
93: if (theta != 0.0) {
94: /* DFP portion of the update */
95: /* U = Real(conj(y) o D o s) */
96: PetscCall(VecCopy(y, DFP));
97: PetscCall(VecConjugate(DFP));
98: PetscCall(VecPointwiseMult(DFP, DFP, W));
99: if (PetscDefined(USE_COMPLEX)) {
100: PetscCall(VecCopy(DFP, W));
101: PetscCall(VecConjugate(W));
102: PetscCall(VecAXPY(DFP, 1.0, W));
103: } else {
104: PetscCall(VecScale(DFP, 2.0));
105: }
107: /* Assemble */
108: PetscCall(VecAXPBY(DFP, stDs / yts, -1.0, V));
109: }
111: if (theta == 0.0) PetscCall(VecAXPY(D, 1.0, BFGS));
112: else if (theta == 1.0) PetscCall(VecAXPY(D, 1.0 / yts, DFP));
113: else {
114: /* Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
115: PetscCall(VecAXPBYPCZ(D, 1.0 - theta, theta / yts, 1.0, BFGS, DFP));
116: }
117: PetscFunctionReturn(PETSC_SUCCESS);
118: }
120: static PetscErrorCode SymBroydenRescaleUpdateDiagonal(Mat B, SymBroydenRescale ldb)
121: {
122: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
123: PetscInt oldest, next;
124: Vec invD, s_last, y_last;
125: LMBasis S, Y;
126: PetscScalar yts;
127: PetscReal sigma;
129: PetscFunctionBegin;
130: next = ldb->k;
131: oldest = PetscMax(0, ldb->k - ldb->sigma_hist);
132: PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_S, next - 1, &yts));
133: PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_S, &S, NULL, NULL));
134: PetscCall(MatLMVMGetUpdatedBasis(B, LMBASIS_Y, &Y, NULL, NULL));
135: PetscCall(LMBasisGetVecRead(S, next - 1, &s_last));
136: PetscCall(LMBasisGetVecRead(Y, next - 1, &y_last));
137: PetscCall(MatLMVMGetJ0InvDiag(B, &invD));
138: if (ldb->forward) {
139: /* We are doing diagonal scaling of the forward Hessian B */
140: /* BFGS = DFP = inv(D); */
141: PetscCall(VecCopy(invD, ldb->invDnew));
142: PetscCall(VecReciprocal(ldb->invDnew));
143: PetscCall(DiagonalUpdate(ldb, ldb->invDnew, s_last, y_last, ldb->V, ldb->W, ldb->BFGS, ldb->DFP, ldb->theta, PetscRealPart(yts)));
144: /* Obtain inverse and ensure positive definite */
145: PetscCall(VecReciprocal(ldb->invDnew));
146: } else {
147: /* Inverse Hessian update instead. */
148: PetscCall(VecCopy(invD, ldb->invDnew));
149: PetscCall(DiagonalUpdate(ldb, ldb->invDnew, y_last, s_last, ldb->V, ldb->W, ldb->DFP, ldb->BFGS, 1.0 - ldb->theta, PetscRealPart(yts)));
150: }
151: PetscCall(VecAbs(ldb->invDnew));
152: PetscCall(LMBasisRestoreVecRead(Y, next - 1, &y_last));
153: PetscCall(LMBasisRestoreVecRead(S, next - 1, &s_last));
155: if (ldb->sigma_hist > 0) {
156: // We are computing the scaling factor sigma that minimizes
157: //
158: // Sum_i || sigma^(alpha) (D^(-beta) o y_i) - sigma^(alpha-1) (D^(1-beta) o s_i) ||_2^2
159: // `-------.-------' `--------.-------'
160: // v_i w_i
161: //
162: // To do this we first have to compute the sums of the dot product terms
163: //
164: // yy_sum = Sum_i v_i^T v_i,
165: // ys_sum = Sum_i v_i^T w_i, and
166: // ss_sum = Sum_i w_i^T w_i.
167: //
168: // These appear in the quadratic equation for the optimality condition for sigma,
169: //
170: // [alpha yy_sum] sigma^2 - [(2 alpha - 1) ys_sum] * sigma + [(alpha - 1) * ss_sum] = 0
171: //
172: // which we solve for sigma.
174: PetscReal yy_sum = 0; /* No safeguard required */
175: PetscReal ys_sum = 0; /* No safeguard required */
176: PetscReal ss_sum = 0; /* No safeguard required */
177: PetscInt start = PetscMax(oldest, lmvm->k - ldb->sigma_hist);
179: Vec D_minus_beta = NULL;
180: Vec D_minus_beta_squared = NULL;
181: Vec D_one_minus_beta = NULL;
182: Vec D_one_minus_beta_squared = NULL;
183: if (ldb->beta == 0.5) {
184: D_minus_beta_squared = ldb->invDnew; // (D^(-0.5))^2 = D^-1
186: PetscCall(VecCopy(ldb->invDnew, ldb->U));
187: PetscCall(VecReciprocal(ldb->U));
188: D_one_minus_beta_squared = ldb->U; // (D^(1-0.5))^2 = D
189: } else if (ldb->beta == 0.0) {
190: PetscCall(VecCopy(ldb->invDnew, ldb->U));
191: PetscCall(VecReciprocal(ldb->U));
192: D_one_minus_beta = ldb->U; // D^1
193: } else if (ldb->beta == 1.0) {
194: D_minus_beta = ldb->invDnew; // D^-1
195: } else {
196: PetscCall(VecCopy(ldb->invDnew, ldb->DFP));
197: PetscCall(VecPow(ldb->DFP, ldb->beta));
198: D_minus_beta = ldb->DFP;
200: PetscCall(VecCopy(ldb->invDnew, ldb->BFGS));
201: PetscCall(VecPow(ldb->BFGS, ldb->beta - 1));
202: D_one_minus_beta = ldb->BFGS;
203: }
204: for (PetscInt i = start - oldest; i < next - oldest; ++i) {
205: Vec s_i, y_i;
207: PetscCall(LMBasisGetVecRead(S, oldest + i, &s_i));
208: PetscCall(LMBasisGetVecRead(Y, oldest + i, &y_i));
209: if (ldb->beta == 0.5) {
210: PetscReal ytDinvy, stDs;
212: PetscCall(VecPointwiseMult(ldb->V, y_i, D_minus_beta_squared));
213: PetscCall(VecPointwiseMult(ldb->W, s_i, D_one_minus_beta_squared));
214: PetscCall(VecDotRealPart(ldb->W, s_i, &stDs));
215: PetscCall(VecDotRealPart(ldb->V, y_i, &ytDinvy));
216: ss_sum += stDs; // ||s||_{D^(2*(1-beta))}^2
217: ys_sum += ldb->yts[i]; // s^T D^(1 - 2*beta) y
218: yy_sum += ytDinvy; // ||y||_{D^(-2*beta)}^2
219: } else if (ldb->beta == 0.0) {
220: PetscScalar ytDs_scalar;
221: PetscReal stDsr;
223: PetscCall(VecPointwiseMult(ldb->W, s_i, D_one_minus_beta));
224: PetscCall(VecDotNorm2(y_i, ldb->W, &ytDs_scalar, &stDsr));
225: ss_sum += stDsr; // ||s||_{D^(2*(1-beta))}^2
226: ys_sum += PetscRealPart(ytDs_scalar); // s^T D^(1 - 2*beta) y
227: yy_sum += ldb->yty[i]; // ||y||_{D^(-2*beta)}^2
228: } else if (ldb->beta == 1.0) {
229: PetscScalar ytDs_scalar;
230: PetscReal ytDyr;
232: PetscCall(VecPointwiseMult(ldb->V, y_i, D_minus_beta));
233: PetscCall(VecDotNorm2(s_i, ldb->V, &ytDs_scalar, &ytDyr));
234: ss_sum += ldb->sts[i]; // ||s||_{D^(2*(1-beta))}^2
235: ys_sum += PetscRealPart(ytDs_scalar); // s^T D^(1 - 2*beta) y
236: yy_sum += ytDyr; // ||y||_{D^(-2*beta)}^2
237: } else {
238: PetscScalar ytDs_scalar;
239: PetscReal ytDyr, stDs;
241: PetscCall(VecPointwiseMult(ldb->V, y_i, D_minus_beta));
242: PetscCall(VecPointwiseMult(ldb->W, s_i, D_one_minus_beta));
243: PetscCall(VecDotNorm2(ldb->W, ldb->V, &ytDs_scalar, &ytDyr));
244: PetscCall(VecDotRealPart(ldb->W, ldb->W, &stDs));
245: ss_sum += stDs; // ||s||_{D^(2*(1-beta))}^2
246: ys_sum += PetscRealPart(ytDs_scalar); // s^T D^(1 - 2*beta) y
247: yy_sum += ytDyr; // ||y||_{D^(-2*beta)}^2
248: }
249: PetscCall(LMBasisRestoreVecRead(Y, oldest + i, &y_i));
250: PetscCall(LMBasisRestoreVecRead(S, oldest + i, &s_i));
251: }
253: if (ldb->alpha == 0.0) {
254: /* Safeguard ys_sum */
255: ys_sum = PetscMax(ldb->tol, ys_sum);
257: sigma = ss_sum / ys_sum;
258: } else if (1.0 == ldb->alpha) {
259: /* yy_sum is never 0; if it were, we'd be at the minimum */
260: sigma = ys_sum / yy_sum;
261: } else {
262: PetscReal a = ldb->alpha * yy_sum;
263: PetscReal b = -(2.0 * ldb->alpha - 1.0) * ys_sum;
264: PetscReal c = (ldb->alpha - 1.0) * ss_sum;
265: PetscReal sqrt_disc = PetscSqrtReal(b * b - 4 * a * c);
267: // numerically stable computation of positive root
268: if (b >= 0.0) {
269: if (a >= 0) {
270: PetscReal denom = PetscMax(-b - sqrt_disc, ldb->tol);
272: sigma = (2 * c) / denom;
273: } else {
274: PetscReal denom = PetscMax(2 * a, ldb->tol);
276: sigma = (-b - sqrt_disc) / denom;
277: }
278: } else {
279: if (a >= 0) {
280: PetscReal denom = PetscMax(2 * a, ldb->tol);
282: sigma = (-b + sqrt_disc) / denom;
283: } else {
284: PetscReal denom = PetscMax(-b + sqrt_disc, ldb->tol);
286: sigma = (2 * c) / denom;
287: }
288: }
289: }
290: } else {
291: sigma = 1.0;
292: }
293: /* If Q has small values, then Q^(r_beta - 1)
294: can have very large values. Hence, ys_sum
295: and ss_sum can be infinity. In this case,
296: sigma can either be not-a-number or infinity. */
297: if (PetscIsNormalReal(sigma)) PetscCall(VecScale(ldb->invDnew, sigma));
298: /* Combine the old diagonal and the new diagonal using a convex limiter */
299: if (ldb->rho == 1.0) PetscCall(VecCopy(ldb->invDnew, invD));
300: else if (ldb->rho) PetscCall(VecAXPBY(invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew));
301: PetscCall(MatLMVMRestoreJ0InvDiag(B, &invD));
302: PetscFunctionReturn(PETSC_SUCCESS);
303: }
305: static PetscErrorCode SymBroydenRescaleUpdateJ0(Mat B, SymBroydenRescale ldb)
306: {
307: PetscFunctionBegin;
308: if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(SymBroydenRescaleUpdateScalar(B, ldb));
309: else if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(SymBroydenRescaleUpdateDiagonal(B, ldb));
310: PetscFunctionReturn(PETSC_SUCCESS);
311: }
313: PETSC_INTERN PetscErrorCode SymBroydenRescaleUpdate(Mat B, SymBroydenRescale ldb)
314: {
315: PetscInt oldest, next;
317: PetscFunctionBegin;
318: PetscCall(SymBroydenRescaleSetUp(B, ldb));
319: if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_NONE || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_USER) PetscFunctionReturn(PETSC_SUCCESS);
320: PetscCall(PetscLogEventBegin(SBRDN_Rescale, NULL, NULL, NULL, NULL));
321: PetscCall(MatLMVMGetRange(B, &oldest, &next));
322: if (next > ldb->k) {
323: PetscInt new_oldest = PetscMax(0, next - ldb->sigma_hist);
324: PetscInt ldb_oldest = PetscMax(0, ldb->k - ldb->sigma_hist);
326: if (new_oldest > ldb_oldest) {
327: for (PetscInt i = new_oldest; i < ldb->k; i++) {
328: ldb->yty[i - new_oldest] = ldb->yty[i - ldb_oldest];
329: ldb->yts[i - new_oldest] = ldb->yts[i - ldb_oldest];
330: ldb->sts[i - new_oldest] = ldb->sts[i - ldb_oldest];
331: }
332: }
333: for (PetscInt i = PetscMax(new_oldest, ldb->k); i < next; i++) {
334: PetscScalar yty, sts, yts;
336: PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_Y, i, &yty));
337: PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_S, i, &yts));
338: PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_S, LMBASIS_S, i, &sts));
339: ldb->yty[i - new_oldest] = PetscRealPart(yty);
340: ldb->yts[i - new_oldest] = PetscRealPart(yts);
341: ldb->sts[i - new_oldest] = PetscRealPart(sts);
342: }
343: ldb->k = next;
344: }
345: PetscCall(SymBroydenRescaleUpdateJ0(B, ldb));
346: PetscCall(PetscLogEventEnd(SBRDN_Rescale, NULL, NULL, NULL, NULL));
347: PetscFunctionReturn(PETSC_SUCCESS);
348: }
350: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetDelta(Mat B, SymBroydenRescale ldb, PetscReal delta)
351: {
352: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
353: PetscBool same;
355: PetscFunctionBegin;
356: same = (delta == ldb->delta) ? PETSC_TRUE : PETSC_FALSE;
357: ldb->delta = delta;
358: ldb->delta = PetscMin(ldb->delta, ldb->delta_max);
359: ldb->delta = PetscMax(ldb->delta, ldb->delta_min);
360: // if there have been no updates yet, propagate delta to J0
361: if (!same && !lmvm->prev_set && (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL)) {
362: ldb->initialized = PETSC_FALSE;
363: B->preallocated = PETSC_FALSE; // make sure SyBroydenInitializeJ0() is called in the next MatCheckPreallocated()
364: }
365: PetscFunctionReturn(PETSC_SUCCESS);
366: }
368: PETSC_INTERN PetscErrorCode SymBroydenRescaleCopy(SymBroydenRescale bctx, SymBroydenRescale mctx)
369: {
370: PetscFunctionBegin;
371: mctx->scale_type = bctx->scale_type;
372: mctx->theta = bctx->theta;
373: mctx->alpha = bctx->alpha;
374: mctx->beta = bctx->beta;
375: mctx->rho = bctx->rho;
376: mctx->delta = bctx->delta;
377: mctx->delta_min = bctx->delta_min;
378: mctx->delta_max = bctx->delta_max;
379: mctx->tol = bctx->tol;
380: mctx->sigma_hist = bctx->sigma_hist;
381: mctx->forward = bctx->forward;
382: PetscFunctionReturn(PETSC_SUCCESS);
383: }
385: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetDiagonalMode(SymBroydenRescale ldb, PetscBool forward)
386: {
387: PetscFunctionBegin;
388: ldb->forward = forward;
389: PetscFunctionReturn(PETSC_SUCCESS);
390: }
392: PETSC_INTERN PetscErrorCode SymBroydenRescaleGetType(SymBroydenRescale ldb, MatLMVMSymBroydenScaleType *stype)
393: {
394: PetscFunctionBegin;
395: *stype = ldb->scale_type;
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetType(SymBroydenRescale ldb, MatLMVMSymBroydenScaleType stype)
400: {
401: PetscFunctionBegin;
402: ldb->scale_type = stype;
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetFromOptions(Mat B, SymBroydenRescale ldb, PetscOptionItems PetscOptionsObject)
407: {
408: MatLMVMSymBroydenScaleType stype = ldb->scale_type;
409: PetscBool flg;
411: PetscFunctionBegin;
412: PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for updating diagonal Jacobian approximation (MATLMVMDIAGBRDN)");
413: PetscCall(PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBroydenScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg));
414: PetscCall(PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL));
415: PetscCall(PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL));
416: PetscCall(PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL));
417: PetscCall(PetscOptionsRangeReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL, 0.0, 1.0));
418: PetscCall(PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL));
419: PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL));
420: PetscCall(PetscOptionsBoundedInt("-mat_lmvm_sigma_hist", "(developer) number of past updates to use in the default J0 scalar", "", ldb->sigma_hist, &ldb->sigma_hist, NULL, 0));
421: PetscOptionsHeadEnd();
422: PetscCheck(!(ldb->theta < 0.0) && !(ldb->theta > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
423: PetscCheck(!(ldb->alpha < 0.0) && !(ldb->alpha > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
424: PetscCheck(!(ldb->rho < 0.0) && !(ldb->rho > 1.0), PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex update limiter in the J0 scaling cannot be outside the range of [0, 1]");
425: PetscCheck(ldb->sigma_hist >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
426: if (flg) PetscCall(SymBroydenRescaleSetType(ldb, stype));
427: PetscFunctionReturn(PETSC_SUCCESS);
428: }
430: static PetscErrorCode SymBroydenRescaleAllocate(Mat B, SymBroydenRescale ldb)
431: {
432: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
434: PetscFunctionBegin;
435: if (!ldb->allocated) {
436: PetscCall(PetscMalloc3(ldb->sigma_hist, &ldb->yty, ldb->sigma_hist, &ldb->yts, ldb->sigma_hist, &ldb->sts));
437: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew));
438: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS));
439: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP));
440: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U));
441: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V));
442: PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W));
443: ldb->allocated = PETSC_TRUE;
444: }
445: PetscFunctionReturn(PETSC_SUCCESS);
446: }
448: PETSC_INTERN PetscErrorCode SymBroydenRescaleSetUp(Mat B, SymBroydenRescale ldb)
449: {
450: PetscFunctionBegin;
451: PetscCall(SymBroydenRescaleAllocate(B, ldb));
452: if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DECIDE) {
453: Mat J0;
454: PetscBool is_constant_or_diagonal;
456: PetscCall(MatLMVMGetJ0(B, &J0));
457: PetscCall(PetscObjectTypeCompareAny((PetscObject)J0, &is_constant_or_diagonal, MATCONSTANTDIAGONAL, MATDIAGONAL, ""));
458: ldb->scale_type = is_constant_or_diagonal ? MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL : MAT_LMVM_SYMBROYDEN_SCALE_NONE;
459: }
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: PETSC_INTERN PetscErrorCode SymBroydenRescaleInitializeJ0(Mat B, SymBroydenRescale ldb)
464: {
465: PetscFunctionBegin;
466: if (!ldb->initialized) {
467: PetscCall(SymBroydenRescaleSetUp(B, ldb));
468: switch (ldb->scale_type) {
469: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL: {
470: Vec invD;
472: PetscCall(MatLMVMGetJ0InvDiag(B, &invD));
473: PetscCall(VecSet(invD, ldb->delta));
474: PetscCall(MatLMVMRestoreJ0InvDiag(B, &invD));
475: break;
476: }
477: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
478: PetscCall(MatLMVMSetJ0Scale(B, 1.0 / ldb->delta));
479: break;
480: default:
481: break;
482: }
483: ldb->initialized = PETSC_TRUE;
484: }
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: PETSC_INTERN PetscErrorCode SymBroydenRescaleView(SymBroydenRescale ldb, PetscViewer pv)
489: {
490: PetscFunctionBegin;
491: PetscBool isascii;
492: PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii));
493: if (isascii) {
494: PetscCall(PetscViewerASCIIPrintf(pv, "Rescale type: %s\n", MatLMVMSymBroydenScaleTypes[ldb->scale_type]));
495: if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
496: PetscCall(PetscViewerASCIIPrintf(pv, "Rescale history: %" PetscInt_FMT "\n", ldb->sigma_hist));
497: PetscCall(PetscViewerASCIIPrintf(pv, "Rescale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho));
498: }
499: if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(PetscViewerASCIIPrintf(pv, "Rescale convex factor: theta=%g\n", (double)ldb->theta));
500: }
501: PetscFunctionReturn(PETSC_SUCCESS);
502: }
504: PETSC_INTERN PetscErrorCode SymBroydenRescaleReset(Mat B, SymBroydenRescale ldb, MatLMVMResetMode mode)
505: {
506: PetscFunctionBegin;
507: ldb->k = 0;
508: if (MatLMVMResetClearsBases(mode)) {
509: if (ldb->allocated) {
510: PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts));
511: PetscCall(VecDestroy(&ldb->invDnew));
512: PetscCall(VecDestroy(&ldb->BFGS));
513: PetscCall(VecDestroy(&ldb->DFP));
514: PetscCall(VecDestroy(&ldb->U));
515: PetscCall(VecDestroy(&ldb->V));
516: PetscCall(VecDestroy(&ldb->W));
517: ldb->allocated = PETSC_FALSE;
518: }
519: }
520: if (B && ldb->initialized && !MatLMVMResetClearsAll(mode)) PetscCall(SymBroydenRescaleInitializeJ0(B, ldb)); // eagerly reset J0 if we are rescaling
521: PetscFunctionReturn(PETSC_SUCCESS);
522: }
524: PETSC_INTERN PetscErrorCode SymBroydenRescaleDestroy(SymBroydenRescale *ldb)
525: {
526: PetscFunctionBegin;
527: PetscCall(SymBroydenRescaleReset(NULL, *ldb, MAT_LMVM_RESET_ALL));
528: PetscCall(PetscFree(*ldb));
529: *ldb = NULL;
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: PETSC_INTERN PetscErrorCode SymBroydenRescaleCreate(SymBroydenRescale *ldb)
534: {
535: PetscFunctionBegin;
536: PetscCall(PetscNew(ldb));
537: (*ldb)->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_DECIDE;
538: (*ldb)->theta = 0.0;
539: (*ldb)->alpha = 1.0;
540: (*ldb)->rho = 1.0;
541: (*ldb)->forward = PETSC_TRUE;
542: (*ldb)->beta = 0.5;
543: (*ldb)->delta = 1.0;
544: (*ldb)->delta_min = 1e-7;
545: (*ldb)->delta_max = 100.0;
546: (*ldb)->tol = 1e-8;
547: (*ldb)->sigma_hist = 1;
548: (*ldb)->allocated = PETSC_FALSE;
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }