Actual source code: badbrdn.c

petsc-3.11.2 2019-05-18
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;
 12: } Mat_BadBrdn;

 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;
 37:   Mat_BadBrdn       *lbb = (Mat_BadBrdn*)lmvm->ctx;
 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;
 90:   Mat_BadBrdn       *lbb = (Mat_BadBrdn*)lmvm->ctx;
 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;
124:   Mat_BadBrdn       *lbb = (Mat_BadBrdn*)lmvm->ctx;
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;
166:   Mat_BadBrdn       *bctx = (Mat_BadBrdn*)bdata->ctx;
167:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
168:   Mat_BadBrdn       *mctx = (Mat_BadBrdn*)mdata->ctx;
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;
189:   Mat_BadBrdn       *lbb = (Mat_BadBrdn*)lmvm->ctx;
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;
209:   Mat_BadBrdn          *lbb = (Mat_BadBrdn*)lmvm->ctx;
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: /*------------------------------------------------------------*/

227: static PetscErrorCode MatDestroy_LMVMBadBrdn(Mat B)
228: {
229:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
230:   Mat_BadBrdn       *lbb = (Mat_BadBrdn*)lmvm->ctx;
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: /*------------------------------------------------------------*/

247: static PetscErrorCode MatSetUp_LMVMBadBrdn(Mat B)
248: {
249:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
250:   Mat_BadBrdn       *lbb = (Mat_BadBrdn*)lmvm->ctx;
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: /*------------------------------------------------------------*/

268: PetscErrorCode MatCreate_LMVMBadBrdn(Mat B)
269: {
270:   Mat_LMVM          *lmvm;
271:   Mat_BadBrdn       *lbb;
272:   PetscErrorCode    ierr;

275:   MatCreate_LMVM(B);
276:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMBADBRDN);
277:   B->ops->setup = MatSetUp_LMVMBadBrdn;
278:   B->ops->destroy = MatDestroy_LMVMBadBrdn;
279:   B->ops->solve = MatSolve_LMVMBadBrdn;

281:   lmvm = (Mat_LMVM*)B->data;
282:   lmvm->square = PETSC_TRUE;
283:   lmvm->ops->allocate = MatAllocate_LMVMBadBrdn;
284:   lmvm->ops->reset = MatReset_LMVMBadBrdn;
285:   lmvm->ops->mult = MatMult_LMVMBadBrdn;
286:   lmvm->ops->update = MatUpdate_LMVMBadBrdn;
287:   lmvm->ops->copy = MatCopy_LMVMBadBrdn;

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: /*@
299:    MatCreateLMVMBadBrdn - Creates a limited-memory modified (aka "bad") Broyden-type 
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()
322:    paradigm instead of this routine directly.

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);
345:   MatSetType(*B, MATLMVMBADBRDN);
346:   MatSetUp(*B);
347:   return(0);
348: }