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