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