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 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 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 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) { 302 PetscCall(VecCopy(ldb->invDnew, invD)); 303 } else if (ldb->rho) PetscCall(VecAXPBY(invD, 1.0 - ldb->rho, ldb->rho, ldb->invDnew)); 304 PetscCall(MatLMVMRestoreJ0InvDiag(B, &invD)); 305 PetscFunctionReturn(PETSC_SUCCESS); 306 } 307 308 static PetscErrorCode SymBroydenRescaleUpdateJ0(Mat B, SymBroydenRescale ldb) 309 { 310 PetscFunctionBegin; 311 if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(SymBroydenRescaleUpdateScalar(B, ldb)); 312 else if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(SymBroydenRescaleUpdateDiagonal(B, ldb)); 313 PetscFunctionReturn(PETSC_SUCCESS); 314 } 315 316 PETSC_INTERN PetscErrorCode SymBroydenRescaleUpdate(Mat B, SymBroydenRescale ldb) 317 { 318 PetscInt oldest, next; 319 320 PetscFunctionBegin; 321 PetscCall(SymBroydenRescaleSetUp(B, ldb)); 322 if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_NONE || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_USER) PetscFunctionReturn(PETSC_SUCCESS); 323 PetscCall(PetscLogEventBegin(SBRDN_Rescale, NULL, NULL, NULL, NULL)); 324 PetscCall(MatLMVMGetRange(B, &oldest, &next)); 325 if (next > ldb->k) { 326 PetscInt new_oldest = PetscMax(0, next - ldb->sigma_hist); 327 PetscInt ldb_oldest = PetscMax(0, ldb->k - ldb->sigma_hist); 328 329 if (new_oldest > ldb_oldest) { 330 for (PetscInt i = new_oldest; i < ldb->k; i++) { 331 ldb->yty[i - new_oldest] = ldb->yty[i - ldb_oldest]; 332 ldb->yts[i - new_oldest] = ldb->yts[i - ldb_oldest]; 333 ldb->sts[i - new_oldest] = ldb->sts[i - ldb_oldest]; 334 } 335 } 336 for (PetscInt i = PetscMax(new_oldest, ldb->k); i < next; i++) { 337 PetscScalar yty, sts, yts; 338 339 PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_Y, i, &yty)); 340 PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_Y, LMBASIS_S, i, &yts)); 341 PetscCall(MatLMVMProductsGetDiagonalValue(B, LMBASIS_S, LMBASIS_S, i, &sts)); 342 ldb->yty[i - new_oldest] = PetscRealPart(yty); 343 ldb->yts[i - new_oldest] = PetscRealPart(yts); 344 ldb->sts[i - new_oldest] = PetscRealPart(sts); 345 } 346 ldb->k = next; 347 } 348 PetscCall(SymBroydenRescaleUpdateJ0(B, ldb)); 349 PetscCall(PetscLogEventEnd(SBRDN_Rescale, NULL, NULL, NULL, NULL)); 350 PetscFunctionReturn(PETSC_SUCCESS); 351 } 352 353 PETSC_INTERN PetscErrorCode SymBroydenRescaleSetDelta(Mat B, SymBroydenRescale ldb, PetscReal delta) 354 { 355 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 356 PetscBool same; 357 358 PetscFunctionBegin; 359 same = (delta == ldb->delta) ? PETSC_TRUE : PETSC_FALSE; 360 ldb->delta = delta; 361 ldb->delta = PetscMin(ldb->delta, ldb->delta_max); 362 ldb->delta = PetscMax(ldb->delta, ldb->delta_min); 363 // if there have been no updates yet, propagate delta to J0 364 if (!same && !lmvm->prev_set && (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL)) { 365 ldb->initialized = PETSC_FALSE; 366 B->preallocated = PETSC_FALSE; // make sure SyBroydenInitializeJ0() is called in the next MatCheckPreallocated() 367 } 368 PetscFunctionReturn(PETSC_SUCCESS); 369 } 370 371 PETSC_INTERN PetscErrorCode SymBroydenRescaleCopy(SymBroydenRescale bctx, SymBroydenRescale mctx) 372 { 373 PetscFunctionBegin; 374 mctx->scale_type = bctx->scale_type; 375 mctx->theta = bctx->theta; 376 mctx->alpha = bctx->alpha; 377 mctx->beta = bctx->beta; 378 mctx->rho = bctx->rho; 379 mctx->delta = bctx->delta; 380 mctx->delta_min = bctx->delta_min; 381 mctx->delta_max = bctx->delta_max; 382 mctx->tol = bctx->tol; 383 mctx->sigma_hist = bctx->sigma_hist; 384 mctx->forward = bctx->forward; 385 PetscFunctionReturn(PETSC_SUCCESS); 386 } 387 388 PETSC_INTERN PetscErrorCode SymBroydenRescaleSetDiagonalMode(SymBroydenRescale ldb, PetscBool forward) 389 { 390 PetscFunctionBegin; 391 ldb->forward = forward; 392 PetscFunctionReturn(PETSC_SUCCESS); 393 } 394 395 PETSC_INTERN PetscErrorCode SymBroydenRescaleGetType(SymBroydenRescale ldb, MatLMVMSymBroydenScaleType *stype) 396 { 397 PetscFunctionBegin; 398 *stype = ldb->scale_type; 399 PetscFunctionReturn(PETSC_SUCCESS); 400 } 401 402 PETSC_INTERN PetscErrorCode SymBroydenRescaleSetType(SymBroydenRescale ldb, MatLMVMSymBroydenScaleType stype) 403 { 404 PetscFunctionBegin; 405 ldb->scale_type = stype; 406 PetscFunctionReturn(PETSC_SUCCESS); 407 } 408 409 PETSC_INTERN PetscErrorCode SymBroydenRescaleSetFromOptions(Mat B, SymBroydenRescale ldb, PetscOptionItems PetscOptionsObject) 410 { 411 MatLMVMSymBroydenScaleType stype = ldb->scale_type; 412 PetscBool flg; 413 414 PetscFunctionBegin; 415 PetscOptionsHeadBegin(PetscOptionsObject, "Restricted Broyden method for updating diagonal Jacobian approximation (MATLMVMDIAGBRDN)"); 416 PetscCall(PetscOptionsEnum("-mat_lmvm_scale_type", "(developer) scaling type applied to J0", "MatLMVMSymBroydenScaleType", MatLMVMSymBroydenScaleTypes, (PetscEnum)stype, (PetscEnum *)&stype, &flg)); 417 PetscCall(PetscOptionsReal("-mat_lmvm_theta", "(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling", "", ldb->theta, &ldb->theta, NULL)); 418 PetscCall(PetscOptionsReal("-mat_lmvm_rho", "(developer) update limiter in the J0 scaling", "", ldb->rho, &ldb->rho, NULL)); 419 PetscCall(PetscOptionsReal("-mat_lmvm_tol", "(developer) tolerance for bounding rescaling denominator", "", ldb->tol, &ldb->tol, NULL)); 420 PetscCall(PetscOptionsRangeReal("-mat_lmvm_alpha", "(developer) convex ratio in the J0 scaling", "", ldb->alpha, &ldb->alpha, NULL, 0.0, 1.0)); 421 PetscCall(PetscOptionsBool("-mat_lmvm_forward", "Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.", "", ldb->forward, &ldb->forward, NULL)); 422 PetscCall(PetscOptionsReal("-mat_lmvm_beta", "(developer) exponential factor in the diagonal J0 scaling", "", ldb->beta, &ldb->beta, NULL)); 423 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)); 424 PetscOptionsHeadEnd(); 425 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]"); 426 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]"); 427 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]"); 428 PetscCheck(ldb->sigma_hist >= 0, PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative"); 429 if (flg) PetscCall(SymBroydenRescaleSetType(ldb, stype)); 430 PetscFunctionReturn(PETSC_SUCCESS); 431 } 432 433 static PetscErrorCode SymBroydenRescaleAllocate(Mat B, SymBroydenRescale ldb) 434 { 435 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 436 437 PetscFunctionBegin; 438 if (!ldb->allocated) { 439 PetscCall(PetscMalloc3(ldb->sigma_hist, &ldb->yty, ldb->sigma_hist, &ldb->yts, ldb->sigma_hist, &ldb->sts)); 440 PetscCall(VecDuplicate(lmvm->Xprev, &ldb->invDnew)); 441 PetscCall(VecDuplicate(lmvm->Xprev, &ldb->BFGS)); 442 PetscCall(VecDuplicate(lmvm->Xprev, &ldb->DFP)); 443 PetscCall(VecDuplicate(lmvm->Xprev, &ldb->U)); 444 PetscCall(VecDuplicate(lmvm->Xprev, &ldb->V)); 445 PetscCall(VecDuplicate(lmvm->Xprev, &ldb->W)); 446 ldb->allocated = PETSC_TRUE; 447 } 448 PetscFunctionReturn(PETSC_SUCCESS); 449 } 450 451 PETSC_INTERN PetscErrorCode SymBroydenRescaleSetUp(Mat B, SymBroydenRescale ldb) 452 { 453 PetscFunctionBegin; 454 PetscCall(SymBroydenRescaleAllocate(B, ldb)); 455 if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DECIDE) { 456 Mat J0; 457 PetscBool is_constant_or_diagonal; 458 459 PetscCall(MatLMVMGetJ0(B, &J0)); 460 PetscCall(PetscObjectTypeCompareAny((PetscObject)J0, &is_constant_or_diagonal, MATCONSTANTDIAGONAL, MATDIAGONAL, "")); 461 ldb->scale_type = is_constant_or_diagonal ? MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL : MAT_LMVM_SYMBROYDEN_SCALE_NONE; 462 } 463 PetscFunctionReturn(PETSC_SUCCESS); 464 } 465 466 PETSC_INTERN PetscErrorCode SymBroydenRescaleInitializeJ0(Mat B, SymBroydenRescale ldb) 467 { 468 PetscFunctionBegin; 469 if (!ldb->initialized) { 470 PetscCall(SymBroydenRescaleSetUp(B, ldb)); 471 switch (ldb->scale_type) { 472 case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL: { 473 Vec invD; 474 475 PetscCall(MatLMVMGetJ0InvDiag(B, &invD)); 476 PetscCall(VecSet(invD, ldb->delta)); 477 PetscCall(MatLMVMRestoreJ0InvDiag(B, &invD)); 478 break; 479 } 480 case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR: 481 PetscCall(MatLMVMSetJ0Scale(B, 1.0 / ldb->delta)); 482 break; 483 default: 484 break; 485 } 486 ldb->initialized = PETSC_TRUE; 487 } 488 PetscFunctionReturn(PETSC_SUCCESS); 489 } 490 491 PETSC_INTERN PetscErrorCode SymBroydenRescaleView(SymBroydenRescale ldb, PetscViewer pv) 492 { 493 PetscFunctionBegin; 494 PetscBool isascii; 495 PetscCall(PetscObjectTypeCompare((PetscObject)pv, PETSCVIEWERASCII, &isascii)); 496 if (isascii) { 497 PetscCall(PetscViewerASCIIPrintf(pv, "Rescale type: %s\n", MatLMVMSymBroydenScaleTypes[ldb->scale_type])); 498 if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR || ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) { 499 PetscCall(PetscViewerASCIIPrintf(pv, "Rescale history: %" PetscInt_FMT "\n", ldb->sigma_hist)); 500 PetscCall(PetscViewerASCIIPrintf(pv, "Rescale params: alpha=%g, beta=%g, rho=%g\n", (double)ldb->alpha, (double)ldb->beta, (double)ldb->rho)); 501 } 502 if (ldb->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(PetscViewerASCIIPrintf(pv, "Rescale convex factor: theta=%g\n", (double)ldb->theta)); 503 } 504 PetscFunctionReturn(PETSC_SUCCESS); 505 } 506 507 PETSC_INTERN PetscErrorCode SymBroydenRescaleReset(Mat B, SymBroydenRescale ldb, MatLMVMResetMode mode) 508 { 509 PetscFunctionBegin; 510 ldb->k = 0; 511 if (MatLMVMResetClearsBases(mode)) { 512 if (ldb->allocated) { 513 PetscCall(PetscFree3(ldb->yty, ldb->yts, ldb->sts)); 514 PetscCall(VecDestroy(&ldb->invDnew)); 515 PetscCall(VecDestroy(&ldb->BFGS)); 516 PetscCall(VecDestroy(&ldb->DFP)); 517 PetscCall(VecDestroy(&ldb->U)); 518 PetscCall(VecDestroy(&ldb->V)); 519 PetscCall(VecDestroy(&ldb->W)); 520 ldb->allocated = PETSC_FALSE; 521 } 522 } 523 if (B && ldb->initialized && !MatLMVMResetClearsAll(mode)) PetscCall(SymBroydenRescaleInitializeJ0(B, ldb)); // eagerly reset J0 if we are rescaling 524 PetscFunctionReturn(PETSC_SUCCESS); 525 } 526 527 PETSC_INTERN PetscErrorCode SymBroydenRescaleDestroy(SymBroydenRescale *ldb) 528 { 529 PetscFunctionBegin; 530 PetscCall(SymBroydenRescaleReset(NULL, *ldb, MAT_LMVM_RESET_ALL)); 531 PetscCall(PetscFree(*ldb)); 532 *ldb = NULL; 533 PetscFunctionReturn(PETSC_SUCCESS); 534 } 535 536 PETSC_INTERN PetscErrorCode SymBroydenRescaleCreate(SymBroydenRescale *ldb) 537 { 538 PetscFunctionBegin; 539 PetscCall(PetscNew(ldb)); 540 (*ldb)->scale_type = MAT_LMVM_SYMBROYDEN_SCALE_DECIDE; 541 (*ldb)->theta = 0.0; 542 (*ldb)->alpha = 1.0; 543 (*ldb)->rho = 1.0; 544 (*ldb)->forward = PETSC_TRUE; 545 (*ldb)->beta = 0.5; 546 (*ldb)->delta = 1.0; 547 (*ldb)->delta_min = 1e-7; 548 (*ldb)->delta_max = 100.0; 549 (*ldb)->tol = 1e-8; 550 (*ldb)->sigma_hist = 1; 551 (*ldb)->allocated = PETSC_FALSE; 552 PetscFunctionReturn(PETSC_SUCCESS); 553 } 554