Actual source code: lmvmmat.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  1:  #include <../src/tao/matrix/lmvmmat.h>
  2:  #include <petsctao.h>
  3:  #include <petscksp.h>
  4:  #include <petsc/private/petscimpl.h>

  6: /* This is a vile hack */
  7: #if defined(PETSC_USE_COMPLEX)
  8: #define VecDot VecDotRealPart
  9: #endif

 11: #define TaoMid(a,b,c)    (((a) < (b)) ?                    \
 12:                            (((b) < (c)) ? (b) :            \
 13:                              (((a) < (c)) ? (c) : (a))) :  \
 14:                            (((a) < (c)) ? (a) :            \
 15:                              (((b) < (c)) ? (c) : (b))))

 17: /* These lists are used for setting options */
 18: static const char *Scale_Table[64] = {"none","scalar","broyden"};

 20: static const char *Rescale_Table[64] = {"none","scalar","gl"};

 22: static const char *Limit_Table[64] = {"none","average","relative","absolute"};

 24: /*@C
 25:   MatCreateLMVM - Creates a limited memory matrix for lmvm algorithms.

 27:   Collective on A

 29:   Input Parameters:
 30: + comm - MPI Communicator
 31: . n - local size of vectors
 32: - N - global size of vectors

 34:   Output Parameters:
 35: . A - New LMVM matrix

 37:   Level: developer

 39: @*/
 40: PetscErrorCode MatCreateLMVM(MPI_Comm comm, PetscInt n, PetscInt N, Mat *A)
 41: {
 42:   MatLMVMCtx     *ctx;
 44:   PetscInt       nhistory;

 47:   /*  create data structure and populate with default values */
 48:   PetscNew(&ctx);
 49:   ctx->lm=5;
 50:   ctx->eps=0.0;
 51:   ctx->limitType=MatLMVM_Limit_None;
 52:   ctx->scaleType=MatLMVM_Scale_Broyden;
 53:   ctx->rScaleType = MatLMVM_Rescale_Scalar;
 54:   ctx->s_alpha = 1.0;
 55:   ctx->r_alpha = 1.0;
 56:   ctx->r_beta = 0.5;
 57:   ctx->mu = 1.0;
 58:   ctx->nu = 100.0;

 60:   ctx->phi = 0.125;

 62:   ctx->scalar_history = 1;
 63:   ctx->rescale_history = 1;

 65:   ctx->delta_min = 1e-7;
 66:   ctx->delta_max = 100.0;
 67: 
 68:   ctx->recycle = PETSC_FALSE;

 70:   /*  Begin configuration */
 71:   PetscOptionsBegin(comm,NULL,NULL,NULL);
 72:   PetscOptionsInt("-tao_lmm_vectors", "vectors to use for approximation", "", ctx->lm, &ctx->lm,NULL);
 73:   PetscOptionsReal("-tao_lmm_limit_mu", "mu limiting factor", "", ctx->mu, &ctx->mu,NULL);
 74:   PetscOptionsReal("-tao_lmm_limit_nu", "nu limiting factor", "", ctx->nu, &ctx->nu,NULL);
 75:   PetscOptionsReal("-tao_lmm_broyden_phi", "phi factor for Broyden scaling", "", ctx->phi, &ctx->phi,NULL);
 76:   PetscOptionsReal("-tao_lmm_scalar_alpha", "alpha factor for scalar scaling", "",ctx->s_alpha, &ctx->s_alpha,NULL);
 77:   PetscOptionsReal("-tao_lmm_rescale_alpha", "alpha factor for rescaling diagonal", "", ctx->r_alpha, &ctx->r_alpha,NULL);
 78:   PetscOptionsReal("-tao_lmm_rescale_beta", "beta factor for rescaling diagonal", "", ctx->r_beta, &ctx->r_beta,NULL);
 79:   PetscOptionsInt("-tao_lmm_scalar_history", "amount of history for scalar scaling", "", ctx->scalar_history, &ctx->scalar_history,NULL);
 80:   PetscOptionsInt("-tao_lmm_rescale_history", "amount of history for rescaling diagonal", "", ctx->rescale_history, &ctx->rescale_history,NULL);
 81:   PetscOptionsReal("-tao_lmm_eps", "rejection tolerance", "", ctx->eps, &ctx->eps,NULL);
 82:   PetscOptionsEList("-tao_lmm_scale_type", "scale type", "", Scale_Table, MatLMVM_Scale_Types, Scale_Table[ctx->scaleType], &ctx->scaleType,NULL);
 83:   PetscOptionsEList("-tao_lmm_rescale_type", "rescale type", "", Rescale_Table, MatLMVM_Rescale_Types, Rescale_Table[ctx->rScaleType], &ctx->rScaleType,NULL);
 84:   PetscOptionsEList("-tao_lmm_limit_type", "limit type", "", Limit_Table, MatLMVM_Limit_Types, Limit_Table[ctx->limitType], &ctx->limitType,NULL);
 85:   PetscOptionsReal("-tao_lmm_delta_min", "minimum delta value", "", ctx->delta_min, &ctx->delta_min,NULL);
 86:   PetscOptionsReal("-tao_lmm_delta_max", "maximum delta value", "", ctx->delta_max, &ctx->delta_max,NULL);
 87:   PetscOptionsBool("-tao_lmm_recycle", "enable recycling of accumulated corrections between subsequent TaoSolve() calls", "", ctx->recycle, &ctx->recycle, NULL);
 88:   PetscOptionsEnd();

 90:   /*  Complete configuration */
 91:   ctx->rescale_history = PetscMin(ctx->rescale_history, ctx->lm);
 92:   PetscMalloc1(ctx->lm+1,&ctx->rho);
 93:   PetscMalloc1(ctx->lm+1,&ctx->beta);

 95:   nhistory = PetscMax(ctx->scalar_history,1);
 96:   PetscMalloc1(nhistory,&ctx->yy_history);
 97:   PetscMalloc1(nhistory,&ctx->ys_history);
 98:   PetscMalloc1(nhistory,&ctx->ss_history);

100:   nhistory = PetscMax(ctx->rescale_history,1);
101:   PetscMalloc1(nhistory,&ctx->yy_rhistory);
102:   PetscMalloc1(nhistory,&ctx->ys_rhistory);
103:   PetscMalloc1(nhistory,&ctx->ss_rhistory);

105:   /*  Finish initializations */
106:   ctx->lmnow = 0;
107:   ctx->iter = 0;
108:   ctx->nupdates = 0;
109:   ctx->nrejects = 0;
110:   ctx->delta = 1.0;

112:   ctx->Gprev = 0;
113:   ctx->Xprev = 0;

115:   ctx->scale = 0;
116:   ctx->useScale = PETSC_FALSE;

118:   ctx->H0_mat = 0;
119:   ctx->H0_ksp = 0;
120:   ctx->H0_norm = 0;
121:   ctx->useDefaultH0 = PETSC_TRUE;

123:   MatCreateShell(comm, n, n, N, N, ctx, A);
124:   MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
125:   MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
126:   return(0);
127: }

129: PetscErrorCode MatLMVMSolve(Mat A, Vec b, Vec x)
130: {
131:   PetscReal      sq, yq, dd;
132:   PetscInt       ll;
133:   PetscBool      scaled;
134:   MatLMVMCtx     *shell;

141:   MatShellGetContext(A,(void**)&shell);
142:   if (shell->lmnow < 1) {
143:     shell->rho[0] = 1.0;
144:     shell->theta = 1.0;
145:   }

147:   VecCopy(b,x);
148:   for (ll = 0; ll < shell->lmnow; ++ll) {
149:     VecDot(x,shell->S[ll],&sq);
150:     shell->beta[ll] = sq * shell->rho[ll];
151:     VecAXPY(x,-shell->beta[ll],shell->Y[ll]);
152:   }

154:   scaled = PETSC_FALSE;
155:   if (!scaled && !shell->useDefaultH0 && shell->H0_mat) {
156:     KSPSolve(shell->H0_ksp,x,shell->U);
157:     VecScale(shell->U, shell->theta);
158:     VecDot(x,shell->U,&dd);
159:     if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
160:       /*  Accept Hessian solve */
161:       VecCopy(shell->U,x);
162:       scaled = PETSC_TRUE;
163:     }
164:   }

166:   if (!scaled && shell->useScale) {
167:     VecPointwiseMult(shell->U,x,shell->scale);
168:     VecDot(x,shell->U,&dd);
169:     if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
170:       /*  Accept scaling */
171:       VecCopy(shell->U,x);
172:       scaled = PETSC_TRUE;
173:     }
174:   }

176:   if (!scaled) {
177:     switch(shell->scaleType) {
178:     case MatLMVM_Scale_None:
179:       break;

181:     case MatLMVM_Scale_Scalar:
182:       VecScale(x,shell->sigma);
183:       break;

185:     case MatLMVM_Scale_Broyden:
186:       VecPointwiseMult(x,x,shell->D);
187:       break;
188:     }
189:   }
190:   for (ll = shell->lmnow-1; ll >= 0; --ll) {
191:     VecDot(x,shell->Y[ll],&yq);
192:     VecAXPY(x,shell->beta[ll]-yq*shell->rho[ll],shell->S[ll]);
193:   }
194:   return(0);
195: }

197: PetscErrorCode MatView_LMVM(Mat A, PetscViewer pv)
198: {
199:   PetscBool      isascii;
201:   MatLMVMCtx     *lmP;

204:   MatShellGetContext(A,(void**)&lmP);
205:   PetscObjectTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii);
206:   if (isascii) {
207:     PetscViewerASCIIPrintf(pv,"LMVM Matrix\n");
208:     PetscViewerASCIIPrintf(pv," Number of vectors: %D\n",lmP->lm);
209:     PetscViewerASCIIPrintf(pv," scale type: %s\n",Scale_Table[lmP->scaleType]);
210:     PetscViewerASCIIPrintf(pv," rescale type: %s\n",Rescale_Table[lmP->rScaleType]);
211:     PetscViewerASCIIPrintf(pv," limit type: %s\n",Limit_Table[lmP->limitType]);
212:     PetscViewerASCIIPrintf(pv," updates: %D\n",lmP->nupdates);
213:     PetscViewerASCIIPrintf(pv," rejects: %D\n",lmP->nrejects);
214:   }
215:   return(0);
216: }

218: PetscErrorCode MatDestroy_LMVM(Mat M)
219: {
220:   MatLMVMCtx     *ctx;

224:   MatShellGetContext(M,(void**)&ctx);
225:   if (ctx->allocated) {
226:     if (ctx->Xprev) {
227:       PetscObjectDereference((PetscObject)ctx->Xprev);
228:     }
229:     if (ctx->Gprev) {
230:       PetscObjectDereference((PetscObject)ctx->Gprev);
231:     }

233:     VecDestroyVecs(ctx->lm+1,&ctx->S);
234:     VecDestroyVecs(ctx->lm+1,&ctx->Y);
235:     VecDestroy(&ctx->D);
236:     VecDestroy(&ctx->U);
237:     VecDestroy(&ctx->V);
238:     VecDestroy(&ctx->W);
239:     VecDestroy(&ctx->P);
240:     VecDestroy(&ctx->Q);
241:     VecDestroy(&ctx->H0_norm);
242:     if (ctx->scale) {
243:       VecDestroy(&ctx->scale);
244:     }
245:   }
246:   if (ctx->H0_mat) {
247:     PetscObjectDereference((PetscObject)ctx->H0_mat);
248:     KSPDestroy(&ctx->H0_ksp);
249:   }
250:   PetscFree(ctx->rho);
251:   PetscFree(ctx->beta);
252:   PetscFree(ctx->yy_history);
253:   PetscFree(ctx->ys_history);
254:   PetscFree(ctx->ss_history);
255:   PetscFree(ctx->yy_rhistory);
256:   PetscFree(ctx->ys_rhistory);
257:   PetscFree(ctx->ss_rhistory);
258:   PetscFree(ctx);
259:   return(0);
260: }

262: PetscErrorCode MatLMVMGetUpdates(Mat M, PetscInt *nupdates)
263: {
265:   MatLMVMCtx     *ctx;
266: 
268:   MatShellGetContext(M, (void**)&ctx);
269:   *nupdates = ctx->nupdates;
270:   return(0);
271: }

273: PetscErrorCode MatLMVMSetRecycleFlag(Mat M, PetscBool flg)
274: {
276:   MatLMVMCtx     *ctx;
277: 
279:   MatShellGetContext(M, (void**)&ctx);
280:   ctx->recycle = flg;
281:   return(0);
282: }

284: PetscErrorCode MatLMVMGetRecycleFlag(Mat M, PetscBool *flg)
285: {
287:   MatLMVMCtx     *ctx;
288: 
290:   MatShellGetContext(M, (void**)&ctx);
291:   *flg = ctx->recycle;
292:   return(0);
293: }

295: PetscErrorCode MatLMVMReset(Mat M)
296: {
298:   MatLMVMCtx     *ctx;
299:   PetscInt       i;

302:   MatShellGetContext(M,(void**)&ctx);
303:   if (ctx->recycle) {
304:     PetscInfo(M, "WARNING: Correction vectors are being reset while recycling is enabled!\n");
305:   }
306:   if (ctx->Gprev) {
307:     PetscObjectDereference((PetscObject)ctx->Gprev);
308:   }
309:   if (ctx->Xprev) {
310:     PetscObjectDereference((PetscObject)ctx->Xprev);
311:   }
312:   ctx->Gprev = ctx->Y[ctx->lm];
313:   ctx->Xprev = ctx->S[ctx->lm];
314:   PetscObjectReference((PetscObject)ctx->Gprev);
315:   PetscObjectReference((PetscObject)ctx->Xprev);
316:   for (i=0; i<ctx->lm; ++i) {
317:     ctx->rho[i] = 0.0;
318:   }
319:   ctx->rho[0] = 1.0;

321:   /*  Set the scaling and diagonal scaling matrix */
322:   switch(ctx->scaleType) {
323:   case MatLMVM_Scale_None:
324:     ctx->sigma = 1.0;
325:     break;
326:   case MatLMVM_Scale_Scalar:
327:     ctx->sigma = ctx->delta;
328:     break;
329:   case MatLMVM_Scale_Broyden:
330:     VecSet(ctx->D,ctx->delta);
331:     break;
332:   }

334:   ctx->iter=0;
335:   ctx->nupdates=0;
336:   ctx->lmnow=0;
337:   return(0);
338: }

340: PetscErrorCode MatLMVMUpdate(Mat M, Vec x, Vec g)
341: {
342:   MatLMVMCtx     *ctx;
343:   PetscReal      rhotemp, rhotol;
344:   PetscReal      y0temp, s0temp;
345:   PetscReal      yDy, yDs, sDs;
346:   PetscReal      sigmanew, denom;
348:   PetscInt       i;
349:   PetscBool      same;
350:   PetscReal      yy_sum=0.0, ys_sum=0.0, ss_sum=0.0;

355:   PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
356:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
357:   MatShellGetContext(M,(void**)&ctx);
358:   if (!ctx->allocated) {
359:     MatLMVMAllocateVectors(M, x);
360:   }

362:   if (0 == ctx->iter) {
363:     MatLMVMReset(M);
364:   }  else {
365:     VecAYPX(ctx->Gprev,-1.0,g);
366:     VecAYPX(ctx->Xprev,-1.0,x);

368:     VecDot(ctx->Gprev,ctx->Xprev,&rhotemp);
369:     VecDot(ctx->Gprev,ctx->Gprev,&y0temp);

371:     rhotol = ctx->eps * y0temp;
372:     if (rhotemp > rhotol) {
373:       ++ctx->nupdates;

375:       ctx->lmnow = PetscMin(ctx->lmnow+1, ctx->lm);
376:       ierr=PetscObjectDereference((PetscObject)ctx->S[ctx->lm]);
377:       ierr=PetscObjectDereference((PetscObject)ctx->Y[ctx->lm]);
378:       for (i = ctx->lm-1; i >= 0; --i) {
379:         ctx->S[i+1] = ctx->S[i];
380:         ctx->Y[i+1] = ctx->Y[i];
381:         ctx->rho[i+1] = ctx->rho[i];
382:       }
383:       ctx->S[0] = ctx->Xprev;
384:       ctx->Y[0] = ctx->Gprev;
385:       PetscObjectReference((PetscObject)ctx->S[0]);
386:       PetscObjectReference((PetscObject)ctx->Y[0]);
387:       ctx->rho[0] = 1.0 / rhotemp;

389:       /*  Compute the scaling */
390:       switch(ctx->scaleType) {
391:       case MatLMVM_Scale_None:
392:         break;

394:       case MatLMVM_Scale_Scalar:
395:         /*  Compute s^T s  */
396:           VecDot(ctx->Xprev,ctx->Xprev,&s0temp);

398:         /*  Scalar is positive; safeguards are not required. */

400:         /*  Save information for scalar scaling */
401:         ctx->yy_history[(ctx->nupdates - 1) % ctx->scalar_history] = y0temp;
402:         ctx->ys_history[(ctx->nupdates - 1) % ctx->scalar_history] = rhotemp;
403:         ctx->ss_history[(ctx->nupdates - 1) % ctx->scalar_history] = s0temp;

405:         /*  Compute summations for scalar scaling */
406:         yy_sum = 0;     /*  No safeguard required; y^T y > 0 */
407:         ys_sum = 0;     /*  No safeguard required; y^T s > 0 */
408:         ss_sum = 0;     /*  No safeguard required; s^T s > 0 */
409:         for (i = 0; i < PetscMin(ctx->nupdates, ctx->scalar_history); ++i) {
410:           yy_sum += ctx->yy_history[i];
411:           ys_sum += ctx->ys_history[i];
412:           ss_sum += ctx->ss_history[i];
413:         }

415:         if (0.0 == ctx->s_alpha) {
416:           /*  Safeguard ys_sum  */
417:           if (0.0 == ys_sum) {
418:             ys_sum = TAO_ZERO_SAFEGUARD;
419:           }

421:           sigmanew = ss_sum / ys_sum;
422:         } else if (1.0 == ctx->s_alpha) {
423:           /*  Safeguard yy_sum  */
424:           if (0.0 == yy_sum) {
425:             yy_sum = TAO_ZERO_SAFEGUARD;
426:           }

428:           sigmanew = ys_sum / yy_sum;
429:         } else {
430:           denom = 2*ctx->s_alpha*yy_sum;

432:           /*  Safeguard denom */
433:           if (0.0 == denom) {
434:             denom = TAO_ZERO_SAFEGUARD;
435:           }

437:           sigmanew = ((2*ctx->s_alpha-1)*ys_sum +  PetscSqrtReal((2*ctx->s_alpha-1)*(2*ctx->s_alpha-1)*ys_sum*ys_sum - 4*(ctx->s_alpha)*(ctx->s_alpha-1)*yy_sum*ss_sum)) / denom;
438:         }

440:         switch(ctx->limitType) {
441:         case MatLMVM_Limit_Average:
442:           if (1.0 == ctx->mu) {
443:             ctx->sigma = sigmanew;
444:           } else if (ctx->mu) {
445:             ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
446:           }
447:           break;

449:         case MatLMVM_Limit_Relative:
450:           if (ctx->mu) {
451:             ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
452:           }
453:           break;

455:         case MatLMVM_Limit_Absolute:
456:           if (ctx->nu) {
457:             ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
458:           }
459:           break;

461:         default:
462:           ctx->sigma = sigmanew;
463:           break;
464:         }
465:         break;

467:       case MatLMVM_Scale_Broyden:
468:         /*  Original version */
469:         /*  Combine DFP and BFGS */

471:         /*  This code appears to be numerically unstable.  We use the */
472:         /*  original version because this was used to generate all of */
473:         /*  the data and because it may be the least unstable of the */
474:         /*  bunch. */

476:         /*  P = Q = inv(D); */
477:         VecCopy(ctx->D,ctx->P);
478:         VecReciprocal(ctx->P);
479:         VecCopy(ctx->P,ctx->Q);

481:         /*  V = y*y */
482:         VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev);

484:         /*  W = inv(D)*s */
485:         VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P);
486:         VecDot(ctx->W,ctx->Xprev,&sDs);

488:         /*  Safeguard rhotemp and sDs */
489:         if (0.0 == rhotemp) {
490:           rhotemp = TAO_ZERO_SAFEGUARD;
491:         }

493:         if (0.0 == sDs) {
494:           sDs = TAO_ZERO_SAFEGUARD;
495:         }

497:         if (1.0 != ctx->phi) {
498:           /*  BFGS portion of the update */
499:           /*  U = (inv(D)*s)*(inv(D)*s) */
500:           VecPointwiseMult(ctx->U,ctx->W,ctx->W);

502:           /*  Assemble */
503:           VecAXPY(ctx->P,1.0/rhotemp,ctx->V);
504:           VecAXPY(ctx->P,-1.0/sDs,ctx->U);
505:         }

507:         if (0.0 != ctx->phi) {
508:           /*  DFP portion of the update */
509:           /*  U = inv(D)*s*y */
510:           VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev);

512:           /*  Assemble */
513:           VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V);
514:           VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U);
515:         }

517:         if (0.0 == ctx->phi) {
518:             VecCopy(ctx->P,ctx->U);
519:         } else if (1.0 == ctx->phi) {
520:             VecCopy(ctx->Q,ctx->U);
521:         } else {
522:           /*  Broyden update U=(1-phi)*P + phi*Q */
523:             VecCopy(ctx->Q,ctx->U);
524:             VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P);
525:         }

527:         /*  Obtain inverse and ensure positive definite */
528:         VecReciprocal(ctx->U);
529:         VecAbs(ctx->U);

531:         switch(ctx->rScaleType) {
532:         case MatLMVM_Rescale_None:
533:             break;

535:         case MatLMVM_Rescale_Scalar:
536:         case MatLMVM_Rescale_GL:
537:           if (ctx->rScaleType == MatLMVM_Rescale_GL) {
538:             /*  Gilbert and Lemarachal use the old diagonal */
539:             VecCopy(ctx->D,ctx->P);
540:           } else {
541:             /*  The default version uses the current diagonal */
542:               VecCopy(ctx->U,ctx->P);
543:           }

545:           /*  Compute s^T s  */
546:           VecDot(ctx->Xprev,ctx->Xprev,&s0temp);

548:           /*  Save information for special cases of scalar rescaling */
549:           ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
550:           ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
551:           ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;

553:           if (0.5 == ctx->r_beta) {
554:             if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
555:               VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P);
556:               VecDot(ctx->V,ctx->Y[0],&yy_sum);

558:               VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);
559:               VecDot(ctx->W,ctx->S[0],&ss_sum);

561:               ys_sum = ctx->ys_rhistory[0];
562:             } else {
563:               VecCopy(ctx->P,ctx->Q);
564:               VecReciprocal(ctx->Q);

566:               /*  Compute summations for scalar scaling */
567:               yy_sum = 0;       /*  No safeguard required */
568:               ys_sum = 0;       /*  No safeguard required */
569:               ss_sum = 0;       /*  No safeguard required */
570:               for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
571:                 VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P);
572:                 VecDot(ctx->V,ctx->Y[i],&yDy);
573:                 yy_sum += yDy;

575:                 VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q);
576:                 VecDot(ctx->W,ctx->S[i],&sDs);
577:                 ss_sum += sDs;
578:                 ys_sum += ctx->ys_rhistory[i];
579:               }
580:             }
581:           } else if (0.0 == ctx->r_beta) {
582:             if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
583:               /*  Compute summations for scalar scaling */
584:               VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);

586:               VecDot(ctx->W, ctx->Y[0], &ys_sum);
587:               VecDot(ctx->W, ctx->W, &ss_sum);
588:               yy_sum += ctx->yy_rhistory[0];
589:             } else {
590:               VecCopy(ctx->Q, ctx->P);
591:               VecReciprocal(ctx->Q);

593:               /*  Compute summations for scalar scaling */
594:               yy_sum = 0;       /*  No safeguard required */
595:               ys_sum = 0;       /*  No safeguard required */
596:               ss_sum = 0;       /*  No safeguard required */
597:               for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
598:                 VecPointwiseMult(ctx->W, ctx->S[i], ctx->Q);
599:                 VecDot(ctx->W, ctx->Y[i], &yDs);
600:                 ys_sum += yDs;

602:                 VecDot(ctx->W, ctx->W, &sDs);
603:                 ss_sum += sDs;

605:                 yy_sum += ctx->yy_rhistory[i];
606:               }
607:             }
608:           } else if (1.0 == ctx->r_beta) {
609:             /*  Compute summations for scalar scaling */
610:             yy_sum = 0; /*  No safeguard required */
611:             ys_sum = 0; /*  No safeguard required */
612:             ss_sum = 0; /*  No safeguard required */
613:             for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
614:               VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P);
615:               VecDot(ctx->V, ctx->S[i], &yDs);
616:               ys_sum += yDs;

618:               VecDot(ctx->V, ctx->V, &yDy);
619:               yy_sum += yDy;

621:               ss_sum += ctx->ss_rhistory[i];
622:             }
623:           } else {
624:             VecCopy(ctx->Q, ctx->P);

626:             VecPow(ctx->P, ctx->r_beta);
627:             VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q);

629:             /*  Compute summations for scalar scaling */
630:             yy_sum = 0; /*  No safeguard required */
631:             ys_sum = 0; /*  No safeguard required */
632:             ss_sum = 0; /*  No safeguard required */
633:             for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
634:               VecPointwiseMult(ctx->V, ctx->P, ctx->Y[i]);
635:               VecPointwiseMult(ctx->W, ctx->Q, ctx->S[i]);

637:               VecDot(ctx->V, ctx->V, &yDy);
638:               VecDot(ctx->V, ctx->W, &yDs);
639:               VecDot(ctx->W, ctx->W, &sDs);

641:               yy_sum += yDy;
642:               ys_sum += yDs;
643:               ss_sum += sDs;
644:             }
645:           }

647:           if (0.0 == ctx->r_alpha) {
648:             /*  Safeguard ys_sum  */
649:             if (0.0 == ys_sum) {
650:               ys_sum = TAO_ZERO_SAFEGUARD;
651:             }

653:             sigmanew = ss_sum / ys_sum;
654:           } else if (1.0 == ctx->r_alpha) {
655:             /*  Safeguard yy_sum  */
656:             if (0.0 == yy_sum) {
657:               ys_sum = TAO_ZERO_SAFEGUARD;
658:             }

660:             sigmanew = ys_sum / yy_sum;
661:           } else {
662:             denom = 2*ctx->r_alpha*yy_sum;

664:             /*  Safeguard denom */
665:             if (0.0 == denom) {
666:               denom = TAO_ZERO_SAFEGUARD;
667:             }

669:             sigmanew = ((2*ctx->r_alpha-1)*ys_sum + PetscSqrtReal((2*ctx->r_alpha-1)*(2*ctx->r_alpha-1)*ys_sum*ys_sum - 4*ctx->r_alpha*(ctx->r_alpha-1)*yy_sum*ss_sum)) / denom;
670:           }

672:           /*  If Q has small values, then Q^(r_beta - 1) */
673:           /*  can have very large values.  Hence, ys_sum */
674:           /*  and ss_sum can be infinity.  In this case, */
675:           /*  sigmanew can either be not-a-number or infinity. */

677:           if (PetscIsInfOrNanReal(sigmanew)) {
678:             /*  sigmanew is not-a-number; skip rescaling */
679:           } else if (!sigmanew) {
680:             /*  sigmanew is zero; this is a bad case; skip rescaling */
681:           } else {
682:             /*  sigmanew is positive */
683:             VecScale(ctx->U, sigmanew);
684:           }
685:           break;
686:         }

688:         /*  Modify for previous information */
689:         switch(ctx->limitType) {
690:         case MatLMVM_Limit_Average:
691:           if (1.0 == ctx->mu) {
692:             VecCopy(ctx->D, ctx->U);
693:           } else if (ctx->mu) {
694:             VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U);
695:           }
696:           break;

698:         case MatLMVM_Limit_Relative:
699:           if (ctx->mu) {
700:             /*  P = (1-mu) * D */
701:             VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D);
702:             /*  Q = (1+mu) * D */
703:             VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D);
704:             VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D);
705:           }
706:           break;

708:         case MatLMVM_Limit_Absolute:
709:           if (ctx->nu) {
710:             VecCopy(ctx->P, ctx->D);
711:             VecShift(ctx->P, -ctx->nu);
712:             VecCopy(ctx->D, ctx->Q);
713:             VecShift(ctx->Q, ctx->nu);
714:             VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P);
715:           }
716:           break;

718:         default:
719:             VecCopy(ctx->U, ctx->D);
720:           break;
721:         }
722:         break;
723:       }
724:       PetscObjectDereference((PetscObject)ctx->Xprev);
725:       PetscObjectDereference((PetscObject)ctx->Gprev);
726:       ctx->Xprev = ctx->S[ctx->lm];
727:       ctx->Gprev = ctx->Y[ctx->lm];
728:       PetscObjectReference((PetscObject)ctx->S[ctx->lm]);
729:       PetscObjectReference((PetscObject)ctx->Y[ctx->lm]);

731:     } else {
732:       ++ctx->nrejects;
733:     }
734:   }

736:   ++ctx->iter;
737:   VecCopy(x, ctx->Xprev);
738:   VecCopy(g, ctx->Gprev);
739:   return(0);
740: }

742: PetscErrorCode MatLMVMSetDelta(Mat m, PetscReal d)
743: {
744:   MatLMVMCtx     *ctx;
746:   PetscBool      same;

750:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
751:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
752:   MatShellGetContext(m,(void**)&ctx);
753:   ctx->delta = PetscAbsReal(d);
754:   ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
755:   ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
756:   return(0);
757: }

759: PetscErrorCode MatLMVMSetScale(Mat m, Vec s)
760: {
761:   MatLMVMCtx     *ctx;
763:   PetscBool      same;

767:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
768:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
769:   MatShellGetContext(m,(void**)&ctx);

771:   VecDestroy(&ctx->scale);
772:   if (s) {
773:     VecDuplicate(s,&ctx->scale);
774:   }
775:   return(0);
776: }

778: PetscErrorCode MatLMVMGetRejects(Mat m, PetscInt *nrejects)
779: {
780:   MatLMVMCtx     *ctx;
782:   PetscBool      same;

786:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
787:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
788:   MatShellGetContext(m,(void**)&ctx);
789:   *nrejects = ctx->nrejects;
790:   return(0);
791: }

793: PetscErrorCode MatLMVMSetH0(Mat m, Mat H0)
794: {
795:   MatLMVMCtx     *ctx;
797:   PetscBool      same;

801:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
802:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
803:   MatShellGetContext(m,(void**)&ctx);

805:   ctx->H0_mat = H0;
806:   PetscObjectReference((PetscObject)ctx->H0_mat);

808:   ctx->useDefaultH0 = PETSC_FALSE;

810:   KSPCreate(PetscObjectComm((PetscObject)H0), &ctx->H0_ksp);
811:   PetscObjectIncrementTabLevel((PetscObject)ctx->H0_ksp, (PetscObject)ctx, 1);
812:   KSPSetOperators(ctx->H0_ksp, H0, H0);
813:   /* its options prefix and setup is handled in TaoSolve_LMVM/TaoSolve_BLMVM */
814:   return(0);
815: }

817: PetscErrorCode MatLMVMGetH0(Mat m, Mat *H0)
818: {
819:   MatLMVMCtx     *ctx;
821:   PetscBool      same;

825:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
826:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");

828:   MatShellGetContext(m,(void**)&ctx);
829:   *H0  = ctx->H0_mat;
830:   return(0);
831: }

833: PetscErrorCode MatLMVMGetH0KSP(Mat m, KSP *H0ksp)
834: {
835:   MatLMVMCtx     *ctx;
837:   PetscBool      same;

841:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
842:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");

844:   MatShellGetContext(m,(void**)&ctx);
845:   *H0ksp  = ctx->H0_ksp;
846:   return(0);
847: }

849: PetscErrorCode MatLMVMGetX0(Mat m, Vec x)
850: {
852:     return(0);
853: }

855: PetscErrorCode MatLMVMSetPrev(Mat M, Vec x, Vec g)
856: {
857:   MatLMVMCtx     *ctx;
859:   PetscBool      same;

864:   PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
865:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
866:   MatShellGetContext(M,(void**)&ctx);
867:   if (ctx->nupdates == 0) {
868:     MatLMVMUpdate(M,x,g);
869:   } else {
870:     VecCopy(x,ctx->Xprev);
871:     VecCopy(g,ctx->Gprev);
872:     /*  TODO scaling specific terms */
873:   }
874:   return(0);
875: }

877: PetscErrorCode MatLMVMRefine(Mat coarse, Mat op, Vec fineX, Vec fineG)
878: {
880:   PetscBool      same;

887:   PetscObjectTypeCompare((PetscObject)coarse,MATSHELL,&same);
888:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
889:   PetscObjectTypeCompare((PetscObject)op,MATSHELL,&same);
890:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
891:   return(0);
892: }

894: PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v)
895: {
897:   MatLMVMCtx     *ctx;
898:   PetscBool      same;

903:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
904:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
905:   MatShellGetContext(m,(void**)&ctx);

907:   /*  Perform allocations */
908:   VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
909:   VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
910:   VecDuplicate(v,&ctx->D);
911:   VecDuplicate(v,&ctx->U);
912:   VecDuplicate(v,&ctx->V);
913:   VecDuplicate(v,&ctx->W);
914:   VecDuplicate(v,&ctx->P);
915:   VecDuplicate(v,&ctx->Q);
916:   VecDuplicate(v,&ctx->H0_norm);
917:   ctx->allocated = PETSC_TRUE;
918:   return(0);
919: }