Actual source code: lmvmmat.c

petsc-3.7.3 2016-08-01
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_mat = 0;
117:   ctx->H0_ksp = 0;
118:   ctx->H0_norm = 0;
119:   ctx->useDefaultH0 = PETSC_TRUE;

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

129: extern 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: }

199: extern PetscErrorCode MatView_LMVM(Mat A, PetscViewer pv)
200: {
201:   PetscBool      isascii;
203:   MatLMVMCtx     *lmP;

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

222: extern PetscErrorCode MatDestroy_LMVM(Mat M)
223: {
224:   MatLMVMCtx     *ctx;

228:   MatShellGetContext(M,(void**)&ctx);
229:   if (ctx->allocated) {
230:     if (ctx->Xprev) {
231:       PetscObjectDereference((PetscObject)ctx->Xprev);
232:     }
233:     if (ctx->Gprev) {
234:       PetscObjectDereference((PetscObject)ctx->Gprev);
235:     }

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

268: extern PetscErrorCode MatLMVMReset(Mat M)
269: {
271:   MatLMVMCtx     *ctx;
272:   PetscInt       i;

275:   MatShellGetContext(M,(void**)&ctx);
276:   if (ctx->Gprev) {
277:     PetscObjectDereference((PetscObject)ctx->Gprev);
278:   }
279:   if (ctx->Xprev) {
280:     PetscObjectDereference((PetscObject)ctx->Xprev);
281:   }
282:   ctx->Gprev = ctx->Y[ctx->lm];
283:   ctx->Xprev = ctx->S[ctx->lm];
284:   PetscObjectReference((PetscObject)ctx->Gprev);
285:   PetscObjectReference((PetscObject)ctx->Xprev);
286:   for (i=0; i<ctx->lm; ++i) {
287:     ctx->rho[i] = 0.0;
288:   }
289:   ctx->rho[0] = 1.0;

291:   /*  Set the scaling and diagonal scaling matrix */
292:   switch(ctx->scaleType) {
293:   case MatLMVM_Scale_None:
294:     ctx->sigma = 1.0;
295:     break;
296:   case MatLMVM_Scale_Scalar:
297:     ctx->sigma = ctx->delta;
298:     break;
299:   case MatLMVM_Scale_Broyden:
300:     VecSet(ctx->D,ctx->delta);
301:     break;
302:   }

304:   ctx->iter=0;
305:   ctx->nupdates=0;
306:   ctx->lmnow=0;
307:   return(0);
308: }

312: extern PetscErrorCode MatLMVMUpdate(Mat M, Vec x, Vec g)
313: {
314:   MatLMVMCtx     *ctx;
315:   PetscReal      rhotemp, rhotol;
316:   PetscReal      y0temp, s0temp;
317:   PetscReal      yDy, yDs, sDs;
318:   PetscReal      sigmanew, denom;
320:   PetscInt       i;
321:   PetscBool      same;
322:   PetscReal      yy_sum=0.0, ys_sum=0.0, ss_sum=0.0;

327:   PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
328:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
329:   MatShellGetContext(M,(void**)&ctx);
330:   if (!ctx->allocated) {
331:     MatLMVMAllocateVectors(M, x);
332:   }

334:   if (0 == ctx->iter) {
335:     MatLMVMReset(M);
336:   }  else {
337:     VecAYPX(ctx->Gprev,-1.0,g);
338:     VecAYPX(ctx->Xprev,-1.0,x);

340:     VecDot(ctx->Gprev,ctx->Xprev,&rhotemp);
341:     VecDot(ctx->Gprev,ctx->Gprev,&y0temp);

343:     rhotol = ctx->eps * y0temp;
344:     if (rhotemp > rhotol) {
345:       ++ctx->nupdates;

347:       ctx->lmnow = PetscMin(ctx->lmnow+1, ctx->lm);
348:       ierr=PetscObjectDereference((PetscObject)ctx->S[ctx->lm]);
349:       ierr=PetscObjectDereference((PetscObject)ctx->Y[ctx->lm]);
350:       for (i = ctx->lm-1; i >= 0; --i) {
351:         ctx->S[i+1] = ctx->S[i];
352:         ctx->Y[i+1] = ctx->Y[i];
353:         ctx->rho[i+1] = ctx->rho[i];
354:       }
355:       ctx->S[0] = ctx->Xprev;
356:       ctx->Y[0] = ctx->Gprev;
357:       PetscObjectReference((PetscObject)ctx->S[0]);
358:       PetscObjectReference((PetscObject)ctx->Y[0]);
359:       ctx->rho[0] = 1.0 / rhotemp;

361:       /*  Compute the scaling */
362:       switch(ctx->scaleType) {
363:       case MatLMVM_Scale_None:
364:         break;

366:       case MatLMVM_Scale_Scalar:
367:         /*  Compute s^T s  */
368:           VecDot(ctx->Xprev,ctx->Xprev,&s0temp);

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

372:         /*  Save information for scalar scaling */
373:         ctx->yy_history[(ctx->nupdates - 1) % ctx->scalar_history] = y0temp;
374:         ctx->ys_history[(ctx->nupdates - 1) % ctx->scalar_history] = rhotemp;
375:         ctx->ss_history[(ctx->nupdates - 1) % ctx->scalar_history] = s0temp;

377:         /*  Compute summations for scalar scaling */
378:         yy_sum = 0;     /*  No safeguard required; y^T y > 0 */
379:         ys_sum = 0;     /*  No safeguard required; y^T s > 0 */
380:         ss_sum = 0;     /*  No safeguard required; s^T s > 0 */
381:         for (i = 0; i < PetscMin(ctx->nupdates, ctx->scalar_history); ++i) {
382:           yy_sum += ctx->yy_history[i];
383:           ys_sum += ctx->ys_history[i];
384:           ss_sum += ctx->ss_history[i];
385:         }

387:         if (0.0 == ctx->s_alpha) {
388:           /*  Safeguard ys_sum  */
389:           if (0.0 == ys_sum) {
390:             ys_sum = TAO_ZERO_SAFEGUARD;
391:           }

393:           sigmanew = ss_sum / ys_sum;
394:         } else if (1.0 == ctx->s_alpha) {
395:           /*  Safeguard yy_sum  */
396:           if (0.0 == yy_sum) {
397:             yy_sum = TAO_ZERO_SAFEGUARD;
398:           }

400:           sigmanew = ys_sum / yy_sum;
401:         } else {
402:           denom = 2*ctx->s_alpha*yy_sum;

404:           /*  Safeguard denom */
405:           if (0.0 == denom) {
406:             denom = TAO_ZERO_SAFEGUARD;
407:           }

409:           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;
410:         }

412:         switch(ctx->limitType) {
413:         case MatLMVM_Limit_Average:
414:           if (1.0 == ctx->mu) {
415:             ctx->sigma = sigmanew;
416:           } else if (ctx->mu) {
417:             ctx->sigma = ctx->mu * sigmanew + (1.0 - ctx->mu) * ctx->sigma;
418:           }
419:           break;

421:         case MatLMVM_Limit_Relative:
422:           if (ctx->mu) {
423:             ctx->sigma = TaoMid((1.0 - ctx->mu) * ctx->sigma, sigmanew, (1.0 + ctx->mu) * ctx->sigma);
424:           }
425:           break;

427:         case MatLMVM_Limit_Absolute:
428:           if (ctx->nu) {
429:             ctx->sigma = TaoMid(ctx->sigma - ctx->nu, sigmanew, ctx->sigma + ctx->nu);
430:           }
431:           break;

433:         default:
434:           ctx->sigma = sigmanew;
435:           break;
436:         }
437:         break;

439:       case MatLMVM_Scale_Broyden:
440:         /*  Original version */
441:         /*  Combine DFP and BFGS */

443:         /*  This code appears to be numerically unstable.  We use the */
444:         /*  original version because this was used to generate all of */
445:         /*  the data and because it may be the least unstable of the */
446:         /*  bunch. */

448:         /*  P = Q = inv(D); */
449:         VecCopy(ctx->D,ctx->P);
450:         VecReciprocal(ctx->P);
451:         VecCopy(ctx->P,ctx->Q);

453:         /*  V = y*y */
454:         VecPointwiseMult(ctx->V,ctx->Gprev,ctx->Gprev);

456:         /*  W = inv(D)*s */
457:         VecPointwiseMult(ctx->W,ctx->Xprev,ctx->P);
458:         VecDot(ctx->W,ctx->Xprev,&sDs);

460:         /*  Safeguard rhotemp and sDs */
461:         if (0.0 == rhotemp) {
462:           rhotemp = TAO_ZERO_SAFEGUARD;
463:         }

465:         if (0.0 == sDs) {
466:           sDs = TAO_ZERO_SAFEGUARD;
467:         }

469:         if (1.0 != ctx->phi) {
470:           /*  BFGS portion of the update */
471:           /*  U = (inv(D)*s)*(inv(D)*s) */
472:           VecPointwiseMult(ctx->U,ctx->W,ctx->W);

474:           /*  Assemble */
475:           VecAXPY(ctx->P,1.0/rhotemp,ctx->V);
476:           VecAXPY(ctx->P,-1.0/sDs,ctx->U);
477:         }

479:         if (0.0 != ctx->phi) {
480:           /*  DFP portion of the update */
481:           /*  U = inv(D)*s*y */
482:           VecPointwiseMult(ctx->U, ctx->W, ctx->Gprev);

484:           /*  Assemble */
485:           VecAXPY(ctx->Q,1.0/rhotemp + sDs/(rhotemp*rhotemp), ctx->V);
486:           VecAXPY(ctx->Q,-2.0/rhotemp,ctx->U);
487:         }

489:         if (0.0 == ctx->phi) {
490:             VecCopy(ctx->P,ctx->U);
491:         } else if (1.0 == ctx->phi) {
492:             VecCopy(ctx->Q,ctx->U);
493:         } else {
494:           /*  Broyden update U=(1-phi)*P + phi*Q */
495:             VecCopy(ctx->Q,ctx->U);
496:             VecAXPBY(ctx->U,1.0-ctx->phi, ctx->phi, ctx->P);
497:         }

499:         /*  Obtain inverse and ensure positive definite */
500:         VecReciprocal(ctx->U);
501:         VecAbs(ctx->U);

503:         switch(ctx->rScaleType) {
504:         case MatLMVM_Rescale_None:
505:             break;

507:         case MatLMVM_Rescale_Scalar:
508:         case MatLMVM_Rescale_GL:
509:           if (ctx->rScaleType == MatLMVM_Rescale_GL) {
510:             /*  Gilbert and Lemarachal use the old diagonal */
511:             VecCopy(ctx->D,ctx->P);
512:           } else {
513:             /*  The default version uses the current diagonal */
514:               VecCopy(ctx->U,ctx->P);
515:           }

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

520:           /*  Save information for special cases of scalar rescaling */
521:           ctx->yy_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = y0temp;
522:           ctx->ys_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = rhotemp;
523:           ctx->ss_rhistory[(ctx->nupdates - 1) % ctx->rescale_history] = s0temp;

525:           if (0.5 == ctx->r_beta) {
526:             if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
527:               VecPointwiseMult(ctx->V,ctx->Y[0],ctx->P);
528:               VecDot(ctx->V,ctx->Y[0],&yy_sum);

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

533:               ys_sum = ctx->ys_rhistory[0];
534:             } else {
535:               VecCopy(ctx->P,ctx->Q);
536:               VecReciprocal(ctx->Q);

538:               /*  Compute summations for scalar scaling */
539:               yy_sum = 0;       /*  No safeguard required */
540:               ys_sum = 0;       /*  No safeguard required */
541:               ss_sum = 0;       /*  No safeguard required */
542:               for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
543:                 VecPointwiseMult(ctx->V,ctx->Y[i],ctx->P);
544:                 VecDot(ctx->V,ctx->Y[i],&yDy);
545:                 yy_sum += yDy;

547:                 VecPointwiseMult(ctx->W,ctx->S[i],ctx->Q);
548:                 VecDot(ctx->W,ctx->S[i],&sDs);
549:                 ss_sum += sDs;
550:                 ys_sum += ctx->ys_rhistory[i];
551:               }
552:             }
553:           } else if (0.0 == ctx->r_beta) {
554:             if (1 == PetscMin(ctx->nupdates, ctx->rescale_history)) {
555:               /*  Compute summations for scalar scaling */
556:               VecPointwiseDivide(ctx->W,ctx->S[0],ctx->P);

558:               VecDot(ctx->W, ctx->Y[0], &ys_sum);
559:               VecDot(ctx->W, ctx->W, &ss_sum);
560:               yy_sum += ctx->yy_rhistory[0];
561:             } else {
562:               VecCopy(ctx->Q, ctx->P);
563:               VecReciprocal(ctx->Q);

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

574:                 VecDot(ctx->W, ctx->W, &sDs);
575:                 ss_sum += sDs;

577:                 yy_sum += ctx->yy_rhistory[i];
578:               }
579:             }
580:           } else if (1.0 == ctx->r_beta) {
581:             /*  Compute summations for scalar scaling */
582:             yy_sum = 0; /*  No safeguard required */
583:             ys_sum = 0; /*  No safeguard required */
584:             ss_sum = 0; /*  No safeguard required */
585:             for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
586:               VecPointwiseMult(ctx->V, ctx->Y[i], ctx->P);
587:               VecDot(ctx->V, ctx->S[i], &yDs);
588:               ys_sum += yDs;

590:               VecDot(ctx->V, ctx->V, &yDy);
591:               yy_sum += yDy;

593:               ss_sum += ctx->ss_rhistory[i];
594:             }
595:           } else {
596:             VecCopy(ctx->Q, ctx->P);

598:             VecPow(ctx->P, ctx->r_beta);
599:             VecPointwiseDivide(ctx->Q, ctx->P, ctx->Q);

601:             /*  Compute summations for scalar scaling */
602:             yy_sum = 0; /*  No safeguard required */
603:             ys_sum = 0; /*  No safeguard required */
604:             ss_sum = 0; /*  No safeguard required */
605:             for (i = 0; i < PetscMin(ctx->nupdates, ctx->rescale_history); ++i) {
606:               VecPointwiseMult(ctx->V, ctx->P, ctx->Y[i]);
607:               VecPointwiseMult(ctx->W, ctx->Q, ctx->S[i]);

609:               VecDot(ctx->V, ctx->V, &yDy);
610:               VecDot(ctx->V, ctx->W, &yDs);
611:               VecDot(ctx->W, ctx->W, &sDs);

613:               yy_sum += yDy;
614:               ys_sum += yDs;
615:               ss_sum += sDs;
616:             }
617:           }

619:           if (0.0 == ctx->r_alpha) {
620:             /*  Safeguard ys_sum  */
621:             if (0.0 == ys_sum) {
622:               ys_sum = TAO_ZERO_SAFEGUARD;
623:             }

625:             sigmanew = ss_sum / ys_sum;
626:           } else if (1.0 == ctx->r_alpha) {
627:             /*  Safeguard yy_sum  */
628:             if (0.0 == yy_sum) {
629:               ys_sum = TAO_ZERO_SAFEGUARD;
630:             }

632:             sigmanew = ys_sum / yy_sum;
633:           } else {
634:             denom = 2*ctx->r_alpha*yy_sum;

636:             /*  Safeguard denom */
637:             if (0.0 == denom) {
638:               denom = TAO_ZERO_SAFEGUARD;
639:             }

641:             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;
642:           }

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

649:           if (PetscIsInfOrNanReal(sigmanew)) {
650:             /*  sigmanew is not-a-number; skip rescaling */
651:           } else if (!sigmanew) {
652:             /*  sigmanew is zero; this is a bad case; skip rescaling */
653:           } else {
654:             /*  sigmanew is positive */
655:             VecScale(ctx->U, sigmanew);
656:           }
657:           break;
658:         }

660:         /*  Modify for previous information */
661:         switch(ctx->limitType) {
662:         case MatLMVM_Limit_Average:
663:           if (1.0 == ctx->mu) {
664:             VecCopy(ctx->D, ctx->U);
665:           } else if (ctx->mu) {
666:             VecAXPBY(ctx->D,ctx->mu, 1.0-ctx->mu,ctx->U);
667:           }
668:           break;

670:         case MatLMVM_Limit_Relative:
671:           if (ctx->mu) {
672:             /*  P = (1-mu) * D */
673:             VecAXPBY(ctx->P, 1.0-ctx->mu, 0.0, ctx->D);
674:             /*  Q = (1+mu) * D */
675:             VecAXPBY(ctx->Q, 1.0+ctx->mu, 0.0, ctx->D);
676:             VecMedian(ctx->P, ctx->U, ctx->Q, ctx->D);
677:           }
678:           break;

680:         case MatLMVM_Limit_Absolute:
681:           if (ctx->nu) {
682:             VecCopy(ctx->P, ctx->D);
683:             VecShift(ctx->P, -ctx->nu);
684:             VecCopy(ctx->D, ctx->Q);
685:             VecShift(ctx->Q, ctx->nu);
686:             VecMedian(ctx->P, ctx->U, ctx->Q, ctx->P);
687:           }
688:           break;

690:         default:
691:             VecCopy(ctx->U, ctx->D);
692:           break;
693:         }
694:         break;
695:       }
696:       PetscObjectDereference((PetscObject)ctx->Xprev);
697:       PetscObjectDereference((PetscObject)ctx->Gprev);
698:       ctx->Xprev = ctx->S[ctx->lm];
699:       ctx->Gprev = ctx->Y[ctx->lm];
700:       PetscObjectReference((PetscObject)ctx->S[ctx->lm]);
701:       PetscObjectReference((PetscObject)ctx->Y[ctx->lm]);

703:     } else {
704:       ++ctx->nrejects;
705:     }
706:   }

708:   ++ctx->iter;
709:   VecCopy(x, ctx->Xprev);
710:   VecCopy(g, ctx->Gprev);
711:   return(0);
712: }

716: extern PetscErrorCode MatLMVMSetDelta(Mat m, PetscReal d)
717: {
718:   MatLMVMCtx     *ctx;
720:   PetscBool      same;

724:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
725:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
726:   MatShellGetContext(m,(void**)&ctx);
727:   ctx->delta = PetscAbsReal(d);
728:   ctx->delta = PetscMax(ctx->delta_min, ctx->delta);
729:   ctx->delta = PetscMin(ctx->delta_max, ctx->delta);
730:   return(0);
731: }

735: extern PetscErrorCode MatLMVMSetScale(Mat m, Vec s)
736: {
737:   MatLMVMCtx     *ctx;
739:   PetscBool      same;

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

747:   VecDestroy(&ctx->scale);
748:   if (s) {
749:     VecDuplicate(s,&ctx->scale);
750:   }
751:   return(0);
752: }

756: extern PetscErrorCode MatLMVMGetRejects(Mat m, PetscInt *nrejects)
757: {
758:   MatLMVMCtx     *ctx;
760:   PetscBool      same;

764:   PetscObjectTypeCompare((PetscObject)m,MATSHELL,&same);
765:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
766:   MatShellGetContext(m,(void**)&ctx);
767:   *nrejects = ctx->nrejects;
768:   return(0);
769: }

773: extern PetscErrorCode MatLMVMSetH0(Mat m, Mat H0)
774: {
775:   MatLMVMCtx     *ctx;
777:   PetscBool      same;

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

785:   ctx->H0_mat = H0;
786:   PetscObjectReference((PetscObject)ctx->H0_mat);

788:   ctx->useDefaultH0 = PETSC_FALSE;

790:   KSPCreate(PetscObjectComm((PetscObject)H0), &ctx->H0_ksp);
791:   KSPSetOperators(ctx->H0_ksp, H0, H0);
792:   /* its options prefix and setup is handled in TaoSolve_LMVM/TaoSolve_BLMVM */
793:   return(0);
794: }

798: extern PetscErrorCode MatLMVMGetH0(Mat m, Mat *H0)
799: {
800:   MatLMVMCtx     *ctx;
802:   PetscBool      same;

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

809:   MatShellGetContext(m,(void**)&ctx);
810:   *H0  = ctx->H0_mat;
811:   return(0);
812: }

816: extern PetscErrorCode MatLMVMGetH0KSP(Mat m, KSP *H0ksp)
817: {
818:   MatLMVMCtx     *ctx;
820:   PetscBool      same;

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

827:   MatShellGetContext(m,(void**)&ctx);
828:   *H0ksp  = ctx->H0_ksp;
829:   return(0);
830: }

834: extern PetscErrorCode MatLMVMGetX0(Mat m, Vec x)
835: {
837:     return(0);
838: }

842: extern PetscErrorCode MatLMVMSetPrev(Mat M, Vec x, Vec g)
843: {
844:   MatLMVMCtx     *ctx;
846:   PetscBool      same;

851:   PetscObjectTypeCompare((PetscObject)M,MATSHELL,&same);
852:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix M is not type MatLMVM");
853:   MatShellGetContext(M,(void**)&ctx);
854:   if (ctx->nupdates == 0) {
855:     MatLMVMUpdate(M,x,g);
856:   } else {
857:     VecCopy(x,ctx->Xprev);
858:     VecCopy(g,ctx->Gprev);
859:     /*  TODO scaling specific terms */
860:   }
861:   return(0);
862: }

866: extern PetscErrorCode MatLMVMRefine(Mat coarse, Mat op, Vec fineX, Vec fineG)
867: {
869:   PetscBool      same;

876:   PetscObjectTypeCompare((PetscObject)coarse,MATSHELL,&same);
877:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
878:   PetscObjectTypeCompare((PetscObject)op,MATSHELL,&same);
879:   if (!same) SETERRQ(PETSC_COMM_SELF,1,"Matrix m is not type MatLMVM");
880:   return(0);
881: }

885: extern PetscErrorCode MatLMVMAllocateVectors(Mat m, Vec v)
886: {
888:   MatLMVMCtx     *ctx;
889:   PetscBool      same;

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

898:   /*  Perform allocations */
899:   VecDuplicateVecs(v,ctx->lm+1,&ctx->S);
900:   VecDuplicateVecs(v,ctx->lm+1,&ctx->Y);
901:   VecDuplicate(v,&ctx->D);
902:   VecDuplicate(v,&ctx->U);
903:   VecDuplicate(v,&ctx->V);
904:   VecDuplicate(v,&ctx->W);
905:   VecDuplicate(v,&ctx->P);
906:   VecDuplicate(v,&ctx->Q);
907:   VecDuplicate(v,&ctx->H0_norm);
908:   ctx->allocated = PETSC_TRUE;
909:   return(0);
910: }