Actual source code: badbrdn.c
1: #include <../src/ksp/ksp/utils/lmvm/brdn/brdn.h>
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).
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.
12: dX <- J0^{-1} * F
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: */
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;
28: PetscFunctionBegin;
29: VecCheckSameSize(F, 2, dX, 3);
30: VecCheckMatCompatible(B, dX, 3, F, 2);
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: }
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: }
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).
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.
62: Z <- J0 * X
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: */
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;
78: PetscFunctionBegin;
79: VecCheckSameSize(X, 2, Z, 3);
80: VecCheckMatCompatible(B, X, 2, Z, 3);
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: }
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: }
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;
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: }
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;
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: }
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;
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: }
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;
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: }
196: static PetscErrorCode MatDestroy_LMVMBadBrdn(Mat B)
197: {
198: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
199: Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx;
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: }
213: static PetscErrorCode MatSetUp_LMVMBadBrdn(Mat B)
214: {
215: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
216: Mat_Brdn *lbb = (Mat_Brdn *)lmvm->ctx;
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: }
231: PetscErrorCode MatCreate_LMVMBadBrdn(Mat B)
232: {
233: Mat_LMVM *lmvm;
234: Mat_Brdn *lbb;
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;
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->solve = MatSolve_LMVMBadBrdn;
249: lmvm->ops->update = MatUpdate_LMVMBadBrdn;
250: lmvm->ops->copy = MatCopy_LMVMBadBrdn;
252: PetscCall(PetscNew(&lbb));
253: lmvm->ctx = (void *)lbb;
254: lbb->allocated = PETSC_FALSE;
255: lbb->needP = lbb->needQ = PETSC_TRUE;
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }
259: /*@
260: MatCreateLMVMBadBroyden - Creates a limited-memory modified (aka "bad") Broyden-type
261: approximation matrix used for a Jacobian. L-BadBrdn is not guaranteed to be
262: symmetric or positive-definite.
264: To use the L-BadBrdn matrix with other vector types, the matrix must be
265: created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
266: This ensures that the internal storage and work vectors are duplicated from the
267: correct type of vector.
269: Collective
271: Input Parameters:
272: + comm - MPI communicator
273: . n - number of local rows for storage vectors
274: - N - global size of the storage vectors
276: Output Parameter:
277: . B - the matrix
279: Options Database Keys:
280: + -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
281: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
282: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
283: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
284: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
285: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
287: Level: intermediate
289: Note:
290: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
291: paradigm instead of this routine directly.
293: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMBADBRDN`, `MatCreateLMVMDFP()`, `MatCreateLMVMSR1()`,
294: `MatCreateLMVMBFGS()`, `MatCreateLMVMBrdn()`, `MatCreateLMVMSymBrdn()`
295: @*/
296: PetscErrorCode MatCreateLMVMBadBroyden(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
297: {
298: PetscFunctionBegin;
299: PetscCall(KSPInitializePackage());
300: PetscCall(MatCreate(comm, B));
301: PetscCall(MatSetSizes(*B, n, n, N, N));
302: PetscCall(MatSetType(*B, MATLMVMBADBROYDEN));
303: PetscCall(MatSetUp(*B));
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }