1 #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h> /*I "petscksp.h" I*/ 2 3 // Bad Broyden is dual to Broyden: MatSolve routines are dual to Broyden MatMult routines 4 5 static PetscErrorCode MatSolve_LMVMBadBrdn_Recursive(Mat B, Vec F, Vec dX) 6 { 7 PetscFunctionBegin; 8 PetscCall(BroydenKernel_Recursive(B, MATLMVM_MODE_DUAL, F, dX)); 9 PetscFunctionReturn(PETSC_SUCCESS); 10 } 11 12 static PetscErrorCode MatSolve_LMVMBadBrdn_CompactDense(Mat B, Vec F, Vec dX) 13 { 14 PetscFunctionBegin; 15 PetscCall(BroydenKernel_CompactDense(B, MATLMVM_MODE_DUAL, F, dX)); 16 PetscFunctionReturn(PETSC_SUCCESS); 17 } 18 19 PETSC_UNUSED static PetscErrorCode MatSolve_LMVMBadBrdn_Dense(Mat B, Vec F, Vec dX) 20 { 21 PetscFunctionBegin; 22 PetscCall(BroydenKernel_Dense(B, MATLMVM_MODE_DUAL, F, dX)); 23 PetscFunctionReturn(PETSC_SUCCESS); 24 } 25 26 static PetscErrorCode MatSolveHermitianTranspose_LMVMBadBrdn_Recursive(Mat B, Vec F, Vec dX) 27 { 28 PetscFunctionBegin; 29 PetscCall(BroydenKernelHermitianTranspose_Recursive(B, MATLMVM_MODE_DUAL, F, dX)); 30 PetscFunctionReturn(PETSC_SUCCESS); 31 } 32 33 static PetscErrorCode MatSolveHermitianTranspose_LMVMBadBrdn_CompactDense(Mat B, Vec F, Vec dX) 34 { 35 PetscFunctionBegin; 36 PetscCall(BroydenKernelHermitianTranspose_CompactDense(B, MATLMVM_MODE_DUAL, F, dX)); 37 PetscFunctionReturn(PETSC_SUCCESS); 38 } 39 40 PETSC_UNUSED static PetscErrorCode MatSolveHermitianTranspose_LMVMBadBrdn_Dense(Mat B, Vec F, Vec dX) 41 { 42 PetscFunctionBegin; 43 PetscCall(BroydenKernelHermitianTranspose_Dense(B, MATLMVM_MODE_DUAL, F, dX)); 44 PetscFunctionReturn(PETSC_SUCCESS); 45 } 46 47 /* 48 The bad Broyden kernel can be written as 49 50 B_{k+1} x = B_k x + (y_k - B_k s_k)^T * (y_k^T B_k s_k)^{-1} y_k^T B_k x 51 = (I + (y_k - B_k s_k)^T (y_k^T B_k s_k)^{-1} y_k^T) (B_k x) 52 \______________________ _________________________/ 53 V 54 recursive rank-1 update 55 56 When the basis (y_k - B_k s_k) and the products (y_k^T B_k s_k) have been computed, the product can be computed by 57 application of rank-1 updates from oldest to newest 58 */ 59 60 static PetscErrorCode BadBroydenKernel_Recursive_Inner(Mat B, MatLMVMMode mode, PetscInt oldest, PetscInt next, Vec B0X) 61 { 62 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 63 Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx; 64 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode); 65 LMBasis Y_minus_BkS = lbrdn->basis[LMVMModeMap(BROYDEN_BASIS_Y_MINUS_BKS, mode)]; 66 LMBasis Y; 67 LMProducts YtBkS = lbrdn->products[LMVMModeMap(BROYDEN_PRODUCTS_YTBKS, mode)]; 68 69 PetscFunctionBegin; 70 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL)); 71 // These cannot be parallelized, notice the data dependence 72 for (PetscInt i = oldest; i < next; i++) { 73 Vec y_i, yimbisi; 74 PetscScalar yitbix; 75 PetscScalar yitbisi; 76 77 PetscCall(LMBasisGetVecRead(Y, i, &y_i)); 78 PetscCall(VecDot(B0X, y_i, &yitbix)); 79 PetscCall(LMBasisRestoreVecRead(Y, i, &y_i)); 80 PetscCall(LMProductsGetDiagonalValue(YtBkS, i, &yitbisi)); 81 PetscCall(LMBasisGetVecRead(Y_minus_BkS, i, &yimbisi)); 82 PetscCall(VecAXPY(B0X, yitbix / yitbisi, yimbisi)); 83 PetscCall(LMBasisRestoreVecRead(Y_minus_BkS, i, &yimbisi)); 84 } 85 PetscFunctionReturn(PETSC_SUCCESS); 86 } 87 88 /* 89 Compute the basis vectors (y_k - B_k s_k) and dot products (y_k^T B_k s_k) recursively 90 */ 91 92 static PetscErrorCode BadBroydenRecursiveBasisUpdate(Mat B, MatLMVMMode mode) 93 { 94 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 95 Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx; 96 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode); 97 MatLMVMBasisType B0S_t = LMVMModeMap(LMBASIS_B0S, mode); 98 LMBasis Y_minus_BkS; 99 LMProducts YtBkS; 100 BroydenBasisType Y_minus_BkS_t = LMVMModeMap(BROYDEN_BASIS_Y_MINUS_BKS, mode); 101 BroydenProductsType YtBkS_t = LMVMModeMap(BROYDEN_PRODUCTS_YTBKS, mode); 102 PetscInt oldest, next; 103 PetscInt products_oldest; 104 LMBasis Y; 105 PetscInt start; 106 107 PetscFunctionBegin; 108 if (!lbrdn->basis[Y_minus_BkS_t]) PetscCall(LMBasisCreate(Y_t == LMBASIS_Y ? lmvm->Fprev : lmvm->Xprev, lmvm->m, &lbrdn->basis[Y_minus_BkS_t])); 109 Y_minus_BkS = lbrdn->basis[Y_minus_BkS_t]; 110 if (!lbrdn->products[YtBkS_t]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_DIAGONAL, &lbrdn->products[YtBkS_t])); 111 YtBkS = lbrdn->products[YtBkS_t]; 112 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL)); 113 PetscCall(MatLMVMGetRange(B, &oldest, &next)); 114 // invalidate computed values if J0 has changed 115 PetscCall(LMProductsPrepare(YtBkS, lmvm->J0, oldest, next)); 116 products_oldest = PetscMax(0, YtBkS->k - lmvm->m); 117 if (oldest > products_oldest) { 118 // recursion is starting from a different starting index, it must be recomputed 119 YtBkS->k = oldest; 120 } 121 Y_minus_BkS->k = start = YtBkS->k; 122 // recompute each column vector in Y_minus_BkS, and the product y_k^T B_k s_k, in order 123 for (PetscInt j = start; j < next; j++) { 124 Vec p_j, y_j, B0s_j; 125 PetscScalar yjtbjsj, alpha; 126 127 PetscCall(LMBasisGetWorkVec(Y_minus_BkS, &p_j)); 128 129 // p_j starts as B_0 * s_j 130 PetscCall(MatLMVMBasisGetVecRead(B, B0S_t, j, &B0s_j, &alpha)); 131 PetscCall(VecAXPBY(p_j, alpha, 0.0, B0s_j)); 132 PetscCall(MatLMVMBasisRestoreVecRead(B, B0S_t, j, &B0s_j, &alpha)); 133 134 /* Use the matmult kernel to compute p_j = B_j * p_j 135 (if j == oldest, p_j is already correct) */ 136 if (j > oldest) PetscCall(BadBroydenKernel_Recursive_Inner(B, mode, oldest, j, p_j)); 137 PetscCall(LMBasisGetVecRead(Y, j, &y_j)); 138 PetscCall(VecDot(p_j, y_j, &yjtbjsj)); 139 PetscCall(LMProductsInsertNextDiagonalValue(YtBkS, j, yjtbjsj)); 140 PetscCall(VecAYPX(p_j, -1.0, y_j)); 141 PetscCall(LMBasisRestoreVecRead(Y, j, &y_j)); 142 PetscCall(LMBasisSetNextVec(Y_minus_BkS, p_j)); 143 PetscCall(LMBasisRestoreWorkVec(Y_minus_BkS, &p_j)); 144 } 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 148 PETSC_INTERN PetscErrorCode BadBroydenKernel_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec Y) 149 { 150 PetscInt oldest, next; 151 152 PetscFunctionBegin; 153 PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, Y)); 154 PetscCall(MatLMVMGetRange(B, &oldest, &next)); 155 if (next > oldest) { 156 PetscCall(BadBroydenRecursiveBasisUpdate(B, mode)); 157 PetscCall(BadBroydenKernel_Recursive_Inner(B, mode, oldest, next, Y)); 158 } 159 PetscFunctionReturn(PETSC_SUCCESS); 160 } 161 162 /* 163 The adjoint of the recursive bad Broyden kernel is 164 165 B_{k+1}^T x = B_k^T (I + y_k (s_k^T B_k^T y_k)^{-1} (y_k - B_k s_k)^T) x 166 \____________________ ___________________________/ 167 V 168 recursive rank-1 update 169 170 which can be computed by application of rank-1 updates from newest to oldest 171 */ 172 173 static PetscErrorCode BadBroydenKernelHermitianTranspose_Recursive_Inner(Mat B, MatLMVMMode mode, PetscInt oldest, PetscInt next, Vec X) 174 { 175 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode); 176 BroydenBasisType Y_minus_BkS_t = LMVMModeMap(BROYDEN_BASIS_Y_MINUS_BKS, mode); 177 BroydenProductsType YtBkS_t = LMVMModeMap(BROYDEN_PRODUCTS_YTBKS, mode); 178 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 179 Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx; 180 LMBasis Y_minus_BkS = lbrdn->basis[Y_minus_BkS_t]; 181 LMBasis Y; 182 LMProducts YtBkS = lbrdn->products[YtBkS_t]; 183 184 PetscFunctionBegin; 185 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL)); 186 // These cannot be parallelized, notice the data dependence 187 for (PetscInt i = next - 1; i >= oldest; i--) { 188 Vec yimbisi, y_i; 189 PetscScalar yimbisitx; 190 PetscScalar yitbisi; 191 192 PetscCall(LMBasisGetVecRead(Y_minus_BkS, i, &yimbisi)); 193 PetscCall(VecDot(X, yimbisi, &yimbisitx)); 194 PetscCall(LMBasisRestoreVecRead(Y_minus_BkS, i, &yimbisi)); 195 PetscCall(LMProductsGetDiagonalValue(YtBkS, i, &yitbisi)); 196 PetscCall(LMBasisGetVecRead(Y, i, &y_i)); 197 PetscCall(VecAXPY(X, yimbisitx / PetscConj(yitbisi), y_i)); 198 PetscCall(LMBasisRestoreVecRead(Y, i, &y_i)); 199 } 200 PetscFunctionReturn(PETSC_SUCCESS); 201 } 202 203 PETSC_INTERN PetscErrorCode BadBroydenKernelHermitianTranspose_Recursive(Mat B, MatLMVMMode mode, Vec X, Vec BX) 204 { 205 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode); 206 PetscInt oldest, next; 207 Vec G = X; 208 LMBasis Y; 209 210 PetscFunctionBegin; 211 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL)); 212 PetscCall(MatLMVMGetRange(B, &oldest, &next)); 213 if (next > oldest) { 214 PetscCall(LMBasisGetWorkVec(Y, &G)); 215 PetscCall(VecCopy(X, G)); 216 PetscCall(BadBroydenRecursiveBasisUpdate(B, mode)); 217 PetscCall(BadBroydenKernelHermitianTranspose_Recursive_Inner(B, mode, oldest, next, G)); 218 } 219 PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, G, BX)); 220 if (next > oldest) PetscCall(LMBasisRestoreWorkVec(Y, &G)); 221 PetscFunctionReturn(PETSC_SUCCESS); 222 } 223 224 /* 225 The bad Broyden kernel can be written as 226 227 B_k = B_0 + (Y - B_0 S) (Y^H B_0 S - stril(Y^H Y))^{-1} Y^H B_0 228 = (I + (Y - B_0 S) (Y^H B_0 S - stril(Y^H Y))^{-1} Y^H) B_0 229 230 where stril is the strictly lower triangular component. We compute and factorize 231 the small matrix in order to apply a single rank m update 232 */ 233 234 static PetscErrorCode BadBroydenCompactProductsUpdate(Mat B, MatLMVMMode mode) 235 { 236 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode); 237 MatLMVMBasisType B0S_t = LMVMModeMap(LMBASIS_B0S, mode); 238 BroydenProductsType YtB0S_minus_YtY_t = LMVMModeMap(BROYDEN_PRODUCTS_YTB0S_MINUS_YTY, mode); 239 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 240 Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx; 241 LMProducts YtB0S, YtY, YtB0S_minus_YtY; 242 LMBasis Y, B0S; 243 PetscScalar alpha; 244 PetscInt oldest, k, next; 245 PetscBool local_is_nonempty; 246 Mat ytbsmyty_local; 247 248 PetscFunctionBegin; 249 if (!lbrdn->products[YtB0S_minus_YtY_t]) PetscCall(MatLMVMCreateProducts(B, LMBLOCK_FULL, &lbrdn->products[YtB0S_minus_YtY_t])); 250 YtB0S_minus_YtY = lbrdn->products[YtB0S_minus_YtY_t]; 251 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL)); 252 PetscCall(MatLMVMGetUpdatedBasis(B, B0S_t, &B0S, &B0S_t, &alpha)); 253 PetscCall(MatLMVMGetRange(B, &oldest, &next)); 254 // invalidate computed values if J0 has changed 255 PetscCall(LMProductsPrepare(YtB0S_minus_YtY, lmvm->J0, oldest, next)); 256 PetscCall(LMProductsGetLocalMatrix(YtB0S_minus_YtY, &ytbsmyty_local, &k, &local_is_nonempty)); 257 if (k < next) { 258 Mat yty_local, ytbs_local; 259 260 PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, B0S_t, LMBLOCK_FULL, &YtB0S)); 261 PetscCall(MatLMVMGetUpdatedProducts(B, Y_t, Y_t, LMBLOCK_STRICT_UPPER_TRIANGLE, &YtY)); 262 PetscCall(LMProductsGetLocalMatrix(YtB0S, &ytbs_local, NULL, NULL)); 263 PetscCall(LMProductsGetLocalMatrix(YtY, &yty_local, NULL, NULL)); 264 if (local_is_nonempty) { 265 PetscCall(MatSetUnfactored(ytbsmyty_local)); 266 PetscCall(MatCopy(yty_local, ytbsmyty_local, SAME_NONZERO_PATTERN)); 267 PetscCall(MatTranspose(ytbsmyty_local, MAT_INPLACE_MATRIX, &ytbsmyty_local)); 268 if (PetscDefined(USE_COMPLEX)) PetscCall(MatConjugate(ytbsmyty_local)); 269 PetscCall(MatScale(ytbsmyty_local, -1.0)); 270 PetscCall(MatAXPY(ytbsmyty_local, alpha, ytbs_local, SAME_NONZERO_PATTERN)); 271 PetscCall(LMProductsOnesOnUnusedDiagonal(ytbsmyty_local, oldest, next)); 272 PetscCall(MatLUFactor(ytbsmyty_local, NULL, NULL, NULL)); 273 } 274 PetscCall(LMProductsRestoreLocalMatrix(YtY, &yty_local, NULL)); 275 PetscCall(LMProductsRestoreLocalMatrix(YtB0S, &ytbs_local, NULL)); 276 } 277 PetscCall(LMProductsRestoreLocalMatrix(YtB0S_minus_YtY, &ytbsmyty_local, &next)); 278 PetscFunctionReturn(PETSC_SUCCESS); 279 } 280 281 PETSC_INTERN PetscErrorCode BadBroydenKernel_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BX) 282 { 283 PetscInt oldest, next; 284 285 PetscFunctionBegin; 286 PetscCall(MatLMVMApplyJ0Mode(mode)(B, X, BX)); 287 PetscCall(MatLMVMGetRange(B, &oldest, &next)); 288 if (next > oldest) { 289 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 290 Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx; 291 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode); 292 MatLMVMBasisType Y_minus_B0S_t = LMVMModeMap(LMBASIS_Y_MINUS_B0S, mode); 293 BroydenProductsType YtB0S_minus_YtY_t = LMVMModeMap(BROYDEN_PRODUCTS_YTB0S_MINUS_YTY, mode); 294 LMProducts YtB0S_minus_YtY; 295 LMBasis Y; 296 Vec YtB0X, v; 297 298 PetscCall(BadBroydenCompactProductsUpdate(B, mode)); 299 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL)); 300 YtB0S_minus_YtY = lbrdn->products[YtB0S_minus_YtY_t]; 301 PetscCall(MatLMVMGetWorkRow(B, &YtB0X)); 302 PetscCall(MatLMVMGetWorkRow(B, &v)); 303 PetscCall(LMBasisGEMVH(Y, oldest, next, 1.0, BX, 0.0, YtB0X)); 304 PetscCall(LMProductsSolve(YtB0S_minus_YtY, oldest, next, YtB0X, v, /* ^H */ PETSC_FALSE)); 305 PetscCall(MatLMVMBasisGEMV(B, Y_minus_B0S_t, oldest, next, 1.0, v, 1.0, BX)); 306 PetscCall(MatLMVMRestoreWorkRow(B, &v)); 307 PetscCall(MatLMVMRestoreWorkRow(B, &YtB0X)); 308 } 309 PetscFunctionReturn(PETSC_SUCCESS); 310 } 311 312 /* 313 The adjoint of the above formula for the bad Broyden kernel is 314 315 B_k^H = B_0^H + B_0^H Y (Y^H B_0 S - stril(Y^H Y))^{-H} (Y - B_0 S)^H 316 = B_0^H (I + Y (Y^H B_0 S - stril(Y^H Y))^{-H} (Y - B_0 S)^H) 317 */ 318 319 PETSC_INTERN PetscErrorCode BadBroydenKernelHermitianTranspose_CompactDense(Mat B, MatLMVMMode mode, Vec X, Vec BHX) 320 { 321 MatLMVMBasisType Y_t = LMVMModeMap(LMBASIS_Y, mode); 322 PetscInt oldest, next; 323 Vec G = X; 324 LMBasis Y = NULL; 325 326 PetscFunctionBegin; 327 PetscCall(MatLMVMGetRange(B, &oldest, &next)); 328 if (next > oldest) { 329 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 330 Mat_Brdn *lbrdn = (Mat_Brdn *)lmvm->ctx; 331 MatLMVMBasisType Y_minus_B0S_t = LMVMModeMap(LMBASIS_Y_MINUS_B0S, mode); 332 BroydenProductsType YtB0S_minus_YtY_t = LMVMModeMap(BROYDEN_PRODUCTS_YTB0S_MINUS_YTY, mode); 333 LMProducts YtB0S_minus_YtY; 334 Vec YmB0StG, v; 335 336 PetscCall(MatLMVMGetUpdatedBasis(B, Y_t, &Y, NULL, NULL)); 337 PetscCall(LMBasisGetWorkVec(Y, &G)); 338 PetscCall(VecCopy(X, G)); 339 PetscCall(BadBroydenCompactProductsUpdate(B, mode)); 340 YtB0S_minus_YtY = lbrdn->products[YtB0S_minus_YtY_t]; 341 PetscCall(MatLMVMGetWorkRow(B, &YmB0StG)); 342 PetscCall(MatLMVMGetWorkRow(B, &v)); 343 PetscCall(MatLMVMBasisGEMVH(B, Y_minus_B0S_t, oldest, next, 1.0, G, 0.0, YmB0StG)); 344 PetscCall(LMProductsSolve(YtB0S_minus_YtY, oldest, next, YmB0StG, v, /* ^H */ PETSC_TRUE)); 345 PetscCall(LMBasisGEMV(Y, oldest, next, 1.0, v, 1.0, G)); 346 PetscCall(MatLMVMRestoreWorkRow(B, &v)); 347 PetscCall(MatLMVMRestoreWorkRow(B, &YmB0StG)); 348 } 349 PetscCall(MatLMVMApplyJ0HermitianTransposeMode(mode)(B, G, BHX)); 350 if (next > oldest) PetscCall(LMBasisRestoreWorkVec(Y, &G)); 351 PetscFunctionReturn(PETSC_SUCCESS); 352 } 353 354 static PetscErrorCode MatMult_LMVMBadBrdn_Recursive(Mat B, Vec F, Vec dX) 355 { 356 PetscFunctionBegin; 357 PetscCall(BadBroydenKernel_Recursive(B, MATLMVM_MODE_PRIMAL, F, dX)); 358 PetscFunctionReturn(PETSC_SUCCESS); 359 } 360 361 static PetscErrorCode MatMult_LMVMBadBrdn_CompactDense(Mat B, Vec F, Vec dX) 362 { 363 PetscFunctionBegin; 364 PetscCall(BadBroydenKernel_CompactDense(B, MATLMVM_MODE_PRIMAL, F, dX)); 365 PetscFunctionReturn(PETSC_SUCCESS); 366 } 367 368 static PetscErrorCode MatMultHermitianTranspose_LMVMBadBrdn_Recursive(Mat B, Vec F, Vec dX) 369 { 370 PetscFunctionBegin; 371 PetscCall(BadBroydenKernelHermitianTranspose_Recursive(B, MATLMVM_MODE_PRIMAL, F, dX)); 372 PetscFunctionReturn(PETSC_SUCCESS); 373 } 374 375 static PetscErrorCode MatMultHermitianTranspose_LMVMBadBrdn_CompactDense(Mat B, Vec F, Vec dX) 376 { 377 PetscFunctionBegin; 378 PetscCall(BadBroydenKernelHermitianTranspose_CompactDense(B, MATLMVM_MODE_PRIMAL, F, dX)); 379 PetscFunctionReturn(PETSC_SUCCESS); 380 } 381 382 static PetscErrorCode MatLMVMSetMultAlgorithm_BadBrdn(Mat B) 383 { 384 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 385 386 PetscFunctionBegin; 387 switch (lmvm->mult_alg) { 388 case MAT_LMVM_MULT_RECURSIVE: 389 lmvm->ops->mult = MatMult_LMVMBadBrdn_Recursive; 390 lmvm->ops->multht = MatMultHermitianTranspose_LMVMBadBrdn_Recursive; 391 lmvm->ops->solve = MatSolve_LMVMBadBrdn_Recursive; 392 lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBadBrdn_Recursive; 393 break; 394 case MAT_LMVM_MULT_DENSE: 395 lmvm->ops->mult = MatMult_LMVMBadBrdn_CompactDense; 396 lmvm->ops->multht = MatMultHermitianTranspose_LMVMBadBrdn_CompactDense; 397 lmvm->ops->solve = MatSolve_LMVMBadBrdn_Dense; 398 lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBadBrdn_Dense; 399 break; 400 case MAT_LMVM_MULT_COMPACT_DENSE: 401 lmvm->ops->mult = MatMult_LMVMBadBrdn_CompactDense; 402 lmvm->ops->multht = MatMultHermitianTranspose_LMVMBadBrdn_CompactDense; 403 lmvm->ops->solve = MatSolve_LMVMBadBrdn_CompactDense; 404 lmvm->ops->solveht = MatSolveHermitianTranspose_LMVMBadBrdn_CompactDense; 405 break; 406 } 407 PetscFunctionReturn(PETSC_SUCCESS); 408 } 409 410 PetscErrorCode MatCreate_LMVMBadBrdn(Mat B) 411 { 412 Mat_LMVM *lmvm; 413 414 PetscFunctionBegin; 415 PetscCall(MatCreate_LMVMBrdn(B)); 416 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMBADBROYDEN)); 417 lmvm = (Mat_LMVM *)B->data; 418 lmvm->ops->setmultalgorithm = MatLMVMSetMultAlgorithm_BadBrdn; 419 lmvm->cache_gradient_products = PETSC_TRUE; 420 PetscCall(MatLMVMSetMultAlgorithm_BadBrdn(B)); 421 PetscFunctionReturn(PETSC_SUCCESS); 422 } 423 424 /*@ 425 MatCreateLMVMBadBroyden - Creates a limited-memory modified (aka "bad") Broyden-type 426 approximation matrix used for a Jacobian. L-BadBrdn is not guaranteed to be 427 symmetric or positive-definite. 428 429 To use the L-BadBrdn matrix with other vector types, the matrix must be 430 created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`. 431 This ensures that the internal storage and work vectors are duplicated from the 432 correct type of vector. 433 434 Collective 435 436 Input Parameters: 437 + comm - MPI communicator 438 . n - number of local rows for storage vectors 439 - N - global size of the storage vectors 440 441 Output Parameter: 442 . B - the matrix 443 444 Options Database Keys: 445 + -mat_lmvm_hist_size - the number of history vectors to keep 446 . -mat_lmvm_mult_algorithm - the algorithm to use for multiplication (recursive, dense, compact_dense) 447 . -mat_lmvm_cache_J0_products - whether products between the base Jacobian J0 and history vectors should be cached or recomputed 448 - -mat_lmvm_debug - (developer) perform internal debugging checks 449 450 Level: intermediate 451 452 Note: 453 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()` 454 paradigm instead of this routine directly. 455 456 .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBADBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`, 457 `MatCreateLMVMBFGS()`, `MatCreateLMVMBroyden()`, `MatCreateLMVMSymBroyden()` 458 @*/ 459 PetscErrorCode MatCreateLMVMBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B) 460 { 461 PetscFunctionBegin; 462 PetscCall(KSPInitializePackage()); 463 PetscCall(MatCreate(comm, B)); 464 PetscCall(MatSetSizes(*B, n, n, N, N)); 465 PetscCall(MatSetType(*B, MATLMVMBADBROYDEN)); 466 PetscCall(MatSetUp(*B)); 467 PetscFunctionReturn(PETSC_SUCCESS); 468 } 469