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) {
181: MatSymBrdnComputeJ0Scalar(B);
182: }
183: } else {
184: /* Update is bad, skip it */
185: ++lmvm->nrejects;
186: ++ldfp->watchdog;
187: }
188: } else {
189: switch (ldfp->scale_type) {
190: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
191: dbase = (Mat_LMVM*)ldfp->D->data;
192: dctx = (Mat_DiagBrdn*)dbase->ctx;
193: VecSet(dctx->invD, ldfp->delta);
194: break;
195: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
196: ldfp->sigma = ldfp->delta;
197: break;
198: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
199: ldfp->sigma = 1.0;
200: break;
201: default:
202: break;
203: }
204: }
206: /* Update the scaling */
207: if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
208: MatLMVMUpdate(ldfp->D, X, F);
209: }
211: if (ldfp->watchdog > ldfp->max_seq_rejects) {
212: MatLMVMReset(B, PETSC_FALSE);
213: if (ldfp->scale_type == MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL) {
214: MatLMVMReset(ldfp->D, PETSC_FALSE);
215: }
216: }
218: /* Save the solution and function to be used in the next update */
219: VecCopy(X, lmvm->Xprev);
220: VecCopy(F, lmvm->Fprev);
221: lmvm->prev_set = PETSC_TRUE;
222: return 0;
223: }
225: /*------------------------------------------------------------*/
227: static PetscErrorCode MatCopy_LMVMDFP(Mat B, Mat M, MatStructure str)
228: {
229: Mat_LMVM *bdata = (Mat_LMVM*)B->data;
230: Mat_SymBrdn *bctx = (Mat_SymBrdn*)bdata->ctx;
231: Mat_LMVM *mdata = (Mat_LMVM*)M->data;
232: Mat_SymBrdn *mctx = (Mat_SymBrdn*)mdata->ctx;
233: PetscInt i;
235: mctx->needQ = bctx->needQ;
236: for (i=0; i<=bdata->k; ++i) {
237: mctx->ytq[i] = bctx->ytq[i];
238: mctx->yts[i] = bctx->yts[i];
239: VecCopy(bctx->Q[i], mctx->Q[i]);
240: }
241: mctx->scale_type = bctx->scale_type;
242: mctx->alpha = bctx->alpha;
243: mctx->beta = bctx->beta;
244: mctx->rho = bctx->rho;
245: mctx->sigma_hist = bctx->sigma_hist;
246: mctx->watchdog = bctx->watchdog;
247: mctx->max_seq_rejects = bctx->max_seq_rejects;
248: switch (bctx->scale_type) {
249: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
250: mctx->sigma = bctx->sigma;
251: break;
252: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
253: MatCopy(bctx->D, mctx->D, SAME_NONZERO_PATTERN);
254: break;
255: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
256: mctx->sigma = 1.0;
257: break;
258: default:
259: break;
260: }
261: return 0;
262: }
264: /*------------------------------------------------------------*/
266: static PetscErrorCode MatReset_LMVMDFP(Mat B, PetscBool destructive)
267: {
268: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
269: Mat_SymBrdn *ldfp = (Mat_SymBrdn*)lmvm->ctx;
270: Mat_LMVM *dbase;
271: Mat_DiagBrdn *dctx;
273: ldfp->watchdog = 0;
274: ldfp->needQ = PETSC_TRUE;
275: if (ldfp->allocated) {
276: if (destructive) {
277: VecDestroy(&ldfp->work);
278: PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
279: VecDestroyVecs(lmvm->m, &ldfp->Q);
280: switch (ldfp->scale_type) {
281: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
282: MatLMVMReset(ldfp->D, PETSC_TRUE);
283: break;
284: default:
285: break;
286: }
287: ldfp->allocated = PETSC_FALSE;
288: } else {
289: switch (ldfp->scale_type) {
290: case MAT_LMVM_SYMBROYDEN_SCALE_SCALAR:
291: ldfp->sigma = ldfp->delta;
292: break;
293: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
294: MatLMVMReset(ldfp->D, PETSC_FALSE);
295: dbase = (Mat_LMVM*)ldfp->D->data;
296: dctx = (Mat_DiagBrdn*)dbase->ctx;
297: VecSet(dctx->invD, ldfp->delta);
298: break;
299: case MAT_LMVM_SYMBROYDEN_SCALE_NONE:
300: ldfp->sigma = 1.0;
301: break;
302: default:
303: break;
304: }
305: }
306: }
307: MatReset_LMVM(B, destructive);
308: return 0;
309: }
311: /*------------------------------------------------------------*/
313: static PetscErrorCode MatAllocate_LMVMDFP(Mat B, Vec X, Vec F)
314: {
315: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
316: Mat_SymBrdn *ldfp = (Mat_SymBrdn*)lmvm->ctx;
318: MatAllocate_LMVM(B, X, F);
319: if (!ldfp->allocated) {
320: VecDuplicate(X, &ldfp->work);
321: PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
322: if (lmvm->m > 0) {
323: VecDuplicateVecs(X, lmvm->m, &ldfp->Q);
324: }
325: switch (ldfp->scale_type) {
326: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
327: MatLMVMAllocate(ldfp->D, X, F);
328: break;
329: default:
330: break;
331: }
332: ldfp->allocated = PETSC_TRUE;
333: }
334: return 0;
335: }
337: /*------------------------------------------------------------*/
339: static PetscErrorCode MatDestroy_LMVMDFP(Mat B)
340: {
341: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
342: Mat_SymBrdn *ldfp = (Mat_SymBrdn*)lmvm->ctx;
344: if (ldfp->allocated) {
345: VecDestroy(&ldfp->work);
346: PetscFree4(ldfp->ytq, ldfp->yts, ldfp->yty, ldfp->sts);
347: VecDestroyVecs(lmvm->m, &ldfp->Q);
348: ldfp->allocated = PETSC_FALSE;
349: }
350: MatDestroy(&ldfp->D);
351: PetscFree(lmvm->ctx);
352: MatDestroy_LMVM(B);
353: return 0;
354: }
356: /*------------------------------------------------------------*/
358: static PetscErrorCode MatSetUp_LMVMDFP(Mat B)
359: {
360: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
361: Mat_SymBrdn *ldfp = (Mat_SymBrdn*)lmvm->ctx;
362: PetscInt n, N;
364: MatSetUp_LMVM(B);
365: if (!ldfp->allocated) {
366: VecDuplicate(lmvm->Xprev, &ldfp->work);
367: PetscMalloc4(lmvm->m, &ldfp->ytq, lmvm->m, &ldfp->yts, lmvm->m, &ldfp->yty, lmvm->m, &ldfp->sts);
368: if (lmvm->m > 0) {
369: VecDuplicateVecs(lmvm->Xprev, lmvm->m, &ldfp->Q);
370: }
371: switch (ldfp->scale_type) {
372: case MAT_LMVM_SYMBROYDEN_SCALE_DIAGONAL:
373: MatGetLocalSize(B, &n, &n);
374: MatGetSize(B, &N, &N);
375: MatSetSizes(ldfp->D, n, n, N, N);
376: MatSetUp(ldfp->D);
377: break;
378: default:
379: break;
380: }
381: ldfp->allocated = PETSC_TRUE;
382: }
383: return 0;
384: }
386: /*------------------------------------------------------------*/
388: static PetscErrorCode MatSetFromOptions_LMVMDFP(PetscOptionItems *PetscOptionsObject, Mat B)
389: {
390: MatSetFromOptions_LMVM(PetscOptionsObject, B);
391: PetscOptionsHead(PetscOptionsObject,"DFP method for approximating SPD Jacobian actions (MATLMVMDFP)");
392: MatSetFromOptions_LMVMSymBrdn_Private(PetscOptionsObject, B);
393: PetscOptionsTail();
394: return 0;
395: }
397: /*------------------------------------------------------------*/
399: PetscErrorCode MatCreate_LMVMDFP(Mat B)
400: {
401: Mat_LMVM *lmvm;
402: Mat_SymBrdn *ldfp;
404: MatCreate_LMVMSymBrdn(B);
405: PetscObjectChangeTypeName((PetscObject)B, MATLMVMDFP);
406: B->ops->setup = MatSetUp_LMVMDFP;
407: B->ops->destroy = MatDestroy_LMVMDFP;
408: B->ops->setfromoptions = MatSetFromOptions_LMVMDFP;
409: B->ops->solve = MatSolve_LMVMDFP;
411: lmvm = (Mat_LMVM*)B->data;
412: lmvm->ops->allocate = MatAllocate_LMVMDFP;
413: lmvm->ops->reset = MatReset_LMVMDFP;
414: lmvm->ops->update = MatUpdate_LMVMDFP;
415: lmvm->ops->mult = MatMult_LMVMDFP;
416: lmvm->ops->copy = MatCopy_LMVMDFP;
418: ldfp = (Mat_SymBrdn*)lmvm->ctx;
419: ldfp->needP = PETSC_FALSE;
420: ldfp->phi = 1.0;
421: return 0;
422: }
424: /*------------------------------------------------------------*/
426: /*@
427: MatCreateLMVMDFP - Creates a limited-memory Davidon-Fletcher-Powell (DFP) matrix
428: used for approximating Jacobians. L-DFP is symmetric positive-definite by
429: construction, and is the dual of L-BFGS where Y and S vectors swap roles.
431: The provided local and global sizes must match the solution and function vectors
432: used with MatLMVMUpdate() and MatSolve(). The resulting L-DFP matrix will have
433: storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
434: parallel. To use the L-DFP matrix with other vector types, the matrix must be
435: created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
436: This ensures that the internal storage and work vectors are duplicated from the
437: correct type of vector.
439: Collective
441: Input Parameters:
442: + comm - MPI communicator, set to PETSC_COMM_SELF
443: . n - number of local rows for storage vectors
444: - N - global size of the storage vectors
446: Output Parameter:
447: . B - the matrix
449: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
450: paradigm instead of this routine directly.
452: Options Database Keys:
453: + -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
454: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
455: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
456: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
457: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
458: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
459: - -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
461: Level: intermediate
463: .seealso: MatCreate(), MATLMVM, MATLMVMDFP, MatCreateLMVMBFGS(), MatCreateLMVMSR1(),
464: MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn(), MatCreateLMVMSymBrdn()
465: @*/
466: PetscErrorCode MatCreateLMVMDFP(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
467: {
468: MatCreate(comm, B);
469: MatSetSizes(*B, n, n, N, N);
470: MatSetType(*B, MATLMVMDFP);
471: MatSetUp(*B);
472: return 0;
473: }