1 #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h> /*I "petscksp.h" I*/ 2 3 /* 4 The solution method is the matrix-free implementation of the inverse Hessian in 5 Equation 6 on page 312 of Griewank "Broyden Updating, The Good and The Bad!" 6 (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf). 7 8 Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever 9 the matrix is updated with a new (S[i], Y[i]) pair. This allows 10 repeated calls of MatSolve without incurring redundant computation. 11 12 dX <- J0^{-1} * F 13 14 for i=0,1,2,...,k 15 # Q[i] = (B_i)^{-1} * Y[i] 16 tau = (Y[i]^T F) / (Y[i]^T Y[i]) 17 dX <- dX + (tau * (S[i] - Q[i])) 18 end 19 */ 20 21 static PetscErrorCode MatSolve_LMVMBadBrdn(Mat B, Vec F, Vec dX) 22 { 23 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 24 Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx; 25 PetscInt i, j; 26 PetscScalar yjtyi, ytf; 27 28 PetscFunctionBegin; 29 VecCheckSameSize(F, 2, dX, 3); 30 VecCheckMatCompatible(B, dX, 3, F, 2); 31 32 if (lbb->needQ) { 33 /* Pre-compute (Q[i] = (B_i)^{-1} * Y[i]) */ 34 for (i = 0; i <= lmvm->k; ++i) { 35 PetscCall(MatLMVMApplyJ0Inv(B, lmvm->Y[i], lbb->Q[i])); 36 for (j = 0; j <= i - 1; ++j) { 37 PetscCall(VecDot(lmvm->Y[j], lmvm->Y[i], &yjtyi)); 38 PetscCall(VecAXPBYPCZ(lbb->Q[i], PetscRealPart(yjtyi) / lbb->yty[j], -PetscRealPart(yjtyi) / lbb->yty[j], 1.0, lmvm->S[j], lbb->Q[j])); 39 } 40 } 41 lbb->needQ = PETSC_FALSE; 42 } 43 44 PetscCall(MatLMVMApplyJ0Inv(B, F, dX)); 45 for (i = 0; i <= lmvm->k; ++i) { 46 PetscCall(VecDot(lmvm->Y[i], F, &ytf)); 47 PetscCall(VecAXPBYPCZ(dX, PetscRealPart(ytf) / lbb->yty[i], -PetscRealPart(ytf) / lbb->yty[i], 1.0, lmvm->S[i], lbb->Q[i])); 48 } 49 PetscFunctionReturn(PETSC_SUCCESS); 50 } 51 52 /* 53 The forward product is the matrix-free implementation of the direct update in 54 Equation 6 on page 302 of Griewank "Broyden Updating, The Good and The Bad!" 55 (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf). 56 57 P[i] = (B_i)*S[i] terms are computed ahead of time whenever 58 the matrix is updated with a new (S[i], Y[i]) pair. This allows 59 repeated calls of MatMult inside KSP solvers without unnecessarily 60 recomputing P[i] terms in expensive nested-loops. 61 62 Z <- J0 * X 63 64 for i=0,1,2,...,k 65 # P[i] = B_i * S[i] 66 tau = (Y[i]^T X) / (Y[i]^T S[i]) 67 dX <- dX + (tau * (Y[i] - P[i])) 68 end 69 */ 70 71 static PetscErrorCode MatMult_LMVMBadBrdn(Mat B, Vec X, Vec Z) 72 { 73 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 74 Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx; 75 PetscInt i, j; 76 PetscScalar yjtsi, ytx; 77 78 PetscFunctionBegin; 79 VecCheckSameSize(X, 2, Z, 3); 80 VecCheckMatCompatible(B, X, 2, Z, 3); 81 82 if (lbb->needP) { 83 /* Pre-compute (P[i] = (B_i) * S[i]) */ 84 for (i = 0; i <= lmvm->k; ++i) { 85 PetscCall(MatLMVMApplyJ0Fwd(B, lmvm->S[i], lbb->P[i])); 86 for (j = 0; j <= i - 1; ++j) { 87 PetscCall(VecDot(lmvm->Y[j], lmvm->S[i], &yjtsi)); 88 PetscCall(VecAXPBYPCZ(lbb->P[i], PetscRealPart(yjtsi) / lbb->yts[j], -PetscRealPart(yjtsi) / lbb->yts[j], 1.0, lmvm->Y[j], lbb->P[j])); 89 } 90 } 91 lbb->needP = PETSC_FALSE; 92 } 93 94 PetscCall(MatLMVMApplyJ0Fwd(B, X, Z)); 95 for (i = 0; i <= lmvm->k; ++i) { 96 PetscCall(VecDot(lmvm->Y[i], X, &ytx)); 97 PetscCall(VecAXPBYPCZ(Z, PetscRealPart(ytx) / lbb->yts[i], -PetscRealPart(ytx) / lbb->yts[i], 1.0, lmvm->Y[i], lbb->P[i])); 98 } 99 PetscFunctionReturn(PETSC_SUCCESS); 100 } 101 102 static PetscErrorCode MatUpdate_LMVMBadBrdn(Mat B, Vec X, Vec F) 103 { 104 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 105 Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx; 106 PetscInt old_k, i; 107 PetscScalar yty, yts; 108 109 PetscFunctionBegin; 110 if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS); 111 if (lmvm->prev_set) { 112 /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */ 113 PetscCall(VecAYPX(lmvm->Xprev, -1.0, X)); 114 PetscCall(VecAYPX(lmvm->Fprev, -1.0, F)); 115 /* Accept the update */ 116 lbb->needP = lbb->needQ = PETSC_TRUE; 117 old_k = lmvm->k; 118 PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev)); 119 /* If we hit the memory limit, shift the yty and yts arrays */ 120 if (old_k == lmvm->k) { 121 for (i = 0; i <= lmvm->k - 1; ++i) { 122 lbb->yty[i] = lbb->yty[i + 1]; 123 lbb->yts[i] = lbb->yts[i + 1]; 124 } 125 } 126 /* Accumulate the latest yTy and yTs dot products */ 127 PetscCall(VecDotBegin(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty)); 128 PetscCall(VecDotBegin(lmvm->Y[lmvm->k], lmvm->S[lmvm->k], &yts)); 129 PetscCall(VecDotEnd(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty)); 130 PetscCall(VecDotEnd(lmvm->Y[lmvm->k], lmvm->S[lmvm->k], &yts)); 131 lbb->yty[lmvm->k] = PetscRealPart(yty); 132 lbb->yts[lmvm->k] = PetscRealPart(yts); 133 } 134 /* Save the solution and function to be used in the next update */ 135 PetscCall(VecCopy(X, lmvm->Xprev)); 136 PetscCall(VecCopy(F, lmvm->Fprev)); 137 lmvm->prev_set = PETSC_TRUE; 138 PetscFunctionReturn(PETSC_SUCCESS); 139 } 140 141 static PetscErrorCode MatCopy_LMVMBadBrdn(Mat B, Mat M, MatStructure str) 142 { 143 Mat_LMVM *bdata = (Mat_LMVM *)B->data; 144 Mat_Brdn *bctx = (Mat_Brdn *)bdata->ctx; 145 Mat_LMVM *mdata = (Mat_LMVM *)M->data; 146 Mat_Brdn *mctx = (Mat_Brdn *)mdata->ctx; 147 PetscInt i; 148 149 PetscFunctionBegin; 150 mctx->needP = bctx->needP; 151 mctx->needQ = bctx->needQ; 152 for (i = 0; i <= bdata->k; ++i) { 153 mctx->yty[i] = bctx->yty[i]; 154 mctx->yts[i] = bctx->yts[i]; 155 PetscCall(VecCopy(bctx->P[i], mctx->P[i])); 156 PetscCall(VecCopy(bctx->Q[i], mctx->Q[i])); 157 } 158 PetscFunctionReturn(PETSC_SUCCESS); 159 } 160 161 static PetscErrorCode MatReset_LMVMBadBrdn(Mat B, PetscBool destructive) 162 { 163 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 164 Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx; 165 166 PetscFunctionBegin; 167 lbb->needP = lbb->needQ = PETSC_TRUE; 168 if (destructive && lbb->allocated) { 169 PetscCall(PetscFree2(lbb->yty, lbb->yts)); 170 PetscCall(VecDestroyVecs(lmvm->m, &lbb->P)); 171 PetscCall(VecDestroyVecs(lmvm->m, &lbb->Q)); 172 lbb->allocated = PETSC_FALSE; 173 } 174 PetscCall(MatReset_LMVM(B, destructive)); 175 PetscFunctionReturn(PETSC_SUCCESS); 176 } 177 178 static PetscErrorCode MatAllocate_LMVMBadBrdn(Mat B, Vec X, Vec F) 179 { 180 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 181 Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx; 182 183 PetscFunctionBegin; 184 PetscCall(MatAllocate_LMVM(B, X, F)); 185 if (!lbb->allocated) { 186 PetscCall(PetscMalloc2(lmvm->m, &lbb->yty, lmvm->m, &lbb->yts)); 187 if (lmvm->m > 0) { 188 PetscCall(VecDuplicateVecs(X, lmvm->m, &lbb->P)); 189 PetscCall(VecDuplicateVecs(X, lmvm->m, &lbb->Q)); 190 } 191 lbb->allocated = PETSC_TRUE; 192 } 193 PetscFunctionReturn(PETSC_SUCCESS); 194 } 195 196 static PetscErrorCode MatDestroy_LMVMBadBrdn(Mat B) 197 { 198 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 199 Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx; 200 201 PetscFunctionBegin; 202 if (lbb->allocated) { 203 PetscCall(PetscFree2(lbb->yty, lbb->yts)); 204 PetscCall(VecDestroyVecs(lmvm->m, &lbb->P)); 205 PetscCall(VecDestroyVecs(lmvm->m, &lbb->Q)); 206 lbb->allocated = PETSC_FALSE; 207 } 208 PetscCall(PetscFree(lmvm->ctx)); 209 PetscCall(MatDestroy_LMVM(B)); 210 PetscFunctionReturn(PETSC_SUCCESS); 211 } 212 213 static PetscErrorCode MatSetUp_LMVMBadBrdn(Mat B) 214 { 215 Mat_LMVM *lmvm = (Mat_LMVM *)B->data; 216 Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx; 217 218 PetscFunctionBegin; 219 PetscCall(MatSetUp_LMVM(B)); 220 if (!lbb->allocated) { 221 PetscCall(PetscMalloc2(lmvm->m, &lbb->yty, lmvm->m, &lbb->yts)); 222 if (lmvm->m > 0) { 223 PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbb->P)); 224 PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbb->Q)); 225 } 226 lbb->allocated = PETSC_TRUE; 227 } 228 PetscFunctionReturn(PETSC_SUCCESS); 229 } 230 231 PetscErrorCode MatCreate_LMVMBadBrdn(Mat B) 232 { 233 Mat_LMVM *lmvm; 234 Mat_Brdn *lbb; 235 236 PetscFunctionBegin; 237 PetscCall(MatCreate_LMVM(B)); 238 PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMBADBROYDEN)); 239 B->ops->setup = MatSetUp_LMVMBadBrdn; 240 B->ops->destroy = MatDestroy_LMVMBadBrdn; 241 B->ops->solve = MatSolve_LMVMBadBrdn; 242 243 lmvm = (Mat_LMVM *)B->data; 244 lmvm->square = PETSC_TRUE; 245 lmvm->ops->allocate = MatAllocate_LMVMBadBrdn; 246 lmvm->ops->reset = MatReset_LMVMBadBrdn; 247 lmvm->ops->mult = MatMult_LMVMBadBrdn; 248 lmvm->ops->update = MatUpdate_LMVMBadBrdn; 249 lmvm->ops->copy = MatCopy_LMVMBadBrdn; 250 251 PetscCall(PetscNew(&lbb)); 252 lmvm->ctx = (void *)lbb; 253 lbb->allocated = PETSC_FALSE; 254 lbb->needP = lbb->needQ = PETSC_TRUE; 255 PetscFunctionReturn(PETSC_SUCCESS); 256 } 257 258 /*@ 259 MatCreateLMVMBadBroyden - Creates a limited-memory modified (aka "bad") Broyden-type 260 approximation matrix used for a Jacobian. L-BadBrdn is not guaranteed to be 261 symmetric or positive-definite. 262 263 To use the L-BadBrdn matrix with other vector types, the matrix must be 264 created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`. 265 This ensures that the internal storage and work vectors are duplicated from the 266 correct type of vector. 267 268 Collective 269 270 Input Parameters: 271 + comm - MPI communicator 272 . n - number of local rows for storage vectors 273 - N - global size of the storage vectors 274 275 Output Parameter: 276 . B - the matrix 277 278 Options Database Keys: 279 + -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal) 280 . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling 281 . -mat_lmvm_rho - (developer) update limiter for the J0 scaling 282 . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling 283 . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling 284 - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling 285 286 Level: intermediate 287 288 Note: 289 It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()` 290 paradigm instead of this routine directly. 291 292 .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBADBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`, 293 `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMSymBrdn()` 294 @*/ 295 PetscErrorCode MatCreateLMVMBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B) 296 { 297 PetscFunctionBegin; 298 PetscCall(MatCreate(comm, B)); 299 PetscCall(MatSetSizes(*B, n, n, N, N)); 300 PetscCall(MatSetType(*B, MATLMVMBADBROYDEN)); 301 PetscCall(MatSetUp(*B)); 302 PetscFunctionReturn(PETSC_SUCCESS); 303 } 304