Actual source code: dfp.c
1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
4: /*
5: Limited-memory Davidon-Fletcher-Powell method for approximating both
6: the forward product and inverse application of a Jacobian.
7: */
9: /*
10: The solution method (approximate inverse Jacobian application) is
11: matrix-vector product version of the recursive formula given in
12: Equation (6.15) of Nocedal and Wright "Numerical Optimization" 2nd
13: edition, pg 139.
15: Note: Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
16: the matrix is updated with a new (S[i], Y[i]) pair. This allows
17: repeated calls of MatSolve without incurring redundant computation.
19: dX <- J0^{-1} * F
21: for i = 0,1,2,...,k
22: # Q[i] = (B_i)^{-1} * Y[i]
23: gamma = (S[i]^T F) / (Y[i]^T S[i])
24: zeta = (Y[i]^T dX) / (Y[i]^T Q[i])
25: dX <- dX + (gamma * S[i]) - (zeta * Q[i])
26: end
27: */
28: PetscErrorCode MatSolve_LMVMDFP(Mat B, Vec F, Vec dX)
29: {
30: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
31: Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
32: PetscInt i, j;
33: PetscScalar yjtqi, sjtyi, ytx, stf, ytq;
35: PetscFunctionBegin;
38: VecCheckSameSize(F, 2, dX, 3);
39: VecCheckMatCompatible(B, dX, 3, F, 2);
41: if (ldfp->needQ) {
42: /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
43: for (i = 0; i <= lmvm->k; ++i) {
44: PetscCall(MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], ldfp->Q[i]));
45: /* Compute the necessary dot products */
46: PetscCall(VecMDot(lmvm->Y[i], i, lmvm->S, ldfp->workscalar));
47: for (j = 0; j < i; ++j) {
48: sjtyi = ldfp->workscalar[j];
49: PetscCall(VecDot(lmvm->Y[j], ldfp->Q[i], &yjtqi));
50: /* Compute the pure DFP component of the inverse application*/
51: PetscCall(VecAXPBYPCZ(ldfp->Q[i], -PetscRealPart(yjtqi) / ldfp->ytq[j], PetscRealPart(sjtyi) / ldfp->yts[j], 1.0, ldfp->Q[j], lmvm->S[j]));
52: }
53: PetscCall(VecDot(lmvm->Y[i], ldfp->Q[i], &ytq));
54: ldfp->ytq[i] = PetscRealPart(ytq);
55: }
56: ldfp->needQ = PETSC_FALSE;
57: }
59: /* Start the outer loop (i) for the recursive formula */
60: PetscCall(MatSymBrdnApplyJ0Inv(B, F, dX));
61: /* Get all the dot products we need */
62: PetscCall(VecMDot(F, lmvm->k + 1, lmvm->S, ldfp->workscalar));
63: for (i = 0; i <= lmvm->k; ++i) {
64: stf = ldfp->workscalar[i];
65: PetscCall(VecDot(lmvm->Y[i], dX, &ytx));
66: /* Update dX_{i+1} = (B^{-1})_{i+1} * f */
67: PetscCall(VecAXPBYPCZ(dX, -PetscRealPart(ytx) / ldfp->ytq[i], PetscRealPart(stf) / ldfp->yts[i], 1.0, ldfp->Q[i], lmvm->S[i]));
68: }
69: PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: /*
73: The forward product for the approximate Jacobian is the matrix-free
74: implementation of the recursive formula given in Equation 6.13 of
75: Nocedal and Wright "Numerical Optimization" 2nd edition, pg 139.
77: This forward product has a two-loop form similar to the BFGS two-loop
78: formulation for the inverse Jacobian application. However, the S and
79: Y vectors have interchanged roles.
81: work <- X
83: for i = k,k-1,k-2,...,0
84: rho[i] = 1 / (Y[i]^T S[i])
85: alpha[i] = rho[i] * (Y[i]^T work)
86: work <- work - (alpha[i] * S[i])
87: end
89: Z <- J0 * work
91: for i = 0,1,2,...,k
92: beta = rho[i] * (S[i]^T Y)
93: Z <- Z + ((alpha[i] - beta) * Y[i])
94: end
95: */
96: PetscErrorCode MatMult_LMVMDFP(Mat B, Vec X, Vec Z)
97: {
98: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
99: Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
100: PetscInt i;
101: PetscReal *alpha, beta;
102: PetscScalar ytx, stz;
104: PetscFunctionBegin;
105: /* Copy the function into the work vector for the first loop */
106: PetscCall(VecCopy(X, ldfp->work));
108: /* Start the first loop */
109: PetscCall(PetscMalloc1(lmvm->k + 1, &alpha));
110: for (i = lmvm->k; i >= 0; --i) {
111: PetscCall(VecDot(lmvm->Y[i], ldfp->work, &ytx));
112: alpha[i] = PetscRealPart(ytx) / ldfp->yts[i];
113: PetscCall(VecAXPY(ldfp->work, -alpha[i], lmvm->S[i]));
114: }
116: /* Apply the forward product with initial Jacobian */
117: PetscCall(MatSymBrdnApplyJ0Fwd(B, ldfp->work, Z));
119: /* Start the second loop */
120: for (i = 0; i <= lmvm->k; ++i) {
121: PetscCall(VecDot(lmvm->S[i], Z, &stz));
122: beta = PetscRealPart(stz) / ldfp->yts[i];
123: PetscCall(VecAXPY(Z, alpha[i] - beta, lmvm->Y[i]));
124: }
125: PetscCall(PetscFree(alpha));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: static PetscErrorCode MatUpdate_LMVMDFP(Mat B, Vec X, Vec F)
130: {
131: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
132: Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
133: Mat_LMVM *dbase;
134: Mat_DiagBrdn *dctx;
135: PetscInt old_k, i;
136: PetscReal curvtol, ytytmp;
137: PetscScalar curvature, ststmp;
139: PetscFunctionBegin;
140: if (!lmvm->m) PetscFunctionReturn(PETSC_SUCCESS);
141: if (lmvm->prev_set) {
142: /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
143: PetscCall(VecAYPX(lmvm->Xprev, -1.0, X));
144: PetscCall(VecAYPX(lmvm->Fprev, -1.0, F));
146: /* Test if the updates can be accepted */
147: PetscCall(VecDotNorm2(lmvm->Xprev, lmvm->Fprev, &curvature, &ytytmp));
148: if (ytytmp < lmvm->eps) curvtol = 0.0;
149: else curvtol = lmvm->eps * ytytmp;
151: if (PetscRealPart(curvature) > curvtol) {
152: /* Update is good, accept it */
153: ldfp->watchdog = 0;
154: ldfp->needQ = PETSC_TRUE;
155: old_k = lmvm->k;
156: PetscCall(MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev));
157: /* If we hit the memory limit, shift the yts, yty and sts arrays */
158: if (old_k == lmvm->k) {
159: for (i = 0; i <= lmvm->k - 1; ++i) {
160: ldfp->yts[i] = ldfp->yts[i + 1];
161: ldfp->yty[i] = ldfp->yty[i + 1];
162: ldfp->sts[i] = ldfp->sts[i + 1];
163: }
164: }
165: /* Update history of useful scalars */
166: PetscCall(VecDot(lmvm->S[lmvm->k], lmvm->S[lmvm->k], &ststmp));
167: ldfp->yts[lmvm->k] = PetscRealPart(curvature);
168: ldfp->yty[lmvm->k] = ytytmp;
169: ldfp->sts[lmvm->k] = PetscRealPart(ststmp);
170: /* Compute the scalar scale if necessary */
171: if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_SCALAR) PetscCall(MatSymBrdnComputeJ0Scalar(B));
172: } else {
173: /* Update is bad, skip it */
174: ++lmvm->nrejects;
175: ++ldfp->watchdog;
176: }
177: } else {
178: switch (ldfp->scale_type) {
179: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
180: dbase = (Mat_LMVM *)ldfp->D->data;
181: dctx = (Mat_DiagBrdn *)dbase->ctx;
182: PetscCall(VecSet(dctx->invD, ldfp->delta));
183: break;
184: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
185: ldfp->sigma = ldfp->delta;
186: break;
187: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
188: ldfp->sigma = 1.0;
189: break;
190: default:
191: break;
192: }
193: }
195: /* Update the scaling */
196: if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMUpdate(ldfp->D, X, F));
198: if (ldfp->watchdog > ldfp->max_seq_rejects) {
199: PetscCall(MatLMVMReset(B, PETSC_FALSE));
200: if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) PetscCall(MatLMVMReset(ldfp->D, PETSC_FALSE));
201: }
203: /* Save the solution and function to be used in the next update */
204: PetscCall(VecCopy(X, lmvm->Xprev));
205: PetscCall(VecCopy(F, lmvm->Fprev));
206: lmvm->prev_set = PETSC_TRUE;
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: static PetscErrorCode MatCopy_LMVMDFP(Mat B, Mat M, MatStructure str)
211: {
212: Mat_LMVM *bdata = (Mat_LMVM *)B->data;
213: Mat_SymBrdn *bctx = (Mat_SymBrdn *)bdata->ctx;
214: Mat_LMVM *mdata = (Mat_LMVM *)M->data;
215: Mat_SymBrdn *mctx = (Mat_SymBrdn *)mdata->ctx;
216: PetscInt i;
218: PetscFunctionBegin;
219: mctx->needQ = bctx->needQ;
220: for (i = 0; i <= bdata->k; ++i) {
221: mctx->ytq[i] = bctx->ytq[i];
222: mctx->yts[i] = bctx->yts[i];
223: PetscCall(VecCopy(bctx->Q[i], mctx->Q[i]));
224: }
225: mctx->scale_type = bctx->scale_type;
226: mctx->alpha = bctx->alpha;
227: mctx->beta = bctx->beta;
228: mctx->rho = bctx->rho;
229: mctx->sigma_hist = bctx->sigma_hist;
230: mctx->watchdog = bctx->watchdog;
231: mctx->max_seq_rejects = bctx->max_seq_rejects;
232: switch (bctx->scale_type) {
233: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
234: mctx->sigma = bctx->sigma;
235: break;
236: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
237: PetscCall(MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN));
238: break;
239: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
240: mctx->sigma = 1.0;
241: break;
242: default:
243: break;
244: }
245: PetscFunctionReturn(PETSC_SUCCESS);
246: }
248: static PetscErrorCode MatReset_LMVMDFP(Mat B, PetscBool destructive)
249: {
250: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
251: Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
252: Mat_LMVM *dbase;
253: Mat_DiagBrdn *dctx;
255: PetscFunctionBegin;
256: ldfp->watchdog = 0;
257: ldfp->needQ = PETSC_TRUE;
258: if (ldfp->allocated) {
259: if (destructive) {
260: PetscCall(VecDestroy(&ldfp->work));
261: PetscCall(PetscFree5(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts, ldfp->workscalar));
262: PetscCall(VecDestroyVecs(lmvm->m, &ldfp->Q));
263: switch (ldfp->scale_type) {
264: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
265: PetscCall(MatLMVMReset(ldfp->D, PETSC_TRUE));
266: break;
267: default:
268: break;
269: }
270: ldfp->allocated = PETSC_FALSE;
271: } else {
272: switch (ldfp->scale_type) {
273: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
274: ldfp->sigma = ldfp->delta;
275: break;
276: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
277: PetscCall(MatLMVMReset(ldfp->D, PETSC_FALSE));
278: dbase = (Mat_LMVM *)ldfp->D->data;
279: dctx = (Mat_DiagBrdn *)dbase->ctx;
280: PetscCall(VecSet(dctx->invD, ldfp->delta));
281: break;
282: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
283: ldfp->sigma = 1.0;
284: break;
285: default:
286: break;
287: }
288: }
289: }
290: PetscCall(MatReset_LMVM(B, destructive));
291: PetscFunctionReturn(PETSC_SUCCESS);
292: }
294: static PetscErrorCode MatAllocate_LMVMDFP(Mat B, Vec X, Vec F)
295: {
296: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
297: Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
299: PetscFunctionBegin;
300: PetscCall(MatAllocate_LMVM(B, X, F));
301: if (!ldfp->allocated) {
302: PetscCall(VecDuplicate(X, &ldfp->work));
303: PetscCall(PetscMalloc5(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts, lmvm->m, &ldfp->workscalar));
304: if (lmvm->m > 0) PetscCall(VecDuplicateVecs(X, lmvm->m, &ldfp->Q));
305: switch (ldfp->scale_type) {
306: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
307: PetscCall(MatLMVMAllocate(ldfp->D, X, F));
308: break;
309: default:
310: break;
311: }
312: ldfp->allocated = PETSC_TRUE;
313: }
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: static PetscErrorCode MatDestroy_LMVMDFP(Mat B)
318: {
319: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
320: Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
322: PetscFunctionBegin;
323: if (ldfp->allocated) {
324: PetscCall(VecDestroy(&ldfp->work));
325: PetscCall(PetscFree5(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts, ldfp->workscalar));
326: PetscCall(VecDestroyVecs(lmvm->m, &ldfp->Q));
327: ldfp->allocated = PETSC_FALSE;
328: }
329: PetscCall(MatDestroy(&ldfp->D));
330: PetscCall(PetscFree(lmvm->ctx));
331: PetscCall(MatDestroy_LMVM(B));
332: PetscFunctionReturn(PETSC_SUCCESS);
333: }
335: static PetscErrorCode MatSetUp_LMVMDFP(Mat B)
336: {
337: Mat_LMVM *lmvm = (Mat_LMVM *)B->data;
338: Mat_SymBrdn *ldfp = (Mat_SymBrdn *)lmvm->ctx;
339: PetscInt n, N;
341: PetscFunctionBegin;
342: PetscCall(MatSetUp_LMVM(B));
343: if (!ldfp->allocated) {
344: PetscCall(VecDuplicate(lmvm->Xprev, &ldfp->work));
345: PetscCall(PetscMalloc5(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts, lmvm->m, &ldfp->workscalar));
346: if (lmvm->m > 0) PetscCall(VecDuplicateVecs(lmvm->Xprev, lmvm->m, &ldfp->Q));
347: switch (ldfp->scale_type) {
348: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
349: PetscCall(MatGetLocalSize(B, &n, &n));
350: PetscCall(MatGetSize(B, &N, &N));
351: PetscCall(MatSetSizes(ldfp->D, n, n, N, N));
352: PetscCall(MatSetUp(ldfp->D));
353: break;
354: default:
355: break;
356: }
357: ldfp->allocated = PETSC_TRUE;
358: }
359: PetscFunctionReturn(PETSC_SUCCESS);
360: }
362: static PetscErrorCode MatSetFromOptions_LMVMDFP(Mat B, PetscOptionItems *PetscOptionsObject)
363: {
364: PetscFunctionBegin;
365: PetscCall(MatSetFromOptions_LMVM(B, PetscOptionsObject));
366: PetscOptionsHeadBegin(PetscOptionsObject, "DFP method for approximating SPD Jacobian actions (MATLMVMDFP)");
367: PetscCall(MatSetFromOptions_LMVMSymBrdn_Private(B, PetscOptionsObject));
368: PetscOptionsHeadEnd();
369: PetscFunctionReturn(PETSC_SUCCESS);
370: }
372: PetscErrorCode MatCreate_LMVMDFP(Mat B)
373: {
374: Mat_LMVM *lmvm;
375: Mat_SymBrdn *ldfp;
377: PetscFunctionBegin;
378: PetscCall(MatCreate_LMVMSymBrdn(B));
379: PetscCall(PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP));
380: B->ops->setup = MatSetUp_LMVMDFP;
381: B->ops->destroy = MatDestroy_LMVMDFP;
382: B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
384: lmvm = (Mat_LMVM *)B->data;
385: lmvm->ops->allocate = MatAllocate_LMVMDFP;
386: lmvm->ops->reset = MatReset_LMVMDFP;
387: lmvm->ops->update = MatUpdate_LMVMDFP;
388: lmvm->ops->mult = MatMult_LMVMDFP;
389: lmvm->ops->solve = MatSolve_LMVMDFP;
390: lmvm->ops->copy = MatCopy_LMVMDFP;
392: ldfp = (Mat_SymBrdn *)lmvm->ctx;
393: ldfp->needP = PETSC_FALSE;
394: ldfp->phi = 1.0;
395: PetscFunctionReturn(PETSC_SUCCESS);
396: }
398: /*@
399: MatCreateLMVMDFP - Creates a limited-memory Davidon-Fletcher-Powell (DFP) matrix
400: used for approximating Jacobians. L-DFP is symmetric positive-definite by
401: construction, and is the dual of L-BFGS where Y and S vectors swap roles.
403: To use the L-DFP matrix with other vector types, the matrix must be
404: created using `MatCreate()` and `MatSetType()`, followed by `MatLMVMAllocate()`.
405: This ensures that the internal storage and work vectors are duplicated from the
406: correct type of vector.
408: Collective
410: Input Parameters:
411: + comm - MPI communicator
412: . n - number of local rows for storage vectors
413: - N - global size of the storage vectors
415: Output Parameter:
416: . B - the matrix
418: Options Database Keys:
419: + -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
420: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
421: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
422: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
423: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
424: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
426: Level: intermediate
428: Note:
429: It is recommended that one use the `MatCreate()`, `MatSetType()` and/or `MatSetFromOptions()`
430: paradigm instead of this routine directly.
432: .seealso: [](ch_ksp), `MatCreate()`, `MATLMVM`, `MATLMVMDFP`, `MatCreateLMVMBFGS()`, `MatCreateLMVMSR1()`,
433: `MatCreateLMVMBrdn()`, `MatCreateLMVMBadBrdn()`, `MatCreateLMVMSymBrdn()`
434: @*/
435: PetscErrorCode MatCreateLMVMDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
436: {
437: PetscFunctionBegin;
438: PetscCall(KSPInitializePackage());
439: PetscCall(MatCreate(comm, B));
440: PetscCall(MatSetSizes(*B, n, n, N, N));
441: PetscCall(MatSetType(*B, MATLMVMDFP));
442: PetscCall(MatSetUp(*B));
443: PetscFunctionReturn(PETSC_SUCCESS);
444: }