Actual source code: brdn.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 *sts, *stq;
12: } Mat_Brdn;
14: /*------------------------------------------------------------*/
16: /*
17: The solution method is the matrix-free implementation of the inverse Hessian
18: representation in 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 = (S[i]^T dX) / (S[i]^T Q[i])
30: dX <- dX + (tau * (S[i] - Q[i]))
31: end
32: */
34: static PetscErrorCode MatSolve_LMVMBrdn(Mat B, Vec F, Vec dX)
35: {
36: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
37: Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
38: PetscErrorCode ierr;
39: PetscInt i, j;
40: PetscScalar sjtqi, stx, stq;
41:
43: VecCheckSameSize(F, 2, dX, 3);
44: VecCheckMatCompatible(B, dX, 3, F, 2);
45:
46: if (lbrdn->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], lbrdn->Q[i]);
50: for (j = 0; j <= i-1; ++j) {
51: VecDot(lmvm->S[j], lbrdn->Q[i], &sjtqi);
52: VecAXPBYPCZ(lbrdn->Q[i], PetscRealPart(sjtqi)/lbrdn->stq[j], -PetscRealPart(sjtqi)/lbrdn->stq[j], 1.0, lmvm->S[j], lbrdn->Q[j]);
53: }
54: VecDot(lmvm->S[i], lbrdn->Q[i], &stq);
55: lbrdn->stq[i] = PetscRealPart(stq);
56: }
57: lbrdn->needQ = PETSC_FALSE;
58: }
59:
60: MatLMVMApplyJ0Inv(B, F, dX);
61: for (i = 0; i <= lmvm->k-1; ++i) {
62: VecDot(lmvm->S[i], dX, &stx);
63: VecAXPBYPCZ(dX, PetscRealPart(stx)/lbrdn->stq[i], -PetscRealPart(stx)/lbrdn->stq[i], 1.0, lmvm->S[i], lbrdn->Q[i]);
64: }
65: return(0);
66: }
68: /*------------------------------------------------------------*/
70: /*
71: The forward product is the matrix-free implementation of Equation 2 in
72: page 302 of Griewank "Broyden Updating, The Good and The Bad!"
73: (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).
74:
75: P[i] = (B_i)*S[i] terms are computed ahead of time whenever
76: the matrix is updated with a new (S[i], Y[i]) pair. This allows
77: repeated calls of MatMult inside KSP solvers without unnecessarily
78: recomputing P[i] terms in expensive nested-loops.
79:
80: Z <- J0 * X
81:
82: for i=0,1,2,...,k
83: # P[i] = B_i * S[i]
84: tau = (S[i]^T X) / (S[i]^T S[i])
85: dX <- dX + (tau * (Y[i] - P[i]))
86: end
87: */
89: static PetscErrorCode MatMult_LMVMBrdn(Mat B, Vec X, Vec Z)
90: {
91: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
92: Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
93: PetscErrorCode ierr;
94: PetscInt i, j;
95: PetscScalar sjtsi, stx;
96:
98: VecCheckSameSize(X, 2, Z, 3);
99: VecCheckMatCompatible(B, X, 2, Z, 3);
100:
101: if (lbrdn->needP) {
102: /* Pre-compute (P[i] = (B_i) * S[i]) */
103: for (i = 0; i <= lmvm->k; ++i) {
104: MatLMVMApplyJ0Fwd(B, lmvm->S[i], lbrdn->P[i]);
105: for (j = 0; j <= i-1; ++j) {
106: VecDot(lmvm->S[j], lmvm->S[i], &sjtsi);
107: VecAXPBYPCZ(lbrdn->P[i], PetscRealPart(sjtsi)/lbrdn->sts[j], -PetscRealPart(sjtsi)/lbrdn->sts[j], 1.0, lmvm->Y[j], lbrdn->P[j]);
108: }
109: }
110: lbrdn->needP = PETSC_FALSE;
111: }
112:
113: MatLMVMApplyJ0Fwd(B, X, Z);
114: for (i = 0; i <= lmvm->k-1; ++i) {
115: VecDot(lmvm->S[i], X, &stx);
116: VecAXPBYPCZ(Z, PetscRealPart(stx)/lbrdn->sts[i], -PetscRealPart(stx)/lbrdn->sts[i], 1.0, lmvm->Y[i], lbrdn->P[i]);
117: }
118: return(0);
119: }
121: /*------------------------------------------------------------*/
123: static PetscErrorCode MatUpdate_LMVMBrdn(Mat B, Vec X, Vec F)
124: {
125: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
126: Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
127: PetscErrorCode ierr;
128: PetscInt old_k, i;
129: PetscScalar sts;
132: if (!lmvm->m) return(0);
133: if (lmvm->prev_set) {
134: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
135: VecAYPX(lmvm->Xprev, -1.0, X);
136: VecAYPX(lmvm->Fprev, -1.0, F);
137: /* Accept the update */
138: lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
139: old_k = lmvm->k;
140: MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
141: /* If we hit the memory limit, shift the sts array */
142: if (old_k == lmvm->k) {
143: for (i = 0; i <= lmvm->k-1; ++i) {
144: lbrdn->sts[i] = lbrdn->sts[i+1];
145: }
146: }
147: VecDot(lmvm->S[lmvm->k], lmvm->S[lmvm->k], &sts);
148: lbrdn->sts[lmvm->k] = PetscRealPart(sts);
149: }
150: /* Save the solution and function to be used in the next update */
151: VecCopy(X, lmvm->Xprev);
152: VecCopy(F, lmvm->Fprev);
153: lmvm->prev_set = PETSC_TRUE;
154: return(0);
155: }
157: /*------------------------------------------------------------*/
159: static PetscErrorCode MatCopy_LMVMBrdn(Mat B, Mat M, MatStructure str)
160: {
161: Mat_LMVM *bdata = (Mat_LMVM*)B->data;
162: Mat_Brdn *bctx = (Mat_Brdn*)bdata->ctx;
163: Mat_LMVM *mdata = (Mat_LMVM*)M->data;
164: Mat_Brdn *mctx = (Mat_Brdn*)mdata->ctx;
165: PetscErrorCode ierr;
166: PetscInt i;
169: mctx->needP = bctx->needP;
170: mctx->needQ = bctx->needQ;
171: for (i=0; i<=bdata->k; ++i) {
172: mctx->sts[i] = bctx->sts[i];
173: mctx->stq[i] = bctx->stq[i];
174: VecCopy(bctx->P[i], mctx->P[i]);
175: VecCopy(bctx->Q[i], mctx->Q[i]);
176: }
177: return(0);
178: }
180: /*------------------------------------------------------------*/
182: static PetscErrorCode MatReset_LMVMBrdn(Mat B, PetscBool destructive)
183: {
184: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
185: Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
186: PetscErrorCode ierr;
187:
189: lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
190: if (destructive && lbrdn->allocated) {
191: PetscFree2(lbrdn->sts, lbrdn->stq);
192: VecDestroyVecs(lmvm->m, &lbrdn->P);
193: VecDestroyVecs(lmvm->m, &lbrdn->Q);
194: lbrdn->allocated = PETSC_FALSE;
195: }
196: MatReset_LMVM(B, destructive);
197: return(0);
198: }
200: /*------------------------------------------------------------*/
202: static PetscErrorCode MatAllocate_LMVMBrdn(Mat B, Vec X, Vec F)
203: {
204: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
205: Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
206: PetscErrorCode ierr;
207:
209: MatAllocate_LMVM(B, X, F);
210: if (!lbrdn->allocated) {
211: PetscMalloc2(lmvm->m, &lbrdn->sts, lmvm->m, &lbrdn->stq);
212: if (lmvm->m > 0) {
213: VecDuplicateVecs(X, lmvm->m, &lbrdn->P);
214: VecDuplicateVecs(X, lmvm->m, &lbrdn->Q);
215: }
216: lbrdn->allocated = PETSC_TRUE;
217: }
218: return(0);
219: }
221: /*------------------------------------------------------------*/
223: static PetscErrorCode MatDestroy_LMVMBrdn(Mat B)
224: {
225: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
226: Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
227: PetscErrorCode ierr;
230: if (lbrdn->allocated) {
231: PetscFree2(lbrdn->sts, lbrdn->stq);
232: VecDestroyVecs(lmvm->m, &lbrdn->P);
233: VecDestroyVecs(lmvm->m, &lbrdn->Q);
234: lbrdn->allocated = PETSC_FALSE;
235: }
236: PetscFree(lmvm->ctx);
237: MatDestroy_LMVM(B);
238: return(0);
239: }
241: /*------------------------------------------------------------*/
243: static PetscErrorCode MatSetUp_LMVMBrdn(Mat B)
244: {
245: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
246: Mat_Brdn *lbrdn = (Mat_Brdn*)lmvm->ctx;
247: PetscErrorCode ierr;
248:
250: MatSetUp_LMVM(B);
251: if (!lbrdn->allocated) {
252: PetscMalloc2(lmvm->m, &lbrdn->sts, lmvm->m, &lbrdn->stq);
253: if (lmvm->m > 0) {
254: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbrdn->P);
255: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbrdn->Q);
256: }
257: lbrdn->allocated = PETSC_TRUE;
258: }
259: return(0);
260: }
262: /*------------------------------------------------------------*/
264: PetscErrorCode MatCreate_LMVMBrdn(Mat B)
265: {
266: Mat_LMVM *lmvm;
267: Mat_Brdn *lbrdn;
268: PetscErrorCode ierr;
271: MatCreate_LMVM(B);
272: PetscObjectChangeTypeName((PetscObject)B, MATLMVMBRDN);
273: B->ops->setup = MatSetUp_LMVMBrdn;
274: B->ops->destroy = MatDestroy_LMVMBrdn;
275: B->ops->solve = MatSolve_LMVMBrdn;
277: lmvm = (Mat_LMVM*)B->data;
278: lmvm->square = PETSC_TRUE;
279: lmvm->ops->allocate = MatAllocate_LMVMBrdn;
280: lmvm->ops->reset = MatReset_LMVMBrdn;
281: lmvm->ops->mult = MatMult_LMVMBrdn;
282: lmvm->ops->update = MatUpdate_LMVMBrdn;
283: lmvm->ops->copy = MatCopy_LMVMBrdn;
285: PetscNewLog(B, &lbrdn);
286: lmvm->ctx = (void*)lbrdn;
287: lbrdn->allocated = PETSC_FALSE;
288: lbrdn->needP = lbrdn->needQ = PETSC_TRUE;
289: return(0);
290: }
292: /*------------------------------------------------------------*/
294: /*@
295: MatCreateLMVMBrdn - Creates a limited-memory "good" Broyden-type approximation
296: matrix used for a Jacobian. L-Brdn is not guaranteed to be symmetric or
297: positive-definite.
298:
299: The provided local and global sizes must match the solution and function vectors
300: used with MatLMVMUpdate() and MatSolve(). The resulting L-Brdn matrix will have
301: storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
302: parallel. To use the L-Brdn matrix with other vector types, the matrix must be
303: created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
304: This ensures that the internal storage and work vectors are duplicated from the
305: correct type of vector.
307: Collective on MPI_Comm
309: Input Parameters:
310: + comm - MPI communicator, set to PETSC_COMM_SELF
311: . n - number of local rows for storage vectors
312: - N - global size of the storage vectors
314: Output Parameter:
315: . B - the matrix
317: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
318: paradigm instead of this routine directly.
320: Options Database Keys:
321: . -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
323: Level: intermediate
325: .seealso: MatCreate(), MATLMVM, MATLMVMBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
326: MatCreateLMVMBFGS(), MatCreateLMVMBadBrdn(), MatCreateLMVMSymBrdn()
327: @*/
328: PetscErrorCode MatCreateLMVMBrdn(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
329: {
330: PetscErrorCode ierr;
331:
333: MatCreate(comm, B);
334: MatSetSizes(*B, n, n, N, N);
335: MatSetType(*B, MATLMVMBRDN);
336: MatSetUp(*B);
337: return(0);
338: }