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