Actual source code: diagbrdn.c

petsc-master 2019-08-18
Report Typos and Errors
  1:  #include <../src/ksp/ksp/utils/lmvm/diagbrdn/diagbrdn.h>

  3: /*------------------------------------------------------------*/

  5: static PetscErrorCode MatSolve_DiagBrdn(Mat B, Vec F, Vec dX)
  6: {
  7:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
  8:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
  9:   PetscErrorCode    ierr;

 12:   VecCheckSameSize(F, 2, dX, 3);
 13:   VecCheckMatCompatible(B, dX, 3, F, 2);
 14:   VecPointwiseMult(dX, ldb->invD, F);
 15:   return(0);
 16: }

 18: /*------------------------------------------------------------*/

 20: static PetscErrorCode MatMult_DiagBrdn(Mat B, Vec X, Vec Z)
 21: {
 22:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
 23:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
 24:   PetscErrorCode    ierr;

 27:   VecCheckSameSize(X, 2, Z, 3);
 28:   VecCheckMatCompatible(B, X, 2, Z, 3);
 29:   VecPointwiseDivide(Z, X, ldb->invD);
 30:   return(0);
 31: }

 33: /*------------------------------------------------------------*/

 35: static PetscErrorCode MatUpdate_DiagBrdn(Mat B, Vec X, Vec F)
 36: {
 37:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
 38:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
 39:   PetscErrorCode    ierr;
 40:   PetscInt          old_k, i, start;
 41:   PetscScalar       yty, ststmp, curvature, ytDy, stDs, ytDs;
 42:   PetscReal         curvtol, sigma, yy_sum, ss_sum, ys_sum, denom;

 45:   if (!lmvm->m) return(0);
 46:   if (lmvm->prev_set) {
 47:     /* Compute the new (S = X - Xprev) and (Y = F - Fprev) vectors */
 48:     VecAYPX(lmvm->Xprev, -1.0, X);
 49:     VecAYPX(lmvm->Fprev, -1.0, F);
 50:     /* Compute tolerance for accepting the update */
 51:     VecDotBegin(lmvm->Xprev, lmvm->Fprev, &curvature);
 52:     VecDotBegin(lmvm->Xprev, lmvm->Xprev, &ststmp);
 53:     VecDotEnd(lmvm->Xprev, lmvm->Fprev, &curvature);
 54:     VecDotEnd(lmvm->Xprev, lmvm->Xprev, &ststmp);
 55:     if (PetscRealPart(ststmp) < lmvm->eps){
 56:       curvtol = 0.0;
 57:     } else {
 58:       curvtol = lmvm->eps * PetscRealPart(ststmp);
 59:     }
 60:     /* Test the curvature for the update */
 61:     if (PetscRealPart(curvature) > curvtol) {
 62:       /* Update is good so we accept it */
 63:       old_k = lmvm->k;
 64:       MatUpdateKernel_LMVM(B, lmvm->Xprev, lmvm->Fprev);
 65:       /* If we hit the memory limit, shift the yty and yts arrays */
 66:       if (old_k == lmvm->k) {
 67:         for (i = 0; i <= lmvm->k-1; ++i) {
 68:           ldb->yty[i] = ldb->yty[i+1];
 69:           ldb->yts[i] = ldb->yts[i+1];
 70:           ldb->sts[i] = ldb->sts[i+1];
 71:         }
 72:       }
 73:       /* Accept dot products into the history */
 74:       VecDot(lmvm->Y[lmvm->k], lmvm->Y[lmvm->k], &yty);
 75:       ldb->yty[lmvm->k] = PetscRealPart(yty);
 76:       ldb->yts[lmvm->k] = PetscRealPart(curvature);
 77:       ldb->sts[lmvm->k] = PetscRealPart(ststmp);
 78:       if (ldb->forward) {
 79:         /* We are doing diagonal scaling of the forward Hessian B */
 80:         /*  BFGS = DFP = inv(D); */
 81:         VecCopy(ldb->invD, ldb->invDnew);
 82:         VecReciprocal(ldb->invDnew);

 84:         /*  V = y*y */
 85:         VecPointwiseMult(ldb->V, lmvm->Y[lmvm->k], lmvm->Y[lmvm->k]);

 87:         /*  W = inv(D)*s */
 88:         VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->S[lmvm->k]);
 89:         VecDot(ldb->W, lmvm->S[lmvm->k], &stDs);

 91:         /*  Safeguard stDs */
 92:         stDs = PetscMax(PetscRealPart(stDs), ldb->tol);

 94:         if (1.0 != ldb->theta) {
 95:           /*  BFGS portion of the update */
 96:           /*  U = (inv(D)*s)*(inv(D)*s) */
 97:           VecPointwiseMult(ldb->U, ldb->W, ldb->W);

 99:           /*  Assemble */
100:           VecAXPBY(ldb->BFGS, -1.0/stDs, 0.0, ldb->U);
101:         }
102:         if (0.0 != ldb->theta) {
103:           /*  DFP portion of the update */
104:           /*  U = inv(D)*s*y */
105:           VecPointwiseMult(ldb->U, ldb->W, lmvm->Y[lmvm->k]);

107:           /*  Assemble */
108:           VecAXPBY(ldb->DFP, stDs/ldb->yts[lmvm->k], 0.0, ldb->V);
109:           VecAXPY(ldb->DFP, -2.0, ldb->U);
110:         }

112:         if (0.0 == ldb->theta) {
113:           VecAXPY(ldb->invDnew, 1.0, ldb->BFGS);
114:         } else if (1.0 == ldb->theta) {
115:           VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->DFP);
116:         } else {
117:           /*  Broyden update Dkp1 = Dk + (1-theta)*P + theta*Q + y_i^2/yts*/
118:           VecAXPBYPCZ(ldb->invDnew, 1.0-ldb->theta, (ldb->theta)/ldb->yts[lmvm->k], 1.0, ldb->BFGS, ldb->DFP);
119:         }

121:         VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
122:         /*  Obtain inverse and ensure positive definite */
123:         VecReciprocal(ldb->invDnew);
124:         VecAbs(ldb->invDnew);

126:       } else {
127:         /* Inverse Hessian update instead. */
128:         VecCopy(ldb->invD, ldb->invDnew);

130:         /*  V = s*s */
131:         VecPointwiseMult(ldb->V, lmvm->S[lmvm->k], lmvm->S[lmvm->k]);

133:         /*  W = D*y */
134:         VecPointwiseMult(ldb->W, ldb->invDnew, lmvm->Y[lmvm->k]);
135:         VecDot(ldb->W, lmvm->Y[lmvm->k], &ytDy);

137:         /*  Safeguard ytDy */
138:         ytDy = PetscMax(PetscRealPart(ytDy), ldb->tol);

140:         if (1.0 != ldb->theta) {
141:           /*  BFGS portion of the update */
142:           /*  U = s*Dy */
143:           VecPointwiseMult(ldb->U, ldb->W, lmvm->S[lmvm->k]);

145:           /*  Assemble */
146:           VecAXPBY(ldb->BFGS, ytDy/ldb->yts[lmvm->k], 0.0, ldb->V);
147:           VecAXPY(ldb->BFGS, -2.0, ldb->U);
148:         }
149:         if (0.0 != ldb->theta) {
150:           /*  DFP portion of the update */

152:           /*  U = (inv(D)*y)*(inv(D)*y) */
153:           VecPointwiseMult(ldb->U, ldb->W, ldb->W);

155:           /*  Assemble */
156:           VecAXPBY(ldb->DFP, -1.0/ytDy, 0.0, ldb->U);
157:         }

159:         if (0.0 == ldb->theta) {
160:           VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->BFGS);
161:         } else if (1.0 == ldb->theta) {
162:           VecAXPY(ldb->invDnew, 1.0, ldb->DFP);
163:         } else {
164:           /*  Broyden update U=(1-theta)*P + theta*Q */
165:           VecAXPBYPCZ(ldb->invDnew, (1.0-ldb->theta)/ldb->yts[lmvm->k], ldb->theta, 1.0, ldb->BFGS, ldb->DFP);
166:         }
167:         VecAXPY(ldb->invDnew, 1.0/ldb->yts[lmvm->k], ldb->V);
168:         /*  Ensure positive definite */
169:         VecAbs(ldb->invDnew);
170:       }
171:       if (ldb->sigma_hist > 0){
172:         /*  Start with re-scaling on the newly computed diagonal */
173:         if (0.5 == ldb->beta) {
174:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
175:             VecPointwiseMult(ldb->V, lmvm->Y[0], ldb->invDnew);
176:             VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);

178:             VecDotBegin(ldb->V, lmvm->Y[0], &ytDy);
179:             VecDotBegin(ldb->W, lmvm->S[0], &stDs);
180:             VecDotEnd(ldb->V, lmvm->Y[0], &ytDy);
181:             VecDotEnd(ldb->W, lmvm->S[0], &stDs);

183:             ss_sum = PetscRealPart(stDs);
184:             yy_sum = PetscRealPart(ytDy);
185:             ys_sum = ldb->yts[0];
186:           } else {
187:             VecCopy(ldb->invDnew, ldb->U);
188:             VecReciprocal(ldb->U);

190:             /*  Compute summations for scalar scaling */
191:             yy_sum = 0;       /*  No safeguard required */
192:             ys_sum = 0;       /*  No safeguard required */
193:             ss_sum = 0;       /*  No safeguard required */
194:             start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
195:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
196:               VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->U);
197:               VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);

199:               VecDotBegin(ldb->W, lmvm->S[i], &stDs);
200:               VecDotBegin(ldb->V, lmvm->Y[i], &ytDy);
201:               VecDotEnd(ldb->W, lmvm->S[i], &stDs);
202:               VecDotEnd(ldb->V, lmvm->Y[i], &ytDy);

204:               ss_sum += PetscRealPart(stDs);
205:               ys_sum += ldb->yts[i];
206:               yy_sum += PetscRealPart(ytDy);
207:             }
208:           }
209:         } else if (0.0 == ldb->beta) {
210:           if (1 == PetscMin(lmvm->nupdates, ldb->sigma_hist)) {
211:             /*  Compute summations for scalar scaling */
212:             VecPointwiseDivide(ldb->W, lmvm->S[0], ldb->invDnew);

214:             VecDotBegin(ldb->W, lmvm->Y[0], &ytDs);
215:             VecDotBegin(ldb->W, ldb->W, &stDs);
216:             VecDotEnd(ldb->W, lmvm->Y[0], &ytDs);
217:             VecDotEnd(ldb->W, ldb->W, &stDs);

219:             ys_sum = PetscRealPart(ytDs);
220:             ss_sum = PetscRealPart(stDs);
221:             yy_sum = ldb->yty[0];
222:           } else {
223:             VecCopy(ldb->invDnew, ldb->U);
224:             VecReciprocal(ldb->U);

226:             /*  Compute summations for scalar scaling */
227:             yy_sum = 0;       /*  No safeguard required */
228:             ys_sum = 0;       /*  No safeguard required */
229:             ss_sum = 0;       /*  No safeguard required */
230:             start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
231:             for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
232:               VecPointwiseMult(ldb->W, lmvm->S[i], ldb->U);

234:               VecDotBegin(ldb->W, lmvm->Y[i], &ytDs);
235:               VecDotBegin(ldb->W, ldb->W, &stDs);
236:               VecDotEnd(ldb->W, lmvm->Y[i], &ytDs);
237:               VecDotEnd(ldb->W, ldb->W, &stDs);

239:               ss_sum += PetscRealPart(stDs);
240:               ys_sum += PetscRealPart(ytDs);
241:               yy_sum += ldb->yty[i];
242:             }
243:           }
244:         } else if (1.0 == ldb->beta) {
245:           /*  Compute summations for scalar scaling */
246:           yy_sum = 0; /*  No safeguard required */
247:           ys_sum = 0; /*  No safeguard required */
248:           ss_sum = 0; /*  No safeguard required */
249:           start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
250:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
251:             VecPointwiseMult(ldb->V, lmvm->Y[i], ldb->invDnew);

253:             VecDotBegin(ldb->V, lmvm->S[i], &ytDs);
254:             VecDotBegin(ldb->V, ldb->V, &ytDy);
255:             VecDotEnd(ldb->V, lmvm->S[i], &ytDs);
256:             VecDotEnd(ldb->V, ldb->V, &ytDy);

258:             yy_sum += PetscRealPart(ytDy);
259:             ys_sum += PetscRealPart(ytDs);
260:             ss_sum += ldb->sts[i];
261:           }
262:         } else {
263:           VecCopy(ldb->invDnew, ldb->U);
264:           VecPow(ldb->U, ldb->beta-1);

266:           /*  Compute summations for scalar scaling */
267:           yy_sum = 0; /*  No safeguard required */
268:           ys_sum = 0; /*  No safeguard required */
269:           ss_sum = 0; /*  No safeguard required */
270:           start = PetscMax(0, lmvm->k-ldb->sigma_hist+1);
271:           for (i = start; i < PetscMin(lmvm->nupdates, ldb->sigma_hist); ++i) {
272:             VecPointwiseMult(ldb->V, ldb->invDnew, lmvm->Y[i]);
273:             VecPointwiseMult(ldb->W, ldb->U, lmvm->S[i]);

275:             VecDotBegin(ldb->V, ldb->W, &ytDs);
276:             VecDotBegin(ldb->V, ldb->V, &ytDy);
277:             VecDotBegin(ldb->W, ldb->W, &stDs);
278:             VecDotEnd(ldb->V, ldb->W, &ytDs);
279:             VecDotEnd(ldb->V, ldb->V, &ytDy);
280:             VecDotEnd(ldb->W, ldb->W, &stDs);

282:             yy_sum += PetscRealPart(ytDy);
283:             ys_sum += PetscRealPart(ytDs);
284:             ss_sum += PetscRealPart(stDs);
285:           }
286:         }

288:         if (0.0 == ldb->alpha) {
289:           /*  Safeguard ys_sum  */
290:           ys_sum = PetscMax(ldb->tol, ys_sum);

292:           sigma = ss_sum / ys_sum;
293:         } else if (1.0 == ldb->alpha) {
294:           /* yy_sum is never 0; if it were, we'd be at the minimum */
295:           sigma = ys_sum / yy_sum;
296:         } else {
297:           denom = 2.0*ldb->alpha*yy_sum;

299:           /*  Safeguard denom */
300:           denom = PetscMax(ldb->tol, denom);

302:           sigma = ((2.0*ldb->alpha-1)*ys_sum + PetscSqrtReal((2.0*ldb->alpha-1)*(2.0*ldb->alpha-1)*ys_sum*ys_sum - 4.0*ldb->alpha*(ldb->alpha-1)*yy_sum*ss_sum)) / denom;
303:         }
304:       } else {
305:         sigma = 1.0;
306:       }
307:       /*  If Q has small values, then Q^(r_beta - 1)
308:        can have very large values.  Hence, ys_sum
309:        and ss_sum can be infinity.  In this case,
310:        sigma can either be not-a-number or infinity. */

312:       if (PetscIsInfOrNanScalar(sigma)) {
313:         /*  sigma is not-a-number; skip rescaling */
314:       } else if (0.0 == sigma) {
315:         /*  sigma is zero; this is a bad case; skip rescaling */
316:       } else {
317:         /*  sigma is positive */
318:         VecScale(ldb->invDnew, sigma);
319:       }

321:       /* Combine the old diagonal and the new diagonal using a convex limiter */
322:       if (1.0 == ldb->rho) {
323:         VecCopy(ldb->invDnew, ldb->invD);
324:       } else if (ldb->rho) {
325:         VecAXPBY(ldb->invD, 1.0-ldb->rho, ldb->rho, ldb->invDnew);
326:       }
327:     } else {
328:       MatLMVMReset(B, PETSC_FALSE);
329:     }
330:     /* End DiagBrdn update */

332:   }
333:   /* Save the solution and function to be used in the next update */
334:   VecCopy(X, lmvm->Xprev);
335:   VecCopy(F, lmvm->Fprev);
336:   lmvm->prev_set = PETSC_TRUE;
337:   return(0);
338: }

340: /*------------------------------------------------------------*/

342: static PetscErrorCode MatCopy_DiagBrdn(Mat B, Mat M, MatStructure str)
343: {
344:   Mat_LMVM          *bdata = (Mat_LMVM*)B->data;
345:   Mat_DiagBrdn      *bctx = (Mat_DiagBrdn*)bdata->ctx;
346:   Mat_LMVM          *mdata = (Mat_LMVM*)M->data;
347:   Mat_DiagBrdn      *mctx = (Mat_DiagBrdn*)mdata->ctx;
348:   PetscErrorCode    ierr;
349:   PetscInt          i;

352:   mctx->theta = bctx->theta;
353:   mctx->alpha = bctx->alpha;
354:   mctx->beta = bctx->beta;
355:   mctx->rho = bctx->rho;
356:   mctx->delta = bctx->delta;
357:   mctx->delta_min = bctx->delta_min;
358:   mctx->delta_max = bctx->delta_max;
359:   mctx->tol = bctx->tol;
360:   mctx->sigma = bctx->sigma;
361:   mctx->sigma_hist = bctx->sigma_hist;
362:   mctx->forward = bctx->forward;
363:   VecCopy(bctx->invD, mctx->invD);
364:   for (i=0; i<=bdata->k; ++i) {
365:     mctx->yty[i] = bctx->yty[i];
366:     mctx->yts[i] = bctx->yts[i];
367:     mctx->sts[i] = mctx->sts[i];
368:   }
369:   return(0);
370: }

372: /*------------------------------------------------------------*/

374: static PetscErrorCode MatView_DiagBrdn(Mat B, PetscViewer pv)
375: {
376:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
377:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
378:   PetscErrorCode    ierr;
379:   PetscBool         isascii;

382:   PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
383:   if (isascii) {
384:     PetscViewerASCIIPrintf(pv,"Scale history: %d\n",ldb->sigma_hist);
385:     PetscViewerASCIIPrintf(pv,"Scale params: alpha=%g, beta=%g, rho=%g\n",(double)ldb->alpha, (double)ldb->beta, (double)ldb->rho);
386:     PetscViewerASCIIPrintf(pv,"Convex factor: theta=%g\n", (double)ldb->theta);
387:   }
388:   MatView_LMVM(B, pv);
389:   return(0);
390: }

392: /*------------------------------------------------------------*/

394: static PetscErrorCode MatSetFromOptions_DiagBrdn(PetscOptionItems *PetscOptionsObject, Mat B)
395: {
396:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
397:   Mat_DiagBrdn       *ldb = (Mat_DiagBrdn*)lmvm->ctx;
398:   PetscErrorCode    ierr;

401:   MatSetFromOptions_LMVM(PetscOptionsObject, B);
402:   PetscOptionsHead(PetscOptionsObject,"Restricted Broyden method for approximating SPD Jacobian actions (MATLMVMDIAGBRDN)");
403:   PetscOptionsReal("-mat_lmvm_theta","(developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling","",ldb->theta,&ldb->theta,NULL);
404:   PetscOptionsReal("-mat_lmvm_rho","(developer) update limiter in the J0 scaling","",ldb->rho,&ldb->rho,NULL);
405:   PetscOptionsReal("-mat_lmvm_tol","(developer) tolerance for bounding rescaling denominator","",ldb->tol,&ldb->tol,NULL);
406:   PetscOptionsReal("-mat_lmvm_alpha","(developer) convex ratio in the J0 scaling","",ldb->alpha,&ldb->alpha,NULL);
407:   PetscOptionsBool("-mat_lmvm_forward","Forward -> Update diagonal scaling for B. Else -> diagonal scaling for H.","",ldb->forward,&ldb->forward,NULL);
408:   PetscOptionsReal("-mat_lmvm_beta","(developer) exponential factor in the diagonal J0 scaling","",ldb->beta,&ldb->beta,NULL);
409:   PetscOptionsInt("-mat_lmvm_sigma_hist","(developer) number of past updates to use in the default J0 scalar","",ldb->sigma_hist,&ldb->sigma_hist,NULL);
410:   PetscOptionsTail();
411:   if ((ldb->theta < 0.0) || (ldb->theta > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio for the diagonal J0 scale cannot be outside the range of [0, 1]");
412:   if ((ldb->alpha < 0.0) || (ldb->alpha > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex ratio in the J0 scaling cannot be outside the range of [0, 1]");
413:   if ((ldb->rho < 0.0) || (ldb->rho > 1.0)) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "convex update limiter in the J0 scaling cannot be outside the range of [0, 1]");
414:   if (ldb->sigma_hist < 0) SETERRQ(PetscObjectComm((PetscObject)B), PETSC_ERR_ARG_OUTOFRANGE, "J0 scaling history length cannot be negative");
415:   return(0);
416: }

418: /*------------------------------------------------------------*/

420: static PetscErrorCode MatReset_DiagBrdn(Mat B, PetscBool destructive)
421: {
422:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
423:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
424:   PetscErrorCode    ierr;

427:   VecSet(ldb->invD, ldb->delta);
428:   if (destructive && ldb->allocated) {
429:     PetscFree3(ldb->yty, ldb->yts, ldb->sts);
430:     VecDestroy(&ldb->invDnew);
431:     VecDestroy(&ldb->invD);
432:     VecDestroy(&ldb->BFGS);
433:     VecDestroy(&ldb->DFP);
434:     VecDestroy(&ldb->U);
435:     VecDestroy(&ldb->V);
436:     VecDestroy(&ldb->W);
437:     ldb->allocated = PETSC_FALSE;
438:   }
439:   MatReset_LMVM(B, destructive);
440:   return(0);
441: }

443: /*------------------------------------------------------------*/

445: static PetscErrorCode MatAllocate_DiagBrdn(Mat B, Vec X, Vec F)
446: {
447:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
448:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
449:   PetscErrorCode    ierr;
450: 
452:   MatAllocate_LMVM(B, X, F);
453:   if (!ldb->allocated) {
454:     PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
455:     VecDuplicate(lmvm->Xprev, &ldb->invDnew);
456:     VecDuplicate(lmvm->Xprev, &ldb->invD);
457:     VecDuplicate(lmvm->Xprev, &ldb->BFGS);
458:     VecDuplicate(lmvm->Xprev, &ldb->DFP);
459:     VecDuplicate(lmvm->Xprev, &ldb->U);
460:     VecDuplicate(lmvm->Xprev, &ldb->V);
461:     VecDuplicate(lmvm->Xprev, &ldb->W);
462:     ldb->allocated = PETSC_TRUE;
463:   }
464:   return(0);
465: }

467: /*------------------------------------------------------------*/

469: static PetscErrorCode MatDestroy_DiagBrdn(Mat B)
470: {
471:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
472:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
473:   PetscErrorCode    ierr;

476:   if (ldb->allocated) {
477:     PetscFree3(ldb->yty, ldb->yts, ldb->sts);
478:     VecDestroy(&ldb->invDnew);
479:     VecDestroy(&ldb->invD);
480:     VecDestroy(&ldb->BFGS);
481:     VecDestroy(&ldb->DFP);
482:     VecDestroy(&ldb->U);
483:     VecDestroy(&ldb->V);
484:     VecDestroy(&ldb->W);
485:     ldb->allocated = PETSC_FALSE;
486:   }
487:   PetscFree(lmvm->ctx);
488:   MatDestroy_LMVM(B);
489:   return(0);
490: }

492: /*------------------------------------------------------------*/

494: static PetscErrorCode MatSetUp_DiagBrdn(Mat B)
495: {
496:   Mat_LMVM          *lmvm = (Mat_LMVM*)B->data;
497:   Mat_DiagBrdn      *ldb = (Mat_DiagBrdn*)lmvm->ctx;
498:   PetscErrorCode    ierr;

501:   MatSetUp_LMVM(B);
502:   if (!ldb->allocated) {
503:     PetscMalloc3(lmvm->m, &ldb->yty, lmvm->m, &ldb->yts, lmvm->m, &ldb->sts);
504:     VecDuplicate(lmvm->Xprev, &ldb->invDnew);
505:     VecDuplicate(lmvm->Xprev, &ldb->invD);
506:     VecDuplicate(lmvm->Xprev, &ldb->BFGS);
507:     VecDuplicate(lmvm->Xprev, &ldb->DFP);
508:     VecDuplicate(lmvm->Xprev, &ldb->U);
509:     VecDuplicate(lmvm->Xprev, &ldb->V);
510:     VecDuplicate(lmvm->Xprev, &ldb->W);
511:     ldb->allocated = PETSC_TRUE;
512:   }
513:   return(0);
514: }

516: /*------------------------------------------------------------*/

518: PetscErrorCode MatCreate_LMVMDiagBrdn(Mat B)
519: {
520:   Mat_LMVM          *lmvm;
521:   Mat_DiagBrdn      *ldb;
522:   PetscErrorCode    ierr;

525:   MatCreate_LMVM(B);
526:   PetscObjectChangeTypeName((PetscObject)B, MATLMVMDIAGBRDN);
527:   B->ops->setup = MatSetUp_DiagBrdn;
528:   B->ops->setfromoptions = MatSetFromOptions_DiagBrdn;
529:   B->ops->destroy = MatDestroy_DiagBrdn;
530:   B->ops->solve = MatSolve_DiagBrdn;
531:   B->ops->view = MatView_DiagBrdn;

533:   lmvm = (Mat_LMVM*)B->data;
534:   lmvm->square = PETSC_TRUE;
535:   lmvm->m = 1;
536:   lmvm->ops->allocate = MatAllocate_DiagBrdn;
537:   lmvm->ops->reset = MatReset_DiagBrdn;
538:   lmvm->ops->mult = MatMult_DiagBrdn;
539:   lmvm->ops->update = MatUpdate_DiagBrdn;
540:   lmvm->ops->copy = MatCopy_DiagBrdn;

542:   PetscNewLog(B, &ldb);
543:   lmvm->ctx = (void*)ldb;
544:   ldb->theta = 0.0;
545:   ldb->alpha = 1.0;
546:   ldb->rho = 1.0;
547:   ldb->forward = PETSC_TRUE;
548:   ldb->beta = 0.5;
549:   ldb->sigma = 1.0;
550:   ldb->delta = 1.0;
551:   ldb->delta_min = 1e-7;
552:   ldb->delta_max = 100.0;
553:   ldb->tol = 1e-8;
554:   ldb->sigma_hist = 1;
555:   ldb->allocated = PETSC_FALSE;
556:   return(0);
557: }

559: /*------------------------------------------------------------*/

561: /*@
562:    MatCreateLMVMDiagBrdn - DiagBrdn creates a symmetric Broyden-type diagonal matrix used 
563:    for approximating Hessians. It consists of a convex combination of DFP and BFGS 
564:    diagonal approximation schemes, such that DiagBrdn = (1-theta)*BFGS + theta*DFP. 
565:    To preserve symmetric positive-definiteness, we restrict theta to be in [0, 1]. 
566:    We also ensure positive definiteness by taking the VecAbs() of the final vector.

568:    There are two ways of approximating the diagonal: using the forward (B) update 
569:    schemes for BFGS and DFP and then taking the inverse, or directly working with 
570:    the inverse (H) update schemes for the BFGS and DFP updates, derived using the 
571:    Sherman-Morrison-Woodbury formula. We have implemented both, controlled by a 
572:    parameter below.

574:    In order to use the DiagBrdn matrix with other vector types, i.e. doing MatMults 
575:    and MatSolves, the matrix must first be created using MatCreate() and MatSetType(), 
576:    followed by MatLMVMAllocate(). Then it will be available for updating 
577:    (via MatLMVMUpdate) in one's favored solver implementation. 
578:    This allows for MPI compatibility.

580:    Collective

582:    Input Parameters:
583: +  comm - MPI communicator, set to PETSC_COMM_SELF
584: .  n - number of local rows for storage vectors
585: -  N - global size of the storage vectors

587:    Output Parameter:
588: .  B - the matrix

590:    It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions()
591:    paradigm instead of this routine directly.

593:    Options Database Keys:
594: +   -mat_lmvm_theta - (developer) convex ratio between BFGS and DFP components of the diagonal J0 scaling
595: .   -mat_lmvm_rho - (developer) update limiter for the J0 scaling
596: .   -mat_lmvm_alpha - (developer) coefficient factor for the quadratic subproblem in J0 scaling
597: .   -mat_lmvm_beta - (developer) exponential factor for the diagonal J0 scaling
598: .   -mat_lmvm_sigma_hist - (developer) number of past updates to use in J0 scaling.
599: .   -mat_lmvm_tol - (developer) tolerance for bounding the denominator of the rescaling away from 0.
600: -   -mat_lmvm_forward - (developer) whether or not to use the forward or backward Broyden update to the diagonal

602:    Level: intermediate

604: .seealso: MatCreate(), MATLMVM, MATLMVMDIAGBRDN, MatCreateLMVMDFP(), MatCreateLMVMSR1(),
605:           MatCreateLMVMBFGS(), MatCreateLMVMBrdn(), MatCreateLMVMSymBrdn()
606: @*/
607: PetscErrorCode MatCreateLMVMDiagBrdn(MPI_Comm comm, PetscInt n, PetscInt N, Mat *B)
608: {
609:   PetscErrorCode    ierr;

612:   MatCreate(comm, B);
613:   MatSetSizes(*B, n, n, N, N);
614:   MatSetType(*B, MATLMVMDIAGBRDN);
615:   MatSetUp(*B);
616:   return(0);
617: }