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