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