Actual source code: lmvmmat.c

  1: #include "lmvmmat.h"   /*I "lmvmmat.h" */
  2: #include "tao_util.h"  /*I "tao_util.h" */
  3: #include "petscksp.h"

  5: #define TaoMid(a,b,c)    (((a) < (b)) ?                    \
  6:                            (((b) < (c)) ? (b) :            \
  7:                              (((a) < (c)) ? (c) : (a))) :  \
  8:                            (((a) < (c)) ? (a) :            \
  9:                              (((b) < (c)) ? (c) : (b))))

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

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

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


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

 30:   Collective on A

 32:   Input Parameters:
 33: + comm - MPI Communicator
 34: . n - local size of vectors
 35: - N - global size of vectors

 37:   Output Parameters:
 38: . A - New LMVM matrix

 40:   Level: developer

 42: @*/
 44: {
 45:     MatLMVMCtx *ctx;
 47:     PetscInt nhistory;


 51:     /*  create data structure and populate with default values */
 52:     PetscMalloc(sizeof(MatLMVMCtx),(void**)&ctx); 
 53:     ctx->lm=5;
 54:     ctx->eps=0.0;
 55:     ctx->limitType=MatLMVM_Limit_None;
 56:     ctx->scaleType=MatLMVM_Scale_Broyden;
 57:     ctx->rScaleType = MatLMVM_Rescale_Scalar;
 58:     ctx->s_alpha = 1.0;
 59:     ctx->r_alpha = 1.0;
 60:     ctx->r_beta = 0.5;
 61:     ctx->mu = 1.0;
 62:     ctx->nu = 100.0;

 64:     ctx->phi = 0.125;                

 66:     ctx->scalar_history = 1;
 67:     ctx->rescale_history = 1;

 69:     ctx->delta_min = 1e-7;
 70:     ctx->delta_max = 100.0;

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

 89:     /*  Complete configuration */
 90:     ctx->rescale_history = PetscMin(ctx->rescale_history, ctx->lm);


 93:     PetscMalloc((ctx->lm+1)*sizeof(PetscReal),(void**)&ctx->rho); 
 94:                        
 95:     PetscMalloc((ctx->lm+1)*sizeof(PetscReal),(void**)&ctx->beta); 
 96:                        

 98:     nhistory = PetscMax(ctx->scalar_history,1);
 99:     PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->yy_history); 
100:                        
101:     PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->ys_history);
102:                        
103:     PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->ss_history); 
104:                        

106:     nhistory = PetscMax(ctx->rescale_history,1);
107:     PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->yy_rhistory);
108:                        
109:     PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->ys_rhistory);
110:                        
111:     PetscMalloc(nhistory*sizeof(PetscReal),(void**)&ctx->ss_rhistory); 
112:                        


115:     /*  Finish initializations */
116:     ctx->lmnow = 0;
117:     ctx->iter = 0;
118:     ctx->nupdates = 0;
119:     ctx->nrejects = 0;
120:     ctx->delta = 1.0;

122:     ctx->Gprev = 0;
123:     ctx->Xprev = 0;

125:     ctx->scale = 0;
126:     ctx->useScale = PETSC_FALSE;

128:     ctx->H0 = 0;
129:     ctx->useDefaultH0=PETSC_TRUE;
130:     
131:     MatCreateShell(comm, n, n, N, N, ctx, A); 
132:     MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
133:     
134:     MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
135:     

137:     return(0);
138: }



142:   
143:   
147: {
148:     PetscReal      sq, yq, dd;
149:     PetscInt       ll;
150:     PetscBool     scaled;
151:     MatLMVMCtx     *shell;

158:     MatShellGetContext(A,(void**)&shell); 
159:     if (shell->lmnow < 1) {
160:         shell->rho[0] = 1.0;
161:     }

163:     VecCopy(b,x); 
164:     for (ll = 0; ll < shell->lmnow; ++ll) {
165:         VecDot(x,shell->S[ll],&sq); 
166:         shell->beta[ll] = sq * shell->rho[ll];
167:         VecAXPY(x,-shell->beta[ll],shell->Y[ll]); 
168:     }

170:     scaled = PETSC_FALSE;
171:     if (!scaled && !shell->useDefaultH0 && shell->H0) {
172:         MatSolve(shell->H0,x,shell->U); 
173:         VecDot(x,shell->U,&dd); 
174:         if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
175:             /*  Accept Hessian solve */
176:             VecCopy(shell->U,x); 
177:             scaled = PETSC_TRUE;
178:         }
179:     }

181:     if (!scaled && shell->useScale) {
182:         VecPointwiseMult(shell->U,x,shell->scale); 
183:         VecDot(x,shell->U,&dd); 
184:         if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
185:             /*  Accept scaling */
186:             VecCopy(shell->U,x);         
187:             scaled = PETSC_TRUE;
188:         }
189:     }
190:   
191:     if (!scaled) {
192:         switch(shell->scaleType) {
193:             case MatLMVM_Scale_None:
194:                 break;

196:             case MatLMVM_Scale_Scalar:
197:                 VecScale(x,shell->sigma); 
198:                 break;
199:   
200:             case MatLMVM_Scale_Broyden:
201:                 VecPointwiseMult(x,x,shell->D); 
202:                 break;
203:         }
204:     } 

206:     for (ll = shell->lmnow-1; ll >= 0; --ll) {
207:         VecDot(x,shell->Y[ll],&yq); 
208:         VecAXPY(x,shell->beta[ll]-yq*shell->rho[ll],shell->S[ll]);
209:         
210:     }
211:     return(0);
212: }



216:   
220: {
221:     PetscBool isascii;
223:     MatLMVMCtx *lmP;

226:     MatShellGetContext(A,(void**)&lmP); 

228:     PetscTypeCompare((PetscObject)pv,PETSCVIEWERASCII,&isascii); 
229:     if (isascii) {
230:         PetscViewerASCIIPrintf(pv,"LMVM Matrix\n"); 
231:         PetscViewerASCIIPrintf(pv," Number of vectors: %D\n",lmP->lm); 
232:         PetscViewerASCIIPrintf(pv," scale type: %s\n",Scale_Table[lmP->scaleType]); 
233:         PetscViewerASCIIPrintf(pv," rescale type: %s\n",Rescale_Table[lmP->rScaleType]); 
234:         PetscViewerASCIIPrintf(pv," limit type: %s\n",Limit_Table[lmP->limitType]); 
235:         PetscViewerASCIIPrintf(pv," updates: %D\n",lmP->nupdates); 
236:         PetscViewerASCIIPrintf(pv," rejects: %D\n",lmP->nrejects); 
237:         
238:     }
239:     else {
240:         SETERRQ(PETSC_COMM_SELF,1,"non-ascii viewers not implemented for MatLMVM\n");
241:     }

243:     return(0);
244:     
245: }

250: {
251:     MatLMVMCtx     *ctx;

255:     MatShellGetContext(M,(void**)&ctx); 
256:     if (ctx->allocated) {
257:       if (ctx->Xprev) {
258:         PetscObjectDereference((PetscObject)ctx->Xprev); 
259:       }
260:       if (ctx->Gprev) {
261:         PetscObjectDereference((PetscObject)ctx->Gprev); 
262:       }

264:       VecDestroyVecs(ctx->lm+1,&ctx->S); 
265:       VecDestroyVecs(ctx->lm+1,&ctx->Y); 
266:       VecDestroy(&ctx->D); 
267:       VecDestroy(&ctx->U); 
268:       VecDestroy(&ctx->V); 
269:       VecDestroy(&ctx->W); 
270:       VecDestroy(&ctx->P); 
271:       VecDestroy(&ctx->Q); 
272:       if (ctx->scale) {
273:         VecDestroy(&ctx->scale); 
274:       }
275:     }
276:     PetscFree(ctx->rho); 
277:     PetscFree(ctx->beta); 
278:     PetscFree(ctx->yy_history); 
279:     PetscFree(ctx->ys_history); 
280:     PetscFree(ctx->ss_history); 
281:     PetscFree(ctx->yy_rhistory); 
282:     PetscFree(ctx->ys_rhistory); 
283:     PetscFree(ctx->ss_rhistory); 

285:     PetscFree(ctx); 

287:     return(0);
288:         
289: }


295: {
297:     MatLMVMCtx *ctx;
298:     PetscInt i;
300:     MatShellGetContext(M,(void**)&ctx); 
301:     if (ctx->Gprev) {
302:       PetscObjectDereference((PetscObject)ctx->Gprev); 
303:     }
304:     if (ctx->Xprev) {
305:       PetscObjectDereference((PetscObject)ctx->Xprev); 
306:     }
307:     ctx->Gprev = ctx->Y[ctx->lm];
308:     ctx->Xprev = ctx->S[ctx->lm];
309:     PetscObjectReference((PetscObject)ctx->Gprev); 
310:     PetscObjectReference((PetscObject)ctx->Xprev); 
311:     for (i=0; i<ctx->lm; ++i) {
312:       ctx->rho[i] = 0.0;
313:     }
314:     ctx->rho[0] = 1.0;
315:     
316:     /*  Set the scaling and diagonal scaling matrix */
317:     switch(ctx->scaleType) {
318:       case MatLMVM_Scale_None:
319:         ctx->sigma = 1.0;
320:         break;
321:       case MatLMVM_Scale_Scalar:
322:         ctx->sigma = ctx->delta;
323:         break;
324:       case MatLMVM_Scale_Broyden:
325:         VecSet(ctx->D,ctx->delta); 
326:         break;
327:     }

329:     ctx->iter=0;
330:     ctx->nupdates=0;
331:     ctx->lmnow=0;

333:     
334:     return(0);
335: }


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;


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

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

370:     VecDot(ctx->Gprev,ctx->Xprev,&rhotemp); 
371:     VecDot(ctx->Gprev,ctx->Gprev,&y0temp); 

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

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

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

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

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

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

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

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

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

431:           sigmanew = ys_sum / yy_sum;
432:         }
433:         else {
434:           denom = 2*ctx->s_alpha*yy_sum;

436:           /*  Safeguard denom */
437:           if (0.0 == denom) {
438:             denom = TAO_ZERO_SAFEGUARD;
439:           }

441:           sigmanew = ((2*ctx->s_alpha-1)*ys_sum + 
442:                     PetscSqrtScalar((2*ctx->s_alpha-1)*(2*ctx->s_alpha-1)*ys_sum*ys_sum - 
443:                  4*(ctx->s_alpha)*(ctx->s_alpha-1)*yy_sum*ss_sum)) / denom;
444:         }

446:         switch(ctx->limitType) {
447:         case MatLMVM_Limit_Average:
448:           if (1.0 == ctx->mu) {
449:             ctx->sigma = sigmanew;
450:           }
451:           else if (ctx->mu) {
452:             ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
453:           }
454:           break;

456:         case MatLMVM_Limit_Relative:
457:           if (ctx->mu) {
458:             ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
459:           }
460:           break;

462:         case MatLMVM_Limit_Absolute:
463:           if (ctx->nu) {
464:             ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
465:           }
466:           break;

468:         default:
469:           ctx->sigma = sigmanew;
470:           break;
471:         }
472:         break;

474:       case MatLMVM_Scale_Broyden:
475:         /*  Original version */
476:         /*  Combine DFP and BFGS */

478:         /*  This code appears to be numerically unstable.  We use the */
479:         /*  original version because this was used to generate all of */
480:         /*  the data and because it may be the least unstable of the */
481:         /*  bunch. */

483:         /*  P = Q = inv(D); */
484:         VecCopy(ctx->D,ctx->P); 
485:         VecReciprocal(ctx->P); 
486:         VecCopy(ctx->P,ctx->Q); 

488:         /*  V = y*y */
489:         VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev); 

491:         /*  W = inv(D)*s */
492:         VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P); 
493:         VecDot(ctx->W,ctx->Xprev,&sDs); 

495:         /*  Safeguard rhotemp and sDs */
496:         if (0.0 == rhotemp) {
497:           rhotemp = TAO_ZERO_SAFEGUARD;
498:         }

500:         if (0.0 == sDs) {
501:           sDs = TAO_ZERO_SAFEGUARD;
502:         }

504:         if (1.0 != ctx->phi) {
505:           /*  BFGS portion of the update */
506:           /*  U = (inv(D)*s)*(inv(D)*s) */
507:           VecPointwiseMult(ctx->U,ctx->W,ctx->W); 


510:           /*  Assemble */
511:           VecAXPY(ctx->P,1.0/rhotemp,ctx->V); 
512:           VecAXPY(ctx->P,-1.0/sDs,ctx->U); 
513:         }

515:         if (0.0 != ctx->phi) {
516:           /*  DFP portion of the update */
517:           /*  U = inv(D)*s*y */
518:           VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev); 

520:           /*  Assemble */
521:           VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V); 
522:           VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U); 
523:         }

525:         if (0.0 == ctx->phi) {
526:             VecCopy(ctx->P,ctx->U); 
527:         }
528:         else if (1.0 == ctx->phi) {
529:             VecCopy(ctx->Q,ctx->U); 
530:         }
531:         else {
532:           /*  Broyden update U=(1-phi)*P + phi*Q */
533:             VecCopy(ctx->Q,ctx->U);
534:             VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P); 
535:         }

537:         /*  Obtain inverse and ensure positive definite */
538:         VecReciprocal(ctx->U); 
539:         VecAbs(ctx->U); 


542:         switch(ctx->rScaleType) {
543:         case MatLMVM_Rescale_None:
544:             break;

546:         case MatLMVM_Rescale_Scalar:
547:         case MatLMVM_Rescale_GL:
548:           if (ctx->rScaleType == MatLMVM_Rescale_GL) {
549:             /*  Gilbert and Lemarachal use the old diagonal */
550:             VecCopy(ctx->D,ctx->P); 
551:           }
552:           else {
553:             /*  The default version uses the current diagonal */
554:               VecCopy(ctx->U,ctx->P); 
555:           }

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

560:           /*  Save information for special cases of scalar rescaling */
561:           ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
562:           ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
563:           ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;

565:           if (0.5 == ctx->r_beta) {
566:             if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
567:               VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P); 
568:               VecDot(ctx->V,ctx->Y[0],&yy_sum); 
569:               
570:               VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P); 
571:               VecDot(ctx->W,ctx->S[0],&ss_sum); 

573:               ys_sum = ctx->ys_rhistory[0];
574:             }
575:             else {
576:               VecCopy(ctx->P,ctx->Q); 
577:               VecReciprocal(ctx->Q); 

579:               /*  Compute summations for scalar scaling */
580:               yy_sum = 0;        /*  No safeguard required */
581:               ys_sum = 0;        /*  No safeguard required */
582:               ss_sum = 0;        /*  No safeguard required */
583:               for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
584:                 VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P); 
585:                 VecDot(ctx->V,ctx->Y[i],&yDy); 
586:                 yy_sum += yDy;

588:                 VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q); 
589:                 VecDot(ctx->W,ctx->S[i],&sDs); 
590:                 ss_sum += sDs;
591:                 ys_sum += ctx->ys_rhistory[i];
592:               }
593:             }
594:           }
595:           else if (0.0 == ctx->r_beta) {
596:             if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
597:               /*  Compute summations for scalar scaling */
598:               VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P); 

600:               VecDot(ctx->W, ctx->Y[0], &ys_sum); 
601:               VecDot(ctx->W, ctx->W, &ss_sum); 
602:               yy_sum += ctx->yy_rhistory[0];
603:             }
604:             else {
605:               VecCopy(ctx->Q, ctx->P); 
606:               VecReciprocal(ctx->Q); 

608:               /*  Compute summations for scalar scaling */
609:               yy_sum = 0;        /*  No safeguard required */
610:               ys_sum = 0;        /*  No safeguard required */
611:               ss_sum = 0;        /*  No safeguard required */
612:               for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
613:                 VecPointwiseMult(ctx->W, ctx->S[i], ctx->Q); 
614:                 VecDot(ctx->W, ctx->Y[i], &yDs); 
615:                 ys_sum += yDs;

617:                 VecDot(ctx->W, ctx->W, &sDs); 
618:                 ss_sum += sDs;
619:   
620:                 yy_sum += ctx->yy_rhistory[i];
621:               }
622:             }
623:           }
624:           else if (1.0 == ctx->r_beta) {
625:             /*  Compute summations for scalar scaling */
626:             yy_sum = 0;        /*  No safeguard required */
627:             ys_sum = 0;        /*  No safeguard required */
628:             ss_sum = 0;        /*  No safeguard required */
629:             for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
630:               VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P); 
631:               VecDot(ctx->V, ctx->S[i], &yDs); 
632:               ys_sum += yDs;

634:               VecDot(ctx->V, ctx->V, &yDy); 
635:               yy_sum += yDy;

637:               ss_sum += ctx->ss_rhistory[i];
638:             }
639:           }
640:           else {
641:             VecCopy(ctx->Q, ctx->P); 

643:             VecPow(ctx->P, ctx->r_beta); 
644:             VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q); 

646:             /*  Compute summations for scalar scaling */
647:             yy_sum = 0;        /*  No safeguard required */
648:             ys_sum = 0;        /*  No safeguard required */
649:             ss_sum = 0;        /*  No safeguard required */
650:             for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
651:               VecPointwiseMult(ctx->V, ctx->P, ctx->Y[i]); 
652:               VecPointwiseMult(ctx->W, ctx->Q, ctx->S[i]); 

654:               VecDot(ctx->V, ctx->V, &yDy); 
655:               VecDot(ctx->V, ctx->W, &yDs); 
656:               VecDot(ctx->W, ctx->W, &sDs); 

658:               yy_sum += yDy;
659:               ys_sum += yDs;
660:               ss_sum += sDs;
661:             }
662:           }

664:           if (0.0 == ctx->r_alpha) {
665:             /*  Safeguard ys_sum  */
666:             if (0.0 == ys_sum) {
667:               ys_sum = TAO_ZERO_SAFEGUARD;
668:             }

670:             sigmanew = ss_sum / ys_sum;
671:           }
672:           else if (1.0 == ctx->r_alpha) {
673:             /*  Safeguard yy_sum  */
674:             if (0.0 == yy_sum) {
675:               ys_sum = TAO_ZERO_SAFEGUARD;
676:             }

678:             sigmanew = ys_sum / yy_sum;
679:           }
680:           else {
681:             denom = 2*ctx->r_alpha*yy_sum;

683:             /*  Safeguard denom */
684:             if (0.0 == denom) {
685:               denom = TAO_ZERO_SAFEGUARD;
686:             }

688:             sigmanew = ((2*ctx->r_alpha-1)*ys_sum +
689:                         PetscSqrtScalar((2*ctx->r_alpha-1)*(2*ctx->r_alpha-1)*ys_sum*ys_sum -
690:                              4*ctx->r_alpha*(ctx->r_alpha-1)*yy_sum*ss_sum)) / denom;
691:           }

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

698:           if (PetscIsInfOrNanReal(sigmanew)) {
699:             /*  sigmanew is not-a-number; skip rescaling */
700:           }
701:           else if (!sigmanew) {
702:             /*  sigmanew is zero; this is a bad case; skip rescaling */
703:           }
704:           else {
705:             /*  sigmanew is positive */
706:             VecScale(ctx->U, sigmanew); 
707:           }
708:           break;
709:         }

711:         /*  Modify for previous information */
712:         switch(ctx->limitType) {
713:         case MatLMVM_Limit_Average:
714:           if (1.0 == ctx->mu) {
715:             VecCopy(ctx->D, ctx->U); 
716:           }
717:           else if (ctx->mu) {
718:             VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U); 
719:           }
720:           break;
721:  
722:         case MatLMVM_Limit_Relative:
723:           if (ctx->mu) {
724:             /*  P = (1-mu) * D */
725:             VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D); 
726:             /*  Q = (1+mu) * D */
727:             VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D); 
728:             VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D); 
729:           }
730:           break;

732:         case MatLMVM_Limit_Absolute:
733:           if (ctx->nu) {
734:             VecCopy(ctx->P, ctx->D); 
735:             VecShift(ctx->P, -ctx->nu); 
736:             VecCopy(ctx->D, ctx->Q); 
737:             VecShift(ctx->Q, ctx->nu); 
738:             VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P); 
739:           }
740:           break;

742:         default:
743:             VecCopy(ctx->U, ctx->D); 
744:           break;
745:         } 
746:         break;
747:       }
748:       PetscObjectDereference((PetscObject)ctx->Xprev); 
749:       PetscObjectDereference((PetscObject)ctx->Gprev); 
750:       ctx->Xprev = ctx->S[ctx->lm]; 
751:       ctx->Gprev = ctx->Y[ctx->lm];
752:       PetscObjectReference((PetscObject)ctx->S[ctx->lm]); 
753:       PetscObjectReference((PetscObject)ctx->Y[ctx->lm]); 

755:     } 
756:     else { 
757:       ++ctx->nrejects;
758:     }
759:   }
760:   
761:   ++ctx->iter;
762:   VecCopy(x, ctx->Xprev); 
763:   VecCopy(g, ctx->Gprev); 
764:   return(0);
765: }

770: {
771:     MatLMVMCtx *ctx;
773:     PetscBool same;
776:     PetscTypeCompare((PetscObject)m,MATSHELL,&same); 
777:     if (!same) {
778:         SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
779:     }
780:     MatShellGetContext(m,(void**)&ctx); 
781:     ctx->delta = PetscAbsReal(d);
782:     ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
783:     ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
784:     return(0);
785:     
786: }

791: {
792:     MatLMVMCtx *ctx;
794:     PetscBool same;
797:     PetscTypeCompare((PetscObject)m,MATSHELL,&same); 
798:     if (!same) {
799:         SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
800:     }
801:     MatShellGetContext(m,(void**)&ctx); 
802:     
803:     if (ctx->scale) {
804:       VecDestroy(&ctx->scale); 
805:     }
806:     if (s) {
807:       VecDuplicate(s,&ctx->scale); 
808:     } else {
809:       ctx->scale = PETSC_NULL;
810:     }
811:     return(0);
812: }

817: {
818:     MatLMVMCtx *ctx;
820:     PetscBool same;
823:     PetscTypeCompare((PetscObject)m,MATSHELL,&same); 
824:     if (!same) {
825:         SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
826:     }
827:     MatShellGetContext(m,(void**)&ctx); 
828:     *nrejects = ctx->nrejects;
829:     return(0);
830: }

835: {
837:     return(0);
838: }

843: {
845:     return(0);
846: }

851: {
852:   MatLMVMCtx *ctx;
854:   PetscBool same;


860:   PetscTypeCompare((PetscObject)M,MATSHELL,&same); 
861:   if (!same) {SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");}
862:   MatShellGetContext(M,(void**)&ctx); 
863:   if (ctx->nupdates == 0) {
864:     MatLMVMUpdate(M,x,g); 
865:   } else {
866:     VecCopy(x,ctx->Xprev); 
867:     VecCopy(g,ctx->Gprev); 
868:     /*  TODO scaling specific terms */
869:   }
870:   return(0);
871:   
872: }
876: {
878:     PetscBool same;
884:     PetscTypeCompare((PetscObject)coarse,MATSHELL,&same); 
885:     if (!same) {
886:         SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
887:     }
888:     PetscTypeCompare((PetscObject)op,MATSHELL,&same); 
889:     if (!same) {
890:         SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
891:     }

893:     return(0);
894: }

899: {
901:     MatLMVMCtx *ctx;
902:     PetscBool same;

907:     PetscTypeCompare((PetscObject)m,MATSHELL,&same); 
908:     if (!same) {
909:         SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
910:     }
911:     MatShellGetContext(m,(void**)&ctx); 
912:     

914:     /*  Perform allocations */
915:     VecDuplicateVecs(v,ctx->lm+1,&ctx->S); 
916:     VecDuplicateVecs(v,ctx->lm+1,&ctx->Y); 
917:     VecDuplicate(v,&ctx->D); 
918:     VecDuplicate(v,&ctx->U); 
919:     VecDuplicate(v,&ctx->V); 
920:     VecDuplicate(v,&ctx->W); 
921:     VecDuplicate(v,&ctx->P); 
922:     VecDuplicate(v,&ctx->Q); 
923:     ctx->allocated = PETSC_TRUE;

925:     return(0);
926: }