Actual source code: lmvmmat.c

petsc-3.6.4 2016-04-12
Report Typos and Errors
  1: #include <../src/tao/matrix/lmvmmat.h>   /*I "lmvmmat.h" */
  2: #include <petsctao.h>  /*I "petsctao.h" */
  3: #include <petscksp.h>

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

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

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

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

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

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

 28:   Collective on A

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

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

 38:   Level: developer

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

 48:   /*  create data structure and populate with default values */
 49:   PetscNew(&ctx);
 50:   ctx->lm=5;
 51:   ctx->eps=0.0;
 52:   ctx->limitType=MatLMVM_Limit_None;
 53:   ctx->scaleType=MatLMVM_Scale_Broyden;
 54:   ctx->rScaleType = MatLMVM_Rescale_Scalar;
 55:   ctx->s_alpha = 1.0;
 56:   ctx->r_alpha = 1.0;
 57:   ctx->r_beta = 0.5;
 58:   ctx->mu = 1.0;
 59:   ctx->nu = 100.0;
 60: 
 61:   ctx->phi = 0.125;
 62: 
 63:   ctx->scalar_history = 1;
 64:   ctx->rescale_history = 1;
 65: 
 66:   ctx->delta_min = 1e-7;
 67:   ctx->delta_max = 100.0;

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

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

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

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

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

110:   ctx->Gprev = 0;
111:   ctx->Xprev = 0;

113:   ctx->scale = 0;
114:   ctx->useScale = PETSC_FALSE;

116:   ctx->H0 = 0;
117:   ctx->useDefaultH0=PETSC_TRUE;

119:   MatCreateShell(comm, n, n, N, N, ctx, A);
120:   MatShellSetOperation(*A,MATOP_DESTROY,(void(*)(void))MatDestroy_LMVM);
121:   MatShellSetOperation(*A,MATOP_VIEW,(void(*)(void))MatView_LMVM);
122:   return(0);
123: }

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

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

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

151:   scaled = PETSC_FALSE;
152:   if (!scaled && !shell->useDefaultH0 && shell->H0) {
153:     MatSolve(shell->H0,x,shell->U);
154:     VecDot(x,shell->U,&dd);
155:     if ((dd > 0.0) && !PetscIsInfOrNanReal(dd)) {
156:       /*  Accept Hessian solve */
157:       VecCopy(shell->U,x);
158:       scaled = PETSC_TRUE;
159:     }
160:   }

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

172:   if (!scaled) {
173:     switch(shell->scaleType) {
174:     case MatLMVM_Scale_None:
175:       break;

177:     case MatLMVM_Scale_Scalar:
178:       VecScale(x,shell->sigma);
179:       break;

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

195: extern PetscErrorCode MatView_LMVM(Mat A, PetscViewer pv)
196: {
197:   PetscBool      isascii;
199:   MatLMVMCtx     *lmP;

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

218: extern 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:     if (ctx->scale) {
242:       VecDestroy(&ctx->scale);
243:     }
244:   }
245:   PetscFree(ctx->rho);
246:   PetscFree(ctx->beta);
247:   PetscFree(ctx->yy_history);
248:   PetscFree(ctx->ys_history);
249:   PetscFree(ctx->ss_history);
250:   PetscFree(ctx->yy_rhistory);
251:   PetscFree(ctx->ys_rhistory);
252:   PetscFree(ctx->ss_rhistory);
253:   PetscFree(ctx);
254:   return(0);
255: }

259: extern PetscErrorCode MatLMVMReset(Mat M)
260: {
262:   MatLMVMCtx     *ctx;
263:   PetscInt       i;

266:   MatShellGetContext(M,(void**)&ctx);
267:   if (ctx->Gprev) {
268:     PetscObjectDereference((PetscObject)ctx->Gprev);
269:   }
270:   if (ctx->Xprev) {
271:     PetscObjectDereference((PetscObject)ctx->Xprev);
272:   }
273:   ctx->Gprev = ctx->Y[ctx->lm];
274:   ctx->Xprev = ctx->S[ctx->lm];
275:   PetscObjectReference((PetscObject)ctx->Gprev);
276:   PetscObjectReference((PetscObject)ctx->Xprev);
277:   for (i=0; i<ctx->lm; ++i) {
278:     ctx->rho[i] = 0.0;
279:   }
280:   ctx->rho[0] = 1.0;

282:   /*  Set the scaling and diagonal scaling matrix */
283:   switch(ctx->scaleType) {
284:   case MatLMVM_Scale_None:
285:     ctx->sigma = 1.0;
286:     break;
287:   case MatLMVM_Scale_Scalar:
288:     ctx->sigma = ctx->delta;
289:     break;
290:   case MatLMVM_Scale_Broyden:
291:     VecSet(ctx->D,ctx->delta);
292:     break;
293:   }

295:   ctx->iter=0;
296:   ctx->nupdates=0;
297:   ctx->lmnow=0;
298:   return(0);
299: }

303: extern PetscErrorCode MatLMVMUpdate(Mat M, Vec x, Vec g)
304: {
305:   MatLMVMCtx     *ctx;
306:   PetscReal      rhotemp, rhotol;
307:   PetscReal      y0temp, s0temp;
308:   PetscReal      yDy, yDs, sDs;
309:   PetscReal      sigmanew, denom;
311:   PetscInt       i;
312:   PetscBool      same;
313:   PetscReal      yy_sum=0.0, ys_sum=0.0, ss_sum=0.0;

318:   PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
319:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
320:   MatShellGetContext(M,(void**)&ctx);
321:   if (!ctx->allocated) {
322:     MatLMVMAllocateVectors(M, x);
323:   }

325:   if (0 == ctx->iter) {
326:     MatLMVMReset(M);
327:   }  else {
328:     VecAYPX(ctx->Gprev,-1.0,g);
329:     VecAYPX(ctx->Xprev,-1.0,x);

331:     VecDot(ctx->Gprev,ctx->Xprev,&rhotemp);
332:     VecDot(ctx->Gprev,ctx->Gprev,&y0temp);

334:     rhotol = ctx->eps * y0temp;
335:     if (rhotemp > rhotol) {
336:       ++ctx->nupdates;

338:       ctx->lmnow = PetscMin(ctx->lmnow+1, ctx->lm);
339:       ierr=PetscObjectDereference((PetscObject)ctx->S[ctx->lm]);
340:       ierr=PetscObjectDereference((PetscObject)ctx->Y[ctx->lm]);
341:       for (i = ctx->lm-1; i >= 0; --i) {
342:         ctx->S[i+1] = ctx->S[i];
343:         ctx->Y[i+1] = ctx->Y[i];
344:         ctx->rho[i+1] = ctx->rho[i];
345:       }
346:       ctx->S[0] = ctx->Xprev;
347:       ctx->Y[0] = ctx->Gprev;
348:       PetscObjectReference((PetscObject)ctx->S[0]);
349:       PetscObjectReference((PetscObject)ctx->Y[0]);
350:       ctx->rho[0] = 1.0 / rhotemp;

352:       /*  Compute the scaling */
353:       switch(ctx->scaleType) {
354:       case MatLMVM_Scale_None:
355:         break;

357:       case MatLMVM_Scale_Scalar:
358:         /*  Compute s^T s  */
359:           VecDot(ctx->Xprev,ctx->Xprev,&s0temp);

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

363:         /*  Save information for scalar scaling */
364:         ctx->yy_history[(ctx->nupdates - 1) % ctx->scalar_history] = y0temp;
365:         ctx->ys_history[(ctx->nupdates - 1) % ctx->scalar_history] = rhotemp;
366:         ctx->ss_history[(ctx->nupdates - 1) % ctx->scalar_history] = s0temp;

368:         /*  Compute summations for scalar scaling */
369:         yy_sum = 0;     /*  No safeguard required; y^T y > 0 */
370:         ys_sum = 0;     /*  No safeguard required; y^T s > 0 */
371:         ss_sum = 0;     /*  No safeguard required; s^T s > 0 */
372:         for (i = 0; i < PetscMin(ctx->nupdates, ctx->scalar_history); ++i) {
373:           yy_sum += ctx->yy_history[i];
374:           ys_sum += ctx->ys_history[i];
375:           ss_sum += ctx->ss_history[i];
376:         }

378:         if (0.0 == ctx->s_alpha) {
379:           /*  Safeguard ys_sum  */
380:           if (0.0 == ys_sum) {
381:             ys_sum = TAO_ZERO_SAFEGUARD;
382:           }

384:           sigmanew = ss_sum / ys_sum;
385:         } else if (1.0 == ctx->s_alpha) {
386:           /*  Safeguard yy_sum  */
387:           if (0.0 == yy_sum) {
388:             yy_sum = TAO_ZERO_SAFEGUARD;
389:           }

391:           sigmanew = ys_sum / yy_sum;
392:         } else {
393:           denom = 2*ctx->s_alpha*yy_sum;

395:           /*  Safeguard denom */
396:           if (0.0 == denom) {
397:             denom = TAO_ZERO_SAFEGUARD;
398:           }

400:           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;
401:         }

403:         switch(ctx->limitType) {
404:         case MatLMVM_Limit_Average:
405:           if (1.0 == ctx->mu) {
406:             ctx->sigma = sigmanew;
407:           } else if (ctx->mu) {
408:             ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
409:           }
410:           break;

412:         case MatLMVM_Limit_Relative:
413:           if (ctx->mu) {
414:             ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
415:           }
416:           break;

418:         case MatLMVM_Limit_Absolute:
419:           if (ctx->nu) {
420:             ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
421:           }
422:           break;

424:         default:
425:           ctx->sigma = sigmanew;
426:           break;
427:         }
428:         break;

430:       case MatLMVM_Scale_Broyden:
431:         /*  Original version */
432:         /*  Combine DFP and BFGS */

434:         /*  This code appears to be numerically unstable.  We use the */
435:         /*  original version because this was used to generate all of */
436:         /*  the data and because it may be the least unstable of the */
437:         /*  bunch. */

439:         /*  P = Q = inv(D); */
440:         VecCopy(ctx->D,ctx->P);
441:         VecReciprocal(ctx->P);
442:         VecCopy(ctx->P,ctx->Q);

444:         /*  V = y*y */
445:         VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev);

447:         /*  W = inv(D)*s */
448:         VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P);
449:         VecDot(ctx->W,ctx->Xprev,&sDs);

451:         /*  Safeguard rhotemp and sDs */
452:         if (0.0 == rhotemp) {
453:           rhotemp = TAO_ZERO_SAFEGUARD;
454:         }

456:         if (0.0 == sDs) {
457:           sDs = TAO_ZERO_SAFEGUARD;
458:         }

460:         if (1.0 != ctx->phi) {
461:           /*  BFGS portion of the update */
462:           /*  U = (inv(D)*s)*(inv(D)*s) */
463:           VecPointwiseMult(ctx->U,ctx->W,ctx->W);

465:           /*  Assemble */
466:           VecAXPY(ctx->P,1.0/rhotemp,ctx->V);
467:           VecAXPY(ctx->P,-1.0/sDs,ctx->U);
468:         }

470:         if (0.0 != ctx->phi) {
471:           /*  DFP portion of the update */
472:           /*  U = inv(D)*s*y */
473:           VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev);

475:           /*  Assemble */
476:           VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V);
477:           VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U);
478:         }

480:         if (0.0 == ctx->phi) {
481:             VecCopy(ctx->P,ctx->U);
482:         } else if (1.0 == ctx->phi) {
483:             VecCopy(ctx->Q,ctx->U);
484:         } else {
485:           /*  Broyden update U=(1-phi)*P + phi*Q */
486:             VecCopy(ctx->Q,ctx->U);
487:             VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P);
488:         }

490:         /*  Obtain inverse and ensure positive definite */
491:         VecReciprocal(ctx->U);
492:         VecAbs(ctx->U);

494:         switch(ctx->rScaleType) {
495:         case MatLMVM_Rescale_None:
496:             break;

498:         case MatLMVM_Rescale_Scalar:
499:         case MatLMVM_Rescale_GL:
500:           if (ctx->rScaleType == MatLMVM_Rescale_GL) {
501:             /*  Gilbert and Lemarachal use the old diagonal */
502:             VecCopy(ctx->D,ctx->P);
503:           } else {
504:             /*  The default version uses the current diagonal */
505:               VecCopy(ctx->U,ctx->P);
506:           }

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

511:           /*  Save information for special cases of scalar rescaling */
512:           ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
513:           ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
514:           ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;

516:           if (0.5 == ctx->r_beta) {
517:             if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
518:               VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P);
519:               VecDot(ctx->V,ctx->Y[0],&yy_sum);

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

524:               ys_sum = ctx->ys_rhistory[0];
525:             } else {
526:               VecCopy(ctx->P,ctx->Q);
527:               VecReciprocal(ctx->Q);

529:               /*  Compute summations for scalar scaling */
530:               yy_sum = 0;       /*  No safeguard required */
531:               ys_sum = 0;       /*  No safeguard required */
532:               ss_sum = 0;       /*  No safeguard required */
533:               for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
534:                 VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P);
535:                 VecDot(ctx->V,ctx->Y[i],&yDy);
536:                 yy_sum += yDy;

538:                 VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q);
539:                 VecDot(ctx->W,ctx->S[i],&sDs);
540:                 ss_sum += sDs;
541:                 ys_sum += ctx->ys_rhistory[i];
542:               }
543:             }
544:           } else if (0.0 == ctx->r_beta) {
545:             if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
546:               /*  Compute summations for scalar scaling */
547:               VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);

549:               VecDot(ctx->W, ctx->Y[0], &ys_sum);
550:               VecDot(ctx->W, ctx->W, &ss_sum);
551:               yy_sum += ctx->yy_rhistory[0];
552:             } else {
553:               VecCopy(ctx->Q, ctx->P);
554:               VecReciprocal(ctx->Q);

556:               /*  Compute summations for scalar scaling */
557:               yy_sum = 0;       /*  No safeguard required */
558:               ys_sum = 0;       /*  No safeguard required */
559:               ss_sum = 0;       /*  No safeguard required */
560:               for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
561:                 VecPointwiseMult(ctx->W, ctx->S[i], ctx->Q);
562:                 VecDot(ctx->W, ctx->Y[i], &yDs);
563:                 ys_sum += yDs;

565:                 VecDot(ctx->W, ctx->W, &sDs);
566:                 ss_sum += sDs;

568:                 yy_sum += ctx->yy_rhistory[i];
569:               }
570:             }
571:           } else if (1.0 == ctx->r_beta) {
572:             /*  Compute summations for scalar scaling */
573:             yy_sum = 0; /*  No safeguard required */
574:             ys_sum = 0; /*  No safeguard required */
575:             ss_sum = 0; /*  No safeguard required */
576:             for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
577:               VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P);
578:               VecDot(ctx->V, ctx->S[i], &yDs);
579:               ys_sum += yDs;

581:               VecDot(ctx->V, ctx->V, &yDy);
582:               yy_sum += yDy;

584:               ss_sum += ctx->ss_rhistory[i];
585:             }
586:           } else {
587:             VecCopy(ctx->Q, ctx->P);

589:             VecPow(ctx->P, ctx->r_beta);
590:             VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q);

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

600:               VecDot(ctx->V, ctx->V, &yDy);
601:               VecDot(ctx->V, ctx->W, &yDs);
602:               VecDot(ctx->W, ctx->W, &sDs);

604:               yy_sum += yDy;
605:               ys_sum += yDs;
606:               ss_sum += sDs;
607:             }
608:           }

610:           if (0.0 == ctx->r_alpha) {
611:             /*  Safeguard ys_sum  */
612:             if (0.0 == ys_sum) {
613:               ys_sum = TAO_ZERO_SAFEGUARD;
614:             }

616:             sigmanew = ss_sum / ys_sum;
617:           } else if (1.0 == ctx->r_alpha) {
618:             /*  Safeguard yy_sum  */
619:             if (0.0 == yy_sum) {
620:               ys_sum = TAO_ZERO_SAFEGUARD;
621:             }

623:             sigmanew = ys_sum / yy_sum;
624:           } else {
625:             denom = 2*ctx->r_alpha*yy_sum;

627:             /*  Safeguard denom */
628:             if (0.0 == denom) {
629:               denom = TAO_ZERO_SAFEGUARD;
630:             }

632:             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;
633:           }

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

640:           if (PetscIsInfOrNanReal(sigmanew)) {
641:             /*  sigmanew is not-a-number; skip rescaling */
642:           } else if (!sigmanew) {
643:             /*  sigmanew is zero; this is a bad case; skip rescaling */
644:           } else {
645:             /*  sigmanew is positive */
646:             VecScale(ctx->U, sigmanew);
647:           }
648:           break;
649:         }

651:         /*  Modify for previous information */
652:         switch(ctx->limitType) {
653:         case MatLMVM_Limit_Average:
654:           if (1.0 == ctx->mu) {
655:             VecCopy(ctx->D, ctx->U);
656:           } else if (ctx->mu) {
657:             VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U);
658:           }
659:           break;

661:         case MatLMVM_Limit_Relative:
662:           if (ctx->mu) {
663:             /*  P = (1-mu) * D */
664:             VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D);
665:             /*  Q = (1+mu) * D */
666:             VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D);
667:             VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D);
668:           }
669:           break;

671:         case MatLMVM_Limit_Absolute:
672:           if (ctx->nu) {
673:             VecCopy(ctx->P, ctx->D);
674:             VecShift(ctx->P, -ctx->nu);
675:             VecCopy(ctx->D, ctx->Q);
676:             VecShift(ctx->Q, ctx->nu);
677:             VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P);
678:           }
679:           break;

681:         default:
682:             VecCopy(ctx->U, ctx->D);
683:           break;
684:         }
685:         break;
686:       }
687:       PetscObjectDereference((PetscObject)ctx->Xprev);
688:       PetscObjectDereference((PetscObject)ctx->Gprev);
689:       ctx->Xprev = ctx->S[ctx->lm];
690:       ctx->Gprev = ctx->Y[ctx->lm];
691:       PetscObjectReference((PetscObject)ctx->S[ctx->lm]);
692:       PetscObjectReference((PetscObject)ctx->Y[ctx->lm]);

694:     } else {
695:       ++ctx->nrejects;
696:     }
697:   }

699:   ++ctx->iter;
700:   VecCopy(x, ctx->Xprev);
701:   VecCopy(g, ctx->Gprev);
702:   return(0);
703: }

707: extern PetscErrorCode MatLMVMSetDelta(Mat m, PetscReal d)
708: {
709:   MatLMVMCtx     *ctx;
711:   PetscBool      same;

715:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
716:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
717:   MatShellGetContext(m,(void**)&ctx);
718:   ctx->delta = PetscAbsReal(d);
719:   ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
720:   ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
721:   return(0);
722: }

726: extern PetscErrorCode MatLMVMSetScale(Mat m, Vec s)
727: {
728:   MatLMVMCtx     *ctx;
730:   PetscBool      same;

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

738:   VecDestroy(&ctx->scale);
739:   if (s) {
740:     VecDuplicate(s,&ctx->scale);
741:   }
742:   return(0);
743: }

747: extern PetscErrorCode MatLMVMGetRejects(Mat m, PetscInt *nrejects)
748: {
749:   MatLMVMCtx     *ctx;
751:   PetscBool      same;

755:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
756:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
757:   MatShellGetContext(m,(void**)&ctx);
758:   *nrejects = ctx->nrejects;
759:   return(0);
760: }

764: extern PetscErrorCode MatLMVMSetH0(Mat m, Mat A)
765: {
767:     return(0);
768: }

772: extern PetscErrorCode MatLMVMGetX0(Mat m, Vec x)
773: {
775:     return(0);
776: }

780: extern PetscErrorCode MatLMVMSetPrev(Mat M, Vec x, Vec g)
781: {
782:   MatLMVMCtx     *ctx;
784:   PetscBool      same;

789:   PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
790:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
791:   MatShellGetContext(M,(void**)&ctx);
792:   if (ctx->nupdates == 0) {
793:     MatLMVMUpdate(M,x,g);
794:   } else {
795:     VecCopy(x,ctx->Xprev);
796:     VecCopy(g,ctx->Gprev);
797:     /*  TODO scaling specific terms */
798:   }
799:   return(0);
800: }

804: extern PetscErrorCode MatLMVMRefine(Mat coarse, Mat op, Vec fineX, Vec fineG)
805: {
807:   PetscBool      same;

814:   PetscObjectTypeCompare((PetscObject)coarse,MATSHELL,&same);
815:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
816:   PetscObjectTypeCompare((PetscObject)op,MATSHELL,&same);
817:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
818:   return(0);
819: }

823: extern PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v)
824: {
826:   MatLMVMCtx     *ctx;
827:   PetscBool      same;

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

836:   /*  Perform allocations */
837:   VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
838:   VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
839:   VecDuplicate(v,&ctx->D);
840:   VecDuplicate(v,&ctx->U);
841:   VecDuplicate(v,&ctx->V);
842:   VecDuplicate(v,&ctx->W);
843:   VecDuplicate(v,&ctx->P);
844:   VecDuplicate(v,&ctx->Q);
845:   ctx->allocated = PETSC_TRUE;
846:   return(0);
847: }