Actual source code: symbadbrdn.c
petsc-3.11.4 2019-09-28
1: #include <../src/ksp/ksp/utils/lmvm/symbrdn/symbrdn.h>
2: #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>
4: /*------------------------------------------------------------*/
6: static PetscErrorCode MatSolve_LMVMSymBadBrdn(Mat B, Vec F, Vec dX)
7: {
8: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
9: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
10: PetscErrorCode ierr;
11: PetscInt i, j;
12: PetscScalar yjtqi, sjtyi, wtyi, ytx, stf, wtf, ytq;
13:
15: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
16: if (lsb->phi == 0.0) {
17: MatSolve_LMVMBFGS(B, F, dX);
18: return(0);
19: }
20: if (lsb->phi == 1.0) {
21: MatSolve_LMVMDFP(B, F, dX);
22: return(0);
23: }
24:
25: VecCheckSameSize(F, 2, dX, 3);
26: VecCheckMatCompatible(B, dX, 3, F, 2);
27:
28: if (lsb->needQ) {
29: /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
30: for (i = 0; i <= lmvm->k; ++i) {
31: MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]);
32: for (j = 0; j <= i-1; ++j) {
33: /* Compute the necessary dot products */
34: VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi);
35: VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
36: VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi);
37: VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
38: /* Compute the pure DFP component of the inverse application*/
39: VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi)/lsb->ytq[j], PetscRealPart(sjtyi)/lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]);
40: /* Tack on the convexly scaled extras to the inverse application*/
41: if (lsb->psi[j] > 0.0) {
42: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]);
43: VecDot(lsb->work, lmvm->Y[i], &wtyi);
44: VecAXPY(lsb->Q[i], lsb->phi*lsb->ytq[j]*PetscRealPart(wtyi), lsb->work);
45: }
46: }
47: VecDot(lmvm->Y[i], lsb->Q[i], &ytq);
48: lsb->ytq[i] = PetscRealPart(ytq);
49: }
50: lsb->needQ = PETSC_FALSE;
51: }
52:
53: /* Start the outer iterations for ((B^{-1}) * dX) */
54: MatSymBrdnApplyJ0Inv(B, F, dX);
55: for (i = 0; i <= lmvm->k; ++i) {
56: /* Compute the necessary dot products -- store yTs and yTp for inner iterations later */
57: VecDotBegin(lmvm->Y[i], dX, &ytx);
58: VecDotBegin(lmvm->S[i], F, &stf);
59: VecDotEnd(lmvm->Y[i], dX, &ytx);
60: VecDotEnd(lmvm->S[i], F, &stf);
61: /* Compute the pure DFP component */
62: VecAXPBYPCZ(dX, -PetscRealPart(ytx)/lsb->ytq[i], PetscRealPart(stf)/lsb->yts[i], 1.0, lsb->Q[i], lmvm->S[i]);
63: /* Tack on the convexly scaled extras */
64: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->ytq[i], 0.0, lmvm->S[i], lsb->Q[i]);
65: VecDot(lsb->work, F, &wtf);
66: VecAXPY(dX, lsb->phi*lsb->ytq[i]*PetscRealPart(wtf), lsb->work);
67: }
69: return(0);
70: }
72: /*------------------------------------------------------------*/
74: static PetscErrorCode MatMult_LMVMSymBadBrdn(Mat B, Vec X, Vec Z)
75: {
76: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
77: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
78: PetscErrorCode ierr;
79: PetscInt i, j;
80: PetscReal numer;
81: PetscScalar sjtpi, sjtyi, yjtsi, yjtqi, wtsi, wtyi, stz, ytx, ytq, wtx, stp;
82:
83:
85: /* Efficient shortcuts for pure BFGS and pure DFP configurations */
86: if (lsb->phi == 0.0) {
87: MatMult_LMVMBFGS(B, X, Z);
88: return(0);
89: }
90: if (lsb->phi == 1.0) {
91: MatMult_LMVMDFP(B, X, Z);
92: return(0);
93: }
94:
95: VecCheckSameSize(X, 2, Z, 3);
96: VecCheckMatCompatible(B, X, 2, Z, 3);
97:
98: if (lsb->needQ) {
99: /* Start the loop for (Q[k] = (B_k)^{-1} * Y[k]) */
100: for (i = 0; i <= lmvm->k; ++i) {
101: MatSymBrdnApplyJ0Inv(B, lmvm->Y[i], lsb->Q[i]);
102: for (j = 0; j <= i-1; ++j) {
103: /* Compute the necessary dot products */
104: VecDotBegin(lmvm->Y[j], lsb->Q[i], &yjtqi);
105: VecDotBegin(lmvm->S[j], lmvm->Y[i], &sjtyi);
106: VecDotEnd(lmvm->Y[j], lsb->Q[i], &yjtqi);
107: VecDotEnd(lmvm->S[j], lmvm->Y[i], &sjtyi);
108: /* Compute the pure DFP component of the inverse application*/
109: VecAXPBYPCZ(lsb->Q[i], -PetscRealPart(yjtqi)/lsb->ytq[j], PetscRealPart(sjtyi)/lsb->yts[j], 1.0, lsb->Q[j], lmvm->S[j]);
110: /* Tack on the convexly scaled extras to the inverse application*/
111: if (lsb->psi[j] > 0.0) {
112: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->ytq[j], 0.0, lmvm->S[j], lsb->Q[j]);
113: VecDot(lsb->work, lmvm->Y[i], &wtyi);
114: VecAXPY(lsb->Q[i], lsb->phi*lsb->ytq[j]*PetscRealPart(wtyi), lsb->work);
115: }
116: }
117: VecDot(lmvm->Y[i], lsb->Q[i], &ytq);
118: lsb->ytq[i] = PetscRealPart(ytq);
119: }
120: lsb->needQ = PETSC_FALSE;
121: }
122: if (lsb->needP) {
123: /* Start the loop for (P[k] = (B_k) * S[k]) */
124: for (i = 0; i <= lmvm->k; ++i) {
125: MatSymBrdnApplyJ0Fwd(B, lmvm->S[i], lsb->P[i]);
126: for (j = 0; j <= i-1; ++j) {
127: /* Compute the necessary dot products */
128: VecDotBegin(lmvm->S[j], lsb->P[i], &sjtpi);
129: VecDotBegin(lmvm->Y[j], lmvm->S[i], &yjtsi);
130: VecDotEnd(lmvm->S[j], lsb->P[i], &sjtpi);
131: VecDotEnd(lmvm->Y[j], lmvm->S[i], &yjtsi);
132: /* Compute the pure BFGS component of the forward product */
133: VecAXPBYPCZ(lsb->P[i], -PetscRealPart(sjtpi)/lsb->stp[j], PetscRealPart(yjtsi)/lsb->yts[j], 1.0, lsb->P[j], lmvm->Y[j]);
134: /* Tack on the convexly scaled extras to the forward product */
135: if (lsb->phi > 0.0) {
136: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[j], -1.0/lsb->stp[j], 0.0, lmvm->Y[j], lsb->P[j]);
137: VecDot(lsb->work, lmvm->S[i], &wtsi);
138: VecAXPY(lsb->P[i], lsb->psi[j]*lsb->stp[j]*PetscRealPart(wtsi), lsb->work);
139: }
140: }
141: VecDot(lmvm->S[i], lsb->P[i], &stp);
142: lsb->stp[i] = PetscRealPart(stp);
143: if (lsb->phi == 1.0) {
144: lsb->psi[i] = 0.0;
145: } else if (lsb->phi == 0.0) {
146: lsb->psi[i] = 1.0;
147: } else {
148: numer = (1.0 - lsb->phi)*lsb->yts[i]*lsb->yts[i];
149: lsb->psi[i] = numer / (numer + (lsb->phi*lsb->ytq[i]*lsb->stp[i]));
150: }
151: }
152: lsb->needP = PETSC_FALSE;
153: }
154:
155: /* Start the outer iterations for (B * X) */
156: MatSymBrdnApplyJ0Fwd(B, X, Z);
157: for (i = 0; i <= lmvm->k; ++i) {
158: /* Compute the necessary dot products */
159: VecDotBegin(lmvm->S[i], Z, &stz);
160: VecDotBegin(lmvm->Y[i], X, &ytx);
161: VecDotEnd(lmvm->S[i], Z, &stz);
162: VecDotEnd(lmvm->Y[i], X, &ytx);
163: /* Compute the pure BFGS component */
164: VecAXPBYPCZ(Z, -PetscRealPart(stz)/lsb->stp[i], PetscRealPart(ytx)/lsb->yts[i], 1.0, lsb->P[i], lmvm->Y[i]);
165: /* Tack on the convexly scaled extras */
166: VecAXPBYPCZ(lsb->work, 1.0/lsb->yts[i], -1.0/lsb->stp[i], 0.0, lmvm->Y[i], lsb->P[i]);
167: VecDot(lsb->work, X, &wtx);
168: VecAXPY(Z, lsb->psi[i]*lsb->stp[i]*PetscRealPart(wtx), lsb->work);
169: }
170: return(0);
171: }
173: /*------------------------------------------------------------*/
175: static PetscErrorCode MatSetFromOptions_LMVMSymBadBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
176: {
177: Mat_LMVM *lmvm = (Mat_LMVM*)B->data;
178: Mat_SymBrdn *lsb = (Mat_SymBrdn*)lmvm->ctx;
179: Mat_LMVM *dbase;
180: Mat_DiagBrdn *dctx;
181: PetscErrorCode ierr;
184: MatSetFromOptions_LMVMSymBrdn(PetscOptionsObject, B);
185: if (lsb->scale_type == SYMBRDN_SCALE_DIAG) {
186: dbase = (Mat_LMVM*)lsb->D->data;
187: dctx = (Mat_DiagBrdn*)dbase->ctx;
188: dctx->forward = PETSC_FALSE;
189: }
190: if (Scale_Table[0][0] == '#') SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Scale_Table starts with illegal hashtag character!"); /* Dummy use of Scale_Table to prevent unused variable warnings in this translation unit */
191: return(0);
192: }
194: /*------------------------------------------------------------*/
196: PetscErrorCode MatCreate_LMVMSymBadBrdn(Mat B)
197: {
198: Mat_LMVM *lmvm;
199: PetscErrorCode ierr;
202: MatCreate_LMVMSymBrdn(B);
203: PetscObjectChangeTypeName((PetscObject)B, MATLMVMSYMBADBRDN);
204: B->ops->setfromoptions = MatSetFromOptions_LMVMSymBadBrdn;
205: B->ops->solve = MatSolve_LMVMSymBadBrdn;
206:
207: lmvm = (Mat_LMVM*)B->data;
208: lmvm->ops->mult = MatMult_LMVMSymBadBrdn;
209: return(0);
210: }
212: /*------------------------------------------------------------*/
214: /*@
215: MatCreateLMVMSymBadBrdn - Creates a limited-memory Symmetric "Bad" Broyden-type matrix used
216: for approximating Jacobians. L-SymBadBrdn is a convex combination of L-DFP and
217: L-BFGS such that SymBrdn^{-1} = (1 - phi)*BFGS^{-1} + phi*DFP^{-1}. The combination factor
218: phi is restricted to the range [0, 1], where the L-SymBadBrdn matrix is guaranteed
219: to be symmetric positive-definite. Note that this combination is on the inverses and not
220: on the forwards. For forward convex combinations, use the L-SymBrdn matrix.
222: The provided local and global sizes must match the solution and function vectors
223: used with MatLMVMUpdate() and MatSolve(). The resulting L-SymBrdn matrix will have
224: storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
225: parallel. To use the L-SymBrdn matrix with other vector types, the matrix must be
226: created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
227: This ensures that the internal storage and work vectors are duplicated from the
228: correct type of vector.
230: Collective on MPI_Comm
232: Input Parameters:
233: + comm - MPI communicator, set to PETSC_COMM_SELF
234: . n - number of local rows for storage vectors
235: - N - global size of the storage vectors
237: Output Parameter:
238: . B - the matrix
240: It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
241: paradigm instead of this routine directly.
243: Options Database Keys:
244: . -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
245: . -mat_lmvm_phi - (developer) convex ratio between BFGS and DFP components of the update
246: . -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
247: . -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
248: . -mat_lmvm_rho - (developer) update limiter for the J0 scaling
249: . -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
250: . -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
251: . -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling
253: Level: intermediate
255: .seealso: MatCreate(), MATLMVM, MATLMVMSYMBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
256: MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMBadBrdn()
257: @*/
258: PetscErrorCode MatCreateLMVMSymBadBrdn(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
259: {
260: PetscErrorCode ierr;
263: MatCreate(comm, B);
264: MatSetSizes(*B, n, n, N, N);
265: MatSetType(*B, MATLMVMSYMBADBRDN);
266: MatSetUp(*B);
267: return(0);
268: }