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