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) 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 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 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 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 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 387 PETSC_INTERN PetscErrorCode SymBroydenRescaleSetDiagonalMode(SymBroydenRescale ldb, PetscBool forward) 388 { 389 PetscFunctionBegin; 390 ldb->forward = forward; 391 PetscFunctionReturn(PETSC_SUCCESS); 392 } 393 394 PETSC_INTERN PetscErrorCode SymBroydenRescaleGetType(SymBroydenRescale ldb, MatLMVMSymBroydenScaleType *stype) 395 { 396 PetscFunctionBegin; 397 *stype = ldb->scale_type; 398 PetscFunctionReturn(PETSC_SUCCESS); 399 } 400 401 PETSC_INTERN PetscErrorCode SymBroydenRescaleSetType(SymBroydenRescale ldb, MatLMVMSymBroydenScaleType stype) 402 { 403 PetscFunctionBegin; 404 ldb->scale_type = stype; 405 PetscFunctionReturn(PETSC_SUCCESS); 406 } 407 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 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 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 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 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 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 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 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