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