Actual source code: lmvmmat.c

  1: #include "lmvmmat.h"
  2: #include "tao_general.h"
  3: #include "tao_solver.h"
  4: #include "taovec.h"

  6: #define Scale_None                0
  7: #define Scale_Scalar                1
  8: #define Scale_Broyden                2
  9: #define Scale_Types             3

 11: #define Rescale_None                0
 12: #define Rescale_Scalar                1
 13: #define Rescale_GL                2
 14: #define Rescale_Types                  3

 16: #define Limit_None                0
 17: #define Limit_Average                1
 18: #define Limit_Relative                2
 19: #define Limit_Absolute                3
 20: #define Limit_Types                4

 22: #define TAO_ZER_SAFEGUARD        1e-8
 23: #define TAO_INF_SAFEGUARD        1e+8

 25: static const char *Scale_Table[64] = {
 26:   "none", "scalar", "broyden"
 27: };

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

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

 39: TaoLMVMMat::TaoLMVMMat(TaoVec *tv)
 40: {
 41:   // Set default values
 42:   lm = 5;
 43:   eps = 0.0;

 45:   limitType = Limit_None;
 46:   scaleType = Scale_Broyden;
 47:   rScaleType = Rescale_Scalar;

 49:   s_alpha = 1.0;
 50:   r_alpha = 1.0;
 51:   r_beta = 0.5;
 52:   mu = 1.0;
 53:   nu = 100.0;

 55:   phi = 0.125;                

 57:   scalar_history = 1;
 58:   rescale_history = 1;

 60:   delta_min = 1e-7;
 61:   delta_max = 100;

 63:   // Begin configuration
 64:   TaoOptionsHead("Limited memory matrix options");
 65:   TaoOptionInt("-tao_lmm_vectors", "vectors to use for approximation", "", lm, &lm, 0);
 66:   TaoOptionDouble("-tao_lmm_limit_mu", "mu limiting factor", "", mu, &mu, 0);
 67:   TaoOptionDouble("-tao_lmm_limit_nu", "nu limiting factor", "", nu, &nu, 0);
 68:   TaoOptionDouble("-tao_lmm_broyden_phi", "phi factor for Broyden scaling", "", phi, &phi, 0);
 69:   TaoOptionDouble("-tao_lmm_scalar_alpha", "alpha factor for scalar scaling", "", s_alpha, &s_alpha, 0);
 70:   TaoOptionDouble("-tao_lmm_rescale_alpha", "alpha factor for rescaling diagonal", "", r_alpha, &r_alpha, 0);
 71:   TaoOptionDouble("-tao_lmm_rescale_beta", "beta factor for rescaling diagonal", "", r_beta, &r_beta, 0);
 72:   TaoOptionInt("-tao_lmm_scalar_history", "amount of history for scalar scaling", "", scalar_history, &scalar_history, 0);
 73:   TaoOptionInt("-tao_lmm_rescale_history", "amount of history for rescaling diagonal", "", rescale_history, &rescale_history, 0);
 74:   TaoOptionDouble("-tao_lmm_eps", "rejection tolerance", "", eps, &eps, 0);
 75:   TaoOptionList("-tao_lmm_scale_type", "scale type", "", Scale_Table, Scale_Types, Scale_Table[scaleType], &scaleType, 0);
 76:   TaoOptionList("-tao_lmm_rescale_type", "rescale type", "", Rescale_Table, Rescale_Types, Rescale_Table[rScaleType], &rScaleType, 0);
 77:   TaoOptionList("-tao_lmm_limit_type", "limit type", "", Limit_Table, Limit_Types, Limit_Table[limitType], &limitType, 0);
 78:   TaoOptionDouble("-tao_lmm_delta_min", "minimum delta value", "", delta_min, &delta_min, 0);
 79:   TaoOptionDouble("-tao_lmm_delta_max", "maximum delta value", "", delta_max, &delta_max, 0);
 80:   TaoOptionsTail();

 82:   // Complete configuration
 83:   rescale_history = TaoMin(rescale_history, lm);

 85:   // Perform allocations
 86:   tv->CloneVecs(lm+1, &S);
 87:   tv->CloneVecs(lm+1, &Y);
 88:   tv->Clone(&D);
 89:   tv->Clone(&U);
 90:   tv->Clone(&V);
 91:   tv->Clone(&W);
 92:   tv->Clone(&P);
 93:   tv->Clone(&Q);

 95:   rho = new double[lm+1];
 96:   beta = new double[lm+1];

 98:   yy_history = new double[TaoMax(scalar_history, 1)];
 99:   ys_history = new double[TaoMax(scalar_history, 1)];
100:   ss_history = new double[TaoMax(scalar_history, 1)];

102:   yy_rhistory = new double[TaoMax(rescale_history, 1)];
103:   ys_rhistory = new double[TaoMax(rescale_history, 1)];
104:   ss_rhistory = new double[TaoMax(rescale_history, 1)];

106:   // Finish initializations
107:   lmnow = 0;
108:   iter = 0;
109:   updates = 0;
110:   rejects = 0;
111:   delta = 1.0;

113:   Gprev = 0;
114:   Xprev = 0;

116:   scale = 0;

118:   H0 = 0;
119:   H0default = TAO_TRUE;
120: }

124: TaoLMVMMat::~TaoLMVMMat()
125: {
126:   int info;

128:   info = S[0]->DestroyVecs(lm+1, S);
129:   info = Y[0]->DestroyVecs(lm+1, Y);
130:   info = TaoVecDestroy(D);
131:   info = TaoVecDestroy(U);
132:   info = TaoVecDestroy(V);
133:   info = TaoVecDestroy(W);
134:   info = TaoVecDestroy(P);
135:   info = TaoVecDestroy(Q);

137:   delete[] rho;
138:   delete[] beta;
139:   if (info) rho=0;

141:   delete[] yy_history;
142:   delete[] ys_history;
143:   delete[] ss_history;

145:   delete[] yy_rhistory;
146:   delete[] ys_rhistory;
147:   delete[] ss_rhistory;
148: }

152: int TaoLMVMMat::Refine(TaoLMVMMat *coarse, TaoMat *op, TaoVec *fineX, TaoVec *fineG)
153: {
154:   double rhotemp, rhotol;
155:   double y0temp, s0temp;
156:   int  info;
157:   TaoInt i;

159:   TaoFunctionBegin;

161:   info = Reset(); CHKERRQ(info);

163:   for (i = 0; i < coarse->lmnow; ++i) {
164:     // Refine S[i] and Y[i]
165:     info = op->Multiply(coarse->S[i], S[lmnow]); CHKERRQ(info);
166:     info = op->Multiply(coarse->Y[i], Y[lmnow]); CHKERRQ(info);

168:     // Check to see if refined vectors are fine
169:     info = Y[lmnow]->Dot(S[lmnow], &rhotemp); CHKERRQ(info);
170:     info = Y[lmnow]->Dot(Y[lmnow], &y0temp); CHKERRQ(info);

172:     rhotol = eps * y0temp;
173:     if (rhotemp > rhotol) {
174:       rho[lmnow] = 1.0 / rhotemp;

176:       // Add information to the scaling
177:       switch(scaleType) {
178:       case Scale_None:
179:         break;
180:         
181:       case Scale_Scalar:
182:         // Compute s^T s 
183:         info = S[lmnow]->Dot(S[lmnow], &s0temp); CHKERRQ(info);

185:         // Save information for scalar scaling
186:         yy_history[updates % scalar_history] = y0temp;
187:         ys_history[updates % scalar_history] = rhotemp;
188:         ss_history[updates % scalar_history] = s0temp;

190:         sigma = coarse->sigma;
191:         break;

193:       case Scale_Broyden:
194:         switch(rScaleType) {
195:         case Rescale_None:
196:           break;
197:           
198:         case Rescale_Scalar:
199:         case Rescale_GL:
200:           // Compute s^T s 
201:           info = S[lmnow]->Dot(S[lmnow], &s0temp); CHKERRQ(info);

203:           // Save information for special cases of scalar rescaling
204:           yy_rhistory[updates % rescale_history] = y0temp;
205:           ys_rhistory[updates % rescale_history] = rhotemp;
206:           ss_rhistory[updates % rescale_history] = s0temp;
207:           break;
208:         }

210:         // Obtain diagonal scaling and ensure positive definite
211:         info = op->Multiply(coarse->D, D); CHKERRQ(info);
212:         info = D->AbsoluteValue(); CHKERRQ(info);
213:         break;
214:       }

216:       // Successful update
217:       ++updates;
218:       ++lmnow;
219:     }
220:     else {
221:       ++rejects;
222:     }
223:   }

225:   // Save the number of iterations and previous values for x and g
226:   // Need to use the actual value for the gradient here and not
227:   // the refined version of the coarse gradient in order for the
228:   // Wolfe conditions to guarantee the update is acceptable.
229:   iter = coarse->iter;
230:   info = Xprev->CopyFrom(fineX); CHKERRQ(info);
231:   info = Gprev->CopyFrom(fineG); CHKERRQ(info);
232:   TaoFunctionReturn(0);
233: }

237: int TaoLMVMMat::Reset()
238: {
239:   TaoInt i;

241:   TaoFunctionBegin;
242:   Gprev = Y[lm];
243:   Xprev = S[lm];
244:   for (i = 0; i < lm; ++i) {
245:     rho[i] = 0.0;
246:   }
247:   rho[0] = 1.0;

249:   // Set the scaling and diagonal scaling matrix
250:   switch(scaleType) {
251:   case Scale_None:
252:     sigma = 1.0;
253:     break;

255:   case Scale_Scalar:
256:     sigma = delta;
257:     break;

259:   case Scale_Broyden:
260:     D->SetToConstant(delta);
261:     break;
262:   }

264:   iter = 0;
265:   updates = 0;
266:   lmnow = 0;
267:   TaoFunctionReturn(0);
268: }

272: int TaoLMVMMat::Presolve()
273: {
274:   TaoFunctionBegin;
275:   TaoFunctionReturn(0);
276: }

280: int TaoLMVMMat::Solve(TaoVec *tb, TaoVec *dx, TaoTruth *tt)
281: {
282:   double   sq, yq, dd;
283:   int      info;
284:   TaoInt ll;
285:   TaoTruth scaled;

287:   TaoFunctionBegin;
288:   if (lmnow < 1) {
289:     rho[0] = 1.0;
290:   }

292:   info = dx->CopyFrom(tb); CHKERRQ(info);
293:   for (ll = 0; ll < lmnow; ++ll) {
294:     info = dx->Dot(S[ll], &sq); CHKERRQ(info);
295:     beta[ll] = sq * rho[ll];
296:     info = dx->Axpy(-beta[ll], Y[ll]); CHKERRQ(info);
297:   }

299:   scaled = TAO_FALSE;
300:   if (!scaled && H0) {
301:     info = H0->Solve(dx, U, tt); CHKERRQ(info);
302:     if (*tt) {
303:       info = dx->Dot(U, &dd); CHKERRQ(info);
304:       if ((dd > 0.0) && !TaoInfOrNaN(dd)) {
305:         // Accept Hessian solve
306:         info = dx->CopyFrom(U); CHKERRQ(info);
307:         scaled = TAO_TRUE;
308:       }
309:     }
310:   }

312:   if (!scaled && scale) {
313:     info = U->PointwiseMultiply(dx, scale); CHKERRQ(info);
314:     info = dx->Dot(U, &dd); CHKERRQ(info);
315:     if ((dd > 0.0) && !TaoInfOrNaN(dd)) {
316:       // Accept scaling
317:       info = dx->CopyFrom(U); CHKERRQ(info);
318:       scaled = TAO_TRUE;
319:     }
320:   }
321:   
322:   if (!scaled) {
323:     switch(scaleType) {
324:     case Scale_None:
325:       break;

327:     case Scale_Scalar:
328:       info = dx->Scale(sigma); CHKERRQ(info);
329:       break;
330:   
331:     case Scale_Broyden:
332:       info = dx->PointwiseMultiply(dx, D); CHKERRQ(info);
333:       break;
334:     }
335:   } 

337:   for (ll = lmnow-1; ll >= 0; --ll) {
338:     info = dx->Dot(Y[ll], &yq); CHKERRQ(info);
339:     info = dx->Axpy(beta[ll] - yq * rho[ll], S[ll]); CHKERRQ(info);
340:   }

342:   *tt = TAO_TRUE;
343:   TaoFunctionReturn(0);
344: }

348: int TaoLMVMMat::Update(TaoVec *x, TaoVec *g)
349: {
350:   double rhotemp, rhotol;
351:   double y0temp, s0temp;
352:   double yDy, yDs, sDs;
353:   double sigmanew, denom;
354:   int  info;
355:   TaoInt i;

357:   double yy_sum, ys_sum, ss_sum;

359:   TaoFunctionBegin;

361:   if (0 == iter) {
362:     info = Reset(); CHKERRQ(info);
363:   } 
364:   else {
365:     info = Gprev->Aypx(-1.0, g); CHKERRQ(info);
366:     info = Xprev->Aypx(-1.0, x); CHKERRQ(info);

368:     info = Gprev->Dot(Xprev, &rhotemp); CHKERRQ(info);
369:     info = Gprev->Dot(Gprev, &y0temp); CHKERRQ(info);

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

375:       lmnow = TaoMin(lmnow+1, lm);
376:       for (i = lm-1; i >= 0; --i) {
377:         S[i+1] = S[i];
378:         Y[i+1] = Y[i];
379:         rho[i+1] = rho[i];
380:       }
381:       S[0] = Xprev;
382:       Y[0] = Gprev;
383:       rho[0] = 1.0 / rhotemp;

385:       // Compute the scaling
386:       switch(scaleType) {
387:       case Scale_None:
388:         break;

390:       case Scale_Scalar:
391:         // Compute s^T s 
392:         info = Xprev->Dot(Xprev, &s0temp); CHKERRQ(info);

394:         // Scalar is positive; safeguards are not required.

396:         // Save information for scalar scaling
397:         yy_history[(updates - 1) % scalar_history] = y0temp;
398:         ys_history[(updates - 1) % scalar_history] = rhotemp;
399:         ss_history[(updates - 1) % scalar_history] = s0temp;

401:         // Compute summations for scalar scaling
402:         yy_sum = 0;        // No safeguard required; y^T y > 0
403:         ys_sum = 0;        // No safeguard required; y^T s > 0
404:         ss_sum = 0;        // No safeguard required; s^T s > 0
405:         for (i = 0; i < TaoMin(updates, scalar_history); ++i) {
406:           yy_sum += yy_history[i];
407:           ys_sum += ys_history[i];
408:           ss_sum += ss_history[i];
409:         }

411:         if (0.0 == s_alpha) {
412:           // Safeguard ys_sum 
413:           if (0.0 == ys_sum) {
414:             ys_sum = TAO_ZER_SAFEGUARD;
415:           }

417:           sigmanew = ss_sum / ys_sum;
418:         }
419:         else if (0.5 == s_alpha) {
420:           // Safeguard yy_sum 
421:           if (0.0 == yy_sum) {
422:             yy_sum = TAO_ZER_SAFEGUARD;
423:           }

425:           sigmanew = sqrt(ss_sum / yy_sum);
426:         }
427:         else if (1.0 == s_alpha) {
428:           // Safeguard yy_sum 
429:           if (0.0 == yy_sum) {
430:             yy_sum = TAO_ZER_SAFEGUARD;
431:           }

433:           sigmanew = ys_sum / yy_sum;
434:         }
435:         else {
436:           denom = 2*s_alpha*yy_sum;

438:           // Safeguard denom
439:           if (0.0 == denom) {
440:             denom = TAO_ZER_SAFEGUARD;
441:           }

443:           sigmanew = ((2*s_alpha-1)*ys_sum + 
444:                       sqrt((2*s_alpha-1)*(2*s_alpha-1)*ys_sum*ys_sum - 
445:                            4*s_alpha*(s_alpha-1)*yy_sum*ss_sum)) / denom;
446:         }

448:         switch(limitType) {
449:         case Limit_Average:
450:           if (1.0 == mu) {
451:             sigma = sigmanew;
452:           }
453:           else if (mu) {
454:             sigma = mu * sigmanew + (1.0 - mu) * sigma;
455:           }
456:           break;

458:         case Limit_Relative:
459:           if (mu) {
460:             sigma = TaoMid((1.0 - mu) * sigma, sigmanew, (1.0 + mu) * sigma);
461:           }
462:           break;

464:         case Limit_Absolute:
465:           if (nu) {
466:             sigma = TaoMid(sigma - nu, sigmanew, sigma + nu);
467:           }
468:           break;

470:         default:
471:           sigma = sigmanew;
472:           break;
473:         }
474:         break;

476:       case Scale_Broyden:
477:         // Original version
478:         // Combine DFP and BFGS

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

485:         // P = Q = inv(D);
486:         info = P->CopyFrom(D); CHKERRQ(info);
487:         info = P->Reciprocal(); CHKERRQ(info);
488:         info = Q->CopyFrom(P); CHKERRQ(info);

490:         // V = y*y
491:         info = V->PointwiseMultiply(Gprev, Gprev); CHKERRQ(info);

493:         // W = inv(D)*s
494:         info = W->PointwiseMultiply(Xprev, P); CHKERRQ(info);
495:         info = W->Dot(Xprev, &sDs); CHKERRQ(info);

497:         // Safeguard rhotemp and sDs
498:         if (0.0 == rhotemp) {
499:           rhotemp = TAO_ZER_SAFEGUARD;
500:         }

502:         if (0.0 == sDs) {
503:           sDs = TAO_ZER_SAFEGUARD;
504:         }

506:         if (1.0 != phi) {
507:           // BFGS portion of the update
508:           // U = (inv(D)*s)*(inv(D)*s)
509:           info = U->PointwiseMultiply(W, W); CHKERRQ(info);

511:           // Assemble
512:           info = P->Axpy(1.0 / rhotemp, V); CHKERRQ(info);
513:           info = P->Axpy(-1.0 / sDs, U); CHKERRQ(info);
514:         }

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

521:           // Assemble
522:           info = Q->Axpy(1.0 / rhotemp + sDs / (rhotemp * rhotemp), V); CHKERRQ(info);
523:           info = Q->Axpy(-2.0 / rhotemp, U); CHKERRQ(info);
524:         }

526:         if (0.0 == phi) {
527:           info = U->CopyFrom(P); CHKERRQ(info);
528:         }
529:         else if (1.0 == phi) {
530:           info = U->CopyFrom(Q); CHKERRQ(info);
531:         }
532:         else {
533:           // Broyden update
534:           info = U->Waxpby(1.0 - phi, P, phi, Q); CHKERRQ(info);
535:         }

537:         // Obtain inverse and ensure positive definite
538:         info = U->Reciprocal(); CHKERRQ(info);
539:         info = U->AbsoluteValue(); CHKERRQ(info);

541:         // Checking the diagonal scaling for not a number and infinity
542:         // should not be necessary for the Broyden update
543:         // info = U->Dot(U, &sDs); CHKERRQ(info);
544:         // if (sDs != sDs) {
545:         //   // not a number
546:         //   info = U->SetToConstant(TAO_ZER_SAFEGUARD); CHKERRQ(info);
547:         // }
548:         // else if ((sDs - sDs) != 0.0)) {
549:         //   // infinity
550:         //   info = U->SetToConstant(TAO_INF_SAFEGUARD); CHKERRQ(info);
551:         // }

553:         switch(rScaleType) {
554:         case Rescale_None:
555:           break;

557:         case Rescale_Scalar:
558:         case Rescale_GL:
559:           if (rScaleType == Rescale_GL) {
560:             // Gilbert and Lemarachal use the old diagonal
561:             info = P->CopyFrom(D); CHKERRQ(info);
562:           }
563:           else {
564:             // The default version uses the current diagonal
565:             info = P->CopyFrom(U); CHKERRQ(info);
566:           }

568:           // Compute s^T s 
569:           info = Xprev->Dot(Xprev, &s0temp); CHKERRQ(info);

571:           // Save information for special cases of scalar rescaling
572:           yy_rhistory[(updates - 1) % rescale_history] = y0temp;
573:           ys_rhistory[(updates - 1) % rescale_history] = rhotemp;
574:           ss_rhistory[(updates - 1) % rescale_history] = s0temp;

576:           if (0.0 == r_beta) {
577:             if (1 == TaoMin(updates, rescale_history)) {
578:               // Compute summations for scalar scaling
579:               info = W->PointwiseDivide(S[0], P); CHKERRQ(info);
580:   
581:               info = W->Dot(Y[0], &ys_sum); CHKERRQ(info);
582:               info = W->Dot(W, &ss_sum); CHKERRQ(info);
583:   
584:               yy_sum += yy_rhistory[0];
585:             }
586:             else {
587:               info = Q->CopyFrom(P); CHKERRQ(info);
588:               info = Q->Reciprocal(); CHKERRQ(info);

590:               // Compute summations for scalar scaling
591:               yy_sum = 0;        // No safeguard required
592:               ys_sum = 0;        // No safeguard required
593:               ss_sum = 0;        // No safeguard required
594:               for (i = 0; i < TaoMin(updates, rescale_history); ++i) {
595:                 info = W->PointwiseMultiply(S[i], Q); CHKERRQ(info);
596:                 info = W->Dot(Y[i], &yDs); CHKERRQ(info);
597:                 ys_sum += yDs;

599:                 info = W->Dot(W, &sDs); CHKERRQ(info);
600:                 ss_sum += sDs;
601:   
602:                 yy_sum += yy_rhistory[i];
603:               }
604:             }
605:           }
606:           else if (0.5 == r_beta) {
607:             if (1 == TaoMin(updates, rescale_history)) {
608:               info = V->PointwiseMultiply(Y[0], P); CHKERRQ(info);
609:               info = V->Dot(Y[0], &yy_sum); CHKERRQ(info);

611:               info = W->PointwiseDivide(S[0], P); CHKERRQ(info);
612:               info = W->Dot(S[0], &ss_sum); CHKERRQ(info);

614:               ys_sum = ys_rhistory[0];
615:             }
616:             else {
617:               info = Q->CopyFrom(P); CHKERRQ(info);
618:               info = Q->Reciprocal(); CHKERRQ(info);

620:               // Compute summations for scalar scaling
621:               yy_sum = 0;        // No safeguard required
622:               ys_sum = 0;        // No safeguard required
623:               ss_sum = 0;        // No safeguard required
624:               for (i = 0; i < TaoMin(updates, rescale_history); ++i) {
625:                 info = V->PointwiseMultiply(Y[i], P); CHKERRQ(info);
626:                 info = V->Dot(Y[i], &yDy); CHKERRQ(info);
627:                 yy_sum += yDy;

629:                 info = W->PointwiseMultiply(S[i], Q); CHKERRQ(info);
630:                 info = W->Dot(S[i], &sDs); CHKERRQ(info);
631:                 ss_sum += sDs;

633:                 ys_sum += ys_rhistory[i];
634:               }
635:             }
636:           }
637:           else if (1.0 == r_beta) {
638:             // Compute summations for scalar scaling
639:             yy_sum = 0;        // No safeguard required
640:             ys_sum = 0;        // No safeguard required
641:             ss_sum = 0;        // No safeguard required
642:             for (i = 0; i < TaoMin(updates, rescale_history); ++i) {
643:               info = V->PointwiseMultiply(Y[i], P); CHKERRQ(info);
644:               info = V->Dot(S[i], &yDs); CHKERRQ(info);
645:               ys_sum += yDs;

647:               info = V->Dot(V, &yDy); CHKERRQ(info);
648:               yy_sum += yDy;

650:               ss_sum += ss_rhistory[i];
651:             }
652:           }
653:           else {
654:             info = Q->CopyFrom(P); CHKERRQ(info);

656:             info = P->Pow(r_beta); CHKERRQ(info);
657:             info = Q->PointwiseDivide(P, Q); CHKERRQ(info);

659:             // Compute summations for scalar scaling
660:             yy_sum = 0;        // No safeguard required
661:             ys_sum = 0;        // No safeguard required
662:             ss_sum = 0;        // No safeguard required
663:             for (i = 0; i < TaoMin(updates, rescale_history); ++i) {
664:               info = V->PointwiseMultiply(P, Y[i]); CHKERRQ(info);
665:               info = W->PointwiseMultiply(Q, S[i]); CHKERRQ(info);

667:               info = V->Dot(V, &yDy); CHKERRQ(info);
668:               info = V->Dot(W, &yDs); CHKERRQ(info);
669:               info = W->Dot(W, &sDs); CHKERRQ(info);

671:               yy_sum += yDy;
672:               ys_sum += yDs;
673:               ss_sum += sDs;
674:             }
675:           }

677:           if (0.0 == r_alpha) {
678:             // Safeguard ys_sum 
679:             if (0.0 == ys_sum) {
680:               ys_sum = TAO_ZER_SAFEGUARD;
681:             }

683:             sigmanew = ss_sum / ys_sum;
684:           }
685:           else if (0.5 == r_alpha) {
686:             // Safeguard yy_sum 
687:             if (0.0 == yy_sum) {
688:               yy_sum = TAO_ZER_SAFEGUARD;
689:             }

691:             sigmanew = sqrt(ss_sum / yy_sum);
692:           }
693:           else if (1.0 == r_alpha) {
694:             // Safeguard yy_sum 
695:             if (0.0 == yy_sum) {
696:               yy_sum = TAO_ZER_SAFEGUARD;
697:             }

699:             sigmanew = ys_sum / yy_sum;
700:           }
701:           else {
702:             denom = 2*r_alpha*yy_sum;

704:             // Safeguard denom
705:             if (0.0 == denom) {
706:               denom = TAO_ZER_SAFEGUARD;
707:             }

709:             sigmanew = ((2*r_alpha-1)*ys_sum +
710:                         sqrt((2*r_alpha-1)*(2*r_alpha-1)*ys_sum*ys_sum -
711:                              4*r_alpha*(r_alpha-1)*yy_sum*ss_sum)) / denom;
712:           }

714:           // If Q has small values, then Q^(r_beta - 1)
715:           // can have very large values.  Hence, ys_sum
716:           // and ss_sum can be infinity.  In this case,
717:           // sigmanew can either be not-a-number or infinity.

719:           if (TaoInfOrNaN(sigmanew)) {
720:             // sigmanew is not-a-number; skip rescaling
721:           }
722:           else if (!sigmanew) {
723:             // sigmanew is zero; this is a bad case; skip rescaling
724:           }
725:           else {
726:             // sigmanew is positive
727:             info = U->Scale(sigmanew); CHKERRQ(info);
728:           }
729:           break;
730:         }

732:         // Modify for previous information
733:         switch(limitType) {
734:         case Limit_Average:
735:           if (1.0 == mu) {
736:             info = D->CopyFrom(U); CHKERRQ(info);
737:           }
738:           else if (mu) {
739:             info = D->Axpby(mu, U, 1.0 - mu); CHKERRQ(info);
740:           }
741:           break;
742:  
743:         case Limit_Relative:
744:           if (mu) {
745:             info = P->ScaleCopyFrom(1.0 - mu, D); CHKERRQ(info);
746:             info = Q->ScaleCopyFrom(1.0 + mu, D); CHKERRQ(info);
747:             info = D->Median(P, U, Q); CHKERRQ(info);
748:           }
749:           break;

751:         case Limit_Absolute:
752:           if (nu) {
753:             info = P->CopyFrom(D); CHKERRQ(info);
754:             info = P->AddConstant(-nu); CHKERRQ(info);
755:                info = Q->CopyFrom(D); CHKERRQ(info);
756:             info = Q->AddConstant(nu); CHKERRQ(info);
757:             info = D->Median(P, U, Q); CHKERRQ(info);
758:           }
759:           break;

761:         default:
762:           info = D->CopyFrom(U); CHKERRQ(info);
763:           break;
764:         } 
765:         break;
766:       }

768:       Xprev = S[lm]; Gprev = Y[lm];
769:     } 
770:     else { 
771:       ++rejects;
772:     }
773:   }
774:   
775:   ++iter;
776:   info = Xprev->CopyFrom(x); CHKERRQ(info);
777:   info = Gprev->CopyFrom(g); CHKERRQ(info);
778:   TaoFunctionReturn(0);
779: }

783: int TaoLMVMMat::View()
784: {
785:   TaoFunctionBegin;
786:   TaoFunctionReturn(0);
787: }

791: int TaoLMVMMat::SetDelta(double d)
792: {
793:   TaoFunctionBegin;
794:   delta = TaoAbsDouble(d);
795:   delta = TaoMax(delta_min, delta);
796:   delta = TaoMin(delta_max, delta);
797:   TaoFunctionReturn(0);
798: }

802: int TaoLMVMMat::SetScale(TaoVec *s)
803: {
804:   TaoFunctionBegin;
805:   scale = s;
806:   TaoFunctionReturn(0);
807: }

811: int TaoLMVMMat::SetH0(TaoMat *HH0)
812: {
813:   TaoFunctionBegin;
814:   if (H0) { 
815:     H0 = 0;
816:   }
817:   H0default = TAO_TRUE;

819:   if (HH0) {
820:     H0 = HH0;
821:     H0default = TAO_FALSE;
822:   }
823:   TaoFunctionReturn(0);
824: }

828: int TaoLMVMMat::GetX0(TaoVec *x)
829: {
830:   int i,info;

832:   TaoFunctionBegin;
833:   info = x->CopyFrom(Xprev); CHKERRQ(info);
834:   for (i = 0; i < lmnow; ++i) {
835:     info = x->Axpy(-1.0, S[i]); CHKERRQ(info);
836:   }
837:   TaoFunctionReturn(0);
838: }

842: int TaoLMVMMat::InitialApproximation(TaoVec *x)
843: {
844:   TaoFunctionBegin;
845:   TaoFunctionReturn(0);
846: }