petsc-3.11.3 2019-06-26
Report Typos and Errors
```  1:  #include <../src/ksp/ksp/utils/lmvm/lmvm.h>

3: /*
4:   Limited-memory "good" Broyden's method for approximating the inverse of
5:   a Jacobian.
6: */

8: typedef struct {
9:   Vec *P, *Q;
10:   PetscBool allocated, needP, needQ;
11:   PetscReal *yty, *yts;

14: /*------------------------------------------------------------*/

16: /*
17:   The solution method is the matrix-free implementation of the inverse Hessian in
18:   Equation 6 on page 312 of Griewank "Broyden Updating, The Good and The Bad!"
19:   (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).
20:
21:   Q[i] = (B_i)^{-1}*S[i] terms are computed ahead of time whenever
22:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
23:   repeated calls of MatSolve without incurring redundant computation.
24:
25:   dX <- J0^{-1} * F
26:
27:   for i=0,1,2,...,k
28:     # Q[i] = (B_i)^{-1} * Y[i]
29:     tau = (Y[i]^T F) / (Y[i]^T Y[i])
30:     dX <- dX + (tau * (S[i] - Q[i]))
31:   end
32:  */

34: static PetscErrorCode MatSolve_LMVMBadBrdn(Mat B, Vec F, Vec dX)
35: {
36:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
38:   PetscErrorCode    ierr;
39:   PetscInt          i, j;
40:   PetscScalar       yjtyi, ytf;
41:
43:   VecCheckSameSize(F, 2, dX, 3);
44:   VecCheckMatCompatible(B, dX, 3, F, 2);
45:
46:   if (lbb->needQ) {
47:     /* Pre-compute (Q[i] = (B_i)^{-1} * Y[i]) */
48:     for (i = 0; i <= lmvm->k; ++i) {
49:       MatLMVMApplyJ0Inv(B, lmvm->Y[i], lbb->Q[i]);
50:       for (j = 0; j <= i-1; ++j) {
51:         VecDot(lmvm->Y[j], lmvm->Y[i], &yjtyi);
52:         VecAXPBYPCZ(lbb->Q[i], PetscRealPart(yjtyi)/lbb->yty[j], -PetscRealPart(yjtyi)/lbb->yty[j], 1.0, lmvm->S[j], lbb->Q[j]);
53:       }
54:     }
55:     lbb->needQ = PETSC_FALSE;
56:   }
57:
58:   MatLMVMApplyJ0Inv(B, F, dX);
59:   for (i = 0; i <= lmvm->k-1; ++i) {
60:     VecDot(lmvm->Y[i], F, &ytf);
61:     VecAXPBYPCZ(dX, PetscRealPart(ytf)/lbb->yty[i], -PetscRealPart(ytf)/lbb->yty[i], 1.0, lmvm->S[i], lbb->Q[i]);
62:   }
63:   return(0);
64: }

66: /*------------------------------------------------------------*/

68: /*
69:   The forward product is the matrix-free implementation of the direct update in
70:   Equation 6 on page 302 of Griewank "Broyden Updating, The Good and The Bad!"
71:   (http://www.emis.ams.org/journals/DMJDMV/vol-ismp/45_griewank-andreas-broyden.pdf).
72:
73:   P[i] = (B_i)*S[i] terms are computed ahead of time whenever
74:   the matrix is updated with a new (S[i], Y[i]) pair. This allows
75:   repeated calls of MatMult inside KSP solvers without unnecessarily
76:   recomputing P[i] terms in expensive nested-loops.
77:
78:   Z <- J0 * X
79:
80:   for i=0,1,2,...,k
81:     # P[i] = B_i * S[i]
82:     tau = (Y[i]^T X) / (Y[i]^T S[i])
83:     dX <- dX + (tau * (Y[i] - P[i]))
84:   end
85:  */

87: static PetscErrorCode MatMult_LMVMBadBrdn(Mat B, Vec X, Vec Z)
88: {
89:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
91:   PetscErrorCode    ierr;
92:   PetscInt          i, j;
93:   PetscScalar       yjtsi, ytx;
94:
96:   VecCheckSameSize(X, 2, Z, 3);
97:   VecCheckMatCompatible(B, X, 2, Z, 3);
98:
99:   if (lbb->needP) {
100:     /* Pre-compute (P[i] = (B_i) * S[i]) */
101:     for (i = 0; i <= lmvm->k; ++i) {
102:       MatLMVMApplyJ0Fwd(B, lmvm->S[i], lbb->P[i]);
103:       for (j = 0; j <= i-1; ++j) {
104:         VecDot(lmvm->Y[j], lmvm->S[i], &yjtsi);
105:         VecAXPBYPCZ(lbb->P[i], PetscRealPart(yjtsi)/lbb->yts[j], -PetscRealPart(yjtsi)/lbb->yts[j], 1.0, lmvm->Y[j], lbb->P[j]);
106:       }
107:     }
108:     lbb->needP = PETSC_FALSE;
109:   }
110:
111:   MatLMVMApplyJ0Fwd(B, X, Z);
112:   for (i = 0; i <= lmvm->k-1; ++i) {
113:     VecDot(lmvm->Y[i], X, &ytx);
114:     VecAXPBYPCZ(Z, PetscRealPart(ytx)/lbb->yts[i], -PetscRealPart(ytx)/lbb->yts[i], 1.0, lmvm->Y[i], lbb->P[i]);
115:   }
116:   return(0);
117: }

119: /*------------------------------------------------------------*/

121: static PetscErrorCode MatUpdate_LMVMBadBrdn(Mat B, Vec X, Vec F)
122: {
123:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
125:   PetscErrorCode    ierr;
126:   PetscInt          old_k, i;
127:   PetscScalar       yty, yts;

130:   if (!lmvm->m) return(0);
131:   if (lmvm->prev_set) {
132:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
133:     VecAYPX(lmvm->Xprev, -1.0, X);
134:     VecAYPX(lmvm->Fprev, -1.0, F);
135:     /* Accept the update */
136:     lbb->needP = lbb->needQ = PETSC_TRUE;
137:     old_k = lmvm->k;
138:     MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
139:     /* If we hit the memory limit, shift the yty and yts arrays */
140:     if (old_k == lmvm->k) {
141:       for (i = 0; i <= lmvm->k-1; ++i) {
142:         lbb->yty[i] = lbb->yty[i+1];
143:         lbb->yts[i] = lbb->yts[i+1];
144:       }
145:     }
146:     /* Accumulate the latest yTy and yTs dot products */
147:     VecDotBegin(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
148:     VecDotBegin(lmvm->Y[lmvm->k], lmvm->S[lmvm->k], &yts);
149:     VecDotEnd(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
150:     VecDotEnd(lmvm->Y[lmvm->k], lmvm->S[lmvm->k], &yts);
151:     lbb->yty[lmvm->k] = PetscRealPart(yty);
152:     lbb->yts[lmvm->k] = PetscRealPart(yts);
153:   }
154:   /* Save the solution and function to be used in the next update */
155:   VecCopy(X, lmvm->Xprev);
156:   VecCopy(F, lmvm->Fprev);
157:   lmvm->prev_set = PETSC_TRUE;
158:   return(0);
159: }

161: /*------------------------------------------------------------*/

163: static PetscErrorCode MatCopy_LMVMBadBrdn(Mat B, Mat M, MatStructure str)
164: {
165:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
167:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
169:   PetscErrorCode    ierr;
170:   PetscInt          i;

173:   mctx->needP = bctx->needP;
174:   mctx->needQ = bctx->needQ;
175:   for (i=0; i<=bdata->k; ++i) {
176:     mctx->yty[i] = bctx->yty[i];
177:     mctx->yts[i] = bctx->yts[i];
178:     VecCopy(bctx->P[i], mctx->P[i]);
179:     VecCopy(bctx->Q[i], mctx->Q[i]);
180:   }
181:   return(0);
182: }

184: /*------------------------------------------------------------*/

186: static PetscErrorCode MatReset_LMVMBadBrdn(Mat B, PetscBool destructive)
187: {
188:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
190:   PetscErrorCode    ierr;
191:
193:   lbb->needP = lbb->needQ = PETSC_TRUE;
194:   if (destructive && lbb->allocated) {
195:     PetscFree2(lbb->yty, lbb->yts);
196:     VecDestroyVecs(lmvm->m, &lbb->P);
197:     VecDestroyVecs(lmvm->m, &lbb->Q);
198:     lbb->allocated = PETSC_FALSE;
199:   }
200:   MatReset_LMVM(B, destructive);
201:   return(0);
202: }

204: /*------------------------------------------------------------*/

206: static PetscErrorCode MatAllocate_LMVMBadBrdn(Mat B, Vec X, Vec F)
207: {
208:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
210:   PetscErrorCode    ierr;
211:
213:   MatAllocate_LMVM(B, X, F);
214:   if (!lbb->allocated) {
215:     PetscMalloc2(lmvm->m, &lbb->yty, lmvm->m, &lbb->yts);
216:     if (lmvm->m > 0) {
217:       VecDuplicateVecs(X, lmvm->m, &lbb->P);
218:       VecDuplicateVecs(X, lmvm->m, &lbb->Q);
219:     }
220:     lbb->allocated = PETSC_TRUE;
221:   }
222:   return(0);
223: }

225: /*------------------------------------------------------------*/

228: {
229:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
231:   PetscErrorCode    ierr;

234:   if (lbb->allocated) {
235:     PetscFree2(lbb->yty, lbb->yts);
236:     VecDestroyVecs(lmvm->m, &lbb->P);
237:     VecDestroyVecs(lmvm->m, &lbb->Q);
238:     lbb->allocated = PETSC_FALSE;
239:   }
240:   PetscFree(lmvm->ctx);
241:   MatDestroy_LMVM(B);
242:   return(0);
243: }

245: /*------------------------------------------------------------*/

248: {
249:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
251:   PetscErrorCode    ierr;
252:
254:   MatSetUp_LMVM(B);
255:   if (!lbb->allocated) {
256:     PetscMalloc2(lmvm->m, &lbb->yty, lmvm->m, &lbb->yts);
257:     if (lmvm->m > 0) {
258:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbb->P);
259:       VecDuplicateVecs(lmvm->Xprev, lmvm->m, &lbb->Q);
260:     }
261:     lbb->allocated = PETSC_TRUE;
262:   }
263:   return(0);
264: }

266: /*------------------------------------------------------------*/

269: {
270:   Mat_LMVM          *lmvm;
272:   PetscErrorCode    ierr;

275:   MatCreate_LMVM(B);

281:   lmvm = (Mat_LMVM*)B->data;
282:   lmvm->square = PETSC_TRUE;

289:   PetscNewLog(B, &lbb);
290:   lmvm->ctx = (void*)lbb;
291:   lbb->allocated = PETSC_FALSE;
292:   lbb->needP = lbb->needQ = PETSC_TRUE;
293:   return(0);
294: }

296: /*------------------------------------------------------------*/

298: /*@
300:    approximation matrix used for a Jacobian. L-BadBrdn is not guaranteed to be
301:    symmetric or positive-definite.
302:
303:    The provided local and global sizes must match the solution and function vectors
304:    used with MatLMVMUpdate() and MatSolve(). The resulting L-BadBrdn matrix will have
305:    storage vectors allocated with VecCreateSeq() in serial and VecCreateMPI() in
306:    parallel. To use the L-BadBrdn matrix with other vector types, the matrix must be
307:    created using MatCreate() and MatSetType(), followed by MatLMVMAllocate().
308:    This ensures that the internal storage and work vectors are duplicated from the
309:    correct type of vector.

311:    Collective on MPI_Comm

313:    Input Parameters:
314: +  comm - MPI communicator, set to PETSC_COMM_SELF
315: .  n - number of local rows for storage vectors
316: -  N - global size of the storage vectors

318:    Output Parameter:
319: .  B - the matrix

321:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()

324:    Options Database Keys:
325: .   -mat_lmvm_num_vecs - maximum number of correction vectors (i.e.: updates) stored
326: .   -mat_lmvm_scale_type - (developer) type of scaling applied to J0 (none, scalar, diagonal)
327: .   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
328: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
329: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
330: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
331: .   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling

333:    Level: intermediate

335: .seealso: MatCreate(), MATLMVM, MATLMVMBADBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
336:           MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMSymBrdn()
337: @*/
338: PetscErrorCode MatCreateLMVMBadBrdn(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
339: {
340:   PetscErrorCode    ierr;
341:
343:   MatCreate(comm, B);
344:   MatSetSizes(*B, n, n, N, N);