Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   1:  #region Translated by Jose Antonio De Santiago-Castillo.
   2:   
   3:  //Translated by Jose Antonio De Santiago-Castillo. 
   4:  //E-mail:JAntonioDeSantiago@gmail.com
   5:  //Web: www.DotNumerics.com
   6:  //
   7:  //Fortran to C# Translation.
   8:  //Translated by:
   9:  //F2CSharp Version 0.71 (November 10, 2009)
  10:  //Code Optimizations: None
  11:  //
  12:  #endregion
  13:   
  14:  using System;
  15:  using DotNumerics.FortranLibrary;
  16:   
  17:  namespace DotNumerics.CSLapack
  18:  {
  19:      /// <summary>
  20:      /// -- LAPACK auxiliary routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DLAIC1 applies one step of incremental condition estimation in
  27:      /// its simplest version:
  28:      /// 
  29:      /// Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
  30:      /// lower triangular matrix L, such that
  31:      /// twonorm(L*x) = sest
  32:      /// Then DLAIC1 computes sestpr, s, c such that
  33:      /// the vector
  34:      /// [ s*x ]
  35:      /// xhat = [  c  ]
  36:      /// is an approximate singular vector of
  37:      /// [ L     0  ]
  38:      /// Lhat = [ w' gamma ]
  39:      /// in the sense that
  40:      /// twonorm(Lhat*xhat) = sestpr.
  41:      /// 
  42:      /// Depending on JOB, an estimate for the largest or smallest singular
  43:      /// value is computed.
  44:      /// 
  45:      /// Note that [s c]' and sestpr**2 is an eigenpair of the system
  46:      /// 
  47:      /// diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
  48:      /// [ gamma ]
  49:      /// 
  50:      /// where  alpha =  x'*w.
  51:      /// 
  52:      ///</summary>
  53:      public class DLAIC1
  54:      {
  55:      
  56:   
  57:          #region Dependencies
  58:          
  59:          DDOT _ddot; DLAMCH _dlamch; 
  60:   
  61:          #endregion
  62:   
  63:   
  64:          #region Fields
  65:          
  66:          const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; const double HALF = 0.5E0; 
  67:          const double FOUR = 4.0E0;double ABSALP = 0; double ABSEST = 0; double ABSGAM = 0; double ALPHA = 0; double B = 0; 
  68:          double COSINE = 0;double EPS = 0; double NORMA = 0; double S1 = 0; double S2 = 0; double SINE = 0; double T = 0; 
  69:          double TEST = 0;double TMP = 0; double ZETA1 = 0; double ZETA2 = 0; 
  70:   
  71:          #endregion
  72:   
  73:          public DLAIC1(DDOT ddot, DLAMCH dlamch)
  74:          {
  75:      
  76:   
  77:              #region Set Dependencies
  78:              
  79:              this._ddot = ddot; this._dlamch = dlamch; 
  80:   
  81:              #endregion
  82:   
  83:          }
  84:      
  85:          public DLAIC1()
  86:          {
  87:      
  88:   
  89:              #region Dependencies (Initialization)
  90:              
  91:              DDOT ddot = new DDOT();
  92:              LSAME lsame = new LSAME();
  93:              DLAMC3 dlamc3 = new DLAMC3();
  94:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  95:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  96:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  97:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  98:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  99:   
 100:              #endregion
 101:   
 102:   
 103:              #region Set Dependencies
 104:              
 105:              this._ddot = ddot; this._dlamch = dlamch; 
 106:   
 107:              #endregion
 108:   
 109:          }
 110:          /// <summary>
 111:          /// Purpose
 112:          /// =======
 113:          /// 
 114:          /// DLAIC1 applies one step of incremental condition estimation in
 115:          /// its simplest version:
 116:          /// 
 117:          /// Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
 118:          /// lower triangular matrix L, such that
 119:          /// twonorm(L*x) = sest
 120:          /// Then DLAIC1 computes sestpr, s, c such that
 121:          /// the vector
 122:          /// [ s*x ]
 123:          /// xhat = [  c  ]
 124:          /// is an approximate singular vector of
 125:          /// [ L     0  ]
 126:          /// Lhat = [ w' gamma ]
 127:          /// in the sense that
 128:          /// twonorm(Lhat*xhat) = sestpr.
 129:          /// 
 130:          /// Depending on JOB, an estimate for the largest or smallest singular
 131:          /// value is computed.
 132:          /// 
 133:          /// Note that [s c]' and sestpr**2 is an eigenpair of the system
 134:          /// 
 135:          /// diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
 136:          /// [ gamma ]
 137:          /// 
 138:          /// where  alpha =  x'*w.
 139:          /// 
 140:          ///</summary>
 141:          /// <param name="JOB">
 142:          /// (input) INTEGER
 143:          /// = 1: an estimate for the largest singular value is computed.
 144:          /// = 2: an estimate for the smallest singular value is computed.
 145:          ///</param>
 146:          /// <param name="J">
 147:          /// (input) INTEGER
 148:          /// Length of X and W
 149:          ///</param>
 150:          /// <param name="X">
 151:          /// (input) DOUBLE PRECISION array, dimension (J)
 152:          /// The j-vector x.
 153:          ///</param>
 154:          /// <param name="SEST">
 155:          /// (input) DOUBLE PRECISION
 156:          /// Estimated singular value of j by j matrix L
 157:          ///</param>
 158:          /// <param name="W">
 159:          /// (input) DOUBLE PRECISION array, dimension (J)
 160:          /// The j-vector w.
 161:          ///</param>
 162:          /// <param name="GAMMA">
 163:          /// (input) DOUBLE PRECISION
 164:          /// The diagonal element gamma.
 165:          ///</param>
 166:          /// <param name="SESTPR">
 167:          /// (output) DOUBLE PRECISION
 168:          /// Estimated singular value of (j+1) by (j+1) matrix Lhat.
 169:          ///</param>
 170:          /// <param name="S">
 171:          /// (output) DOUBLE PRECISION
 172:          /// Sine needed in forming xhat.
 173:          ///</param>
 174:          /// <param name="C">
 175:          /// (output) DOUBLE PRECISION
 176:          /// Cosine needed in forming xhat.
 177:          ///</param>
 178:          public void Run(int JOB, int J, double[] X, int offset_x, double SEST, double[] W, int offset_w, double GAMMA
 179:                           , ref double SESTPR, ref double S, ref double C)
 180:          {
 181:   
 182:              #region Array Index Correction
 183:              
 184:               int o_x = -1 + offset_x;  int o_w = -1 + offset_w; 
 185:   
 186:              #endregion
 187:   
 188:   
 189:              #region Prolog
 190:              
 191:              // *
 192:              // *  -- LAPACK auxiliary routine (version 3.1) --
 193:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 194:              // *     November 2006
 195:              // *
 196:              // *     .. Scalar Arguments ..
 197:              // *     ..
 198:              // *     .. Array Arguments ..
 199:              // *     ..
 200:              // *
 201:              // *  Purpose
 202:              // *  =======
 203:              // *
 204:              // *  DLAIC1 applies one step of incremental condition estimation in
 205:              // *  its simplest version:
 206:              // *
 207:              // *  Let x, twonorm(x) = 1, be an approximate singular vector of an j-by-j
 208:              // *  lower triangular matrix L, such that
 209:              // *           twonorm(L*x) = sest
 210:              // *  Then DLAIC1 computes sestpr, s, c such that
 211:              // *  the vector
 212:              // *                  [ s*x ]
 213:              // *           xhat = [  c  ]
 214:              // *  is an approximate singular vector of
 215:              // *                  [ L     0  ]
 216:              // *           Lhat = [ w' gamma ]
 217:              // *  in the sense that
 218:              // *           twonorm(Lhat*xhat) = sestpr.
 219:              // *
 220:              // *  Depending on JOB, an estimate for the largest or smallest singular
 221:              // *  value is computed.
 222:              // *
 223:              // *  Note that [s c]' and sestpr**2 is an eigenpair of the system
 224:              // *
 225:              // *      diag(sest*sest, 0) + [alpha  gamma] * [ alpha ]
 226:              // *                                            [ gamma ]
 227:              // *
 228:              // *  where  alpha =  x'*w.
 229:              // *
 230:              // *  Arguments
 231:              // *  =========
 232:              // *
 233:              // *  JOB     (input) INTEGER
 234:              // *          = 1: an estimate for the largest singular value is computed.
 235:              // *          = 2: an estimate for the smallest singular value is computed.
 236:              // *
 237:              // *  J       (input) INTEGER
 238:              // *          Length of X and W
 239:              // *
 240:              // *  X       (input) DOUBLE PRECISION array, dimension (J)
 241:              // *          The j-vector x.
 242:              // *
 243:              // *  SEST    (input) DOUBLE PRECISION
 244:              // *          Estimated singular value of j by j matrix L
 245:              // *
 246:              // *  W       (input) DOUBLE PRECISION array, dimension (J)
 247:              // *          The j-vector w.
 248:              // *
 249:              // *  GAMMA   (input) DOUBLE PRECISION
 250:              // *          The diagonal element gamma.
 251:              // *
 252:              // *  SESTPR  (output) DOUBLE PRECISION
 253:              // *          Estimated singular value of (j+1) by (j+1) matrix Lhat.
 254:              // *
 255:              // *  S       (output) DOUBLE PRECISION
 256:              // *          Sine needed in forming xhat.
 257:              // *
 258:              // *  C       (output) DOUBLE PRECISION
 259:              // *          Cosine needed in forming xhat.
 260:              // *
 261:              // *  =====================================================================
 262:              // *
 263:              // *     .. Parameters ..
 264:              // *     ..
 265:              // *     .. Local Scalars ..
 266:              // *     ..
 267:              // *     .. Intrinsic Functions ..
 268:              //      INTRINSIC          ABS, MAX, SIGN, SQRT;
 269:              // *     ..
 270:              // *     .. External Functions ..
 271:              // *     ..
 272:              // *     .. Executable Statements ..
 273:              // *
 274:   
 275:              #endregion
 276:   
 277:   
 278:              #region Body
 279:              
 280:              EPS = this._dlamch.Run("Epsilon");
 281:              ALPHA = this._ddot.Run(J, X, offset_x, 1, W, offset_w, 1);
 282:              // *
 283:              ABSALP = Math.Abs(ALPHA);
 284:              ABSGAM = Math.Abs(GAMMA);
 285:              ABSEST = Math.Abs(SEST);
 286:              // *
 287:              if (JOB == 1)
 288:              {
 289:                  // *
 290:                  // *        Estimating largest singular value
 291:                  // *
 292:                  // *        special cases
 293:                  // *
 294:                  if (SEST == ZERO)
 295:                  {
 296:                      S1 = Math.Max(ABSGAM, ABSALP);
 297:                      if (S1 == ZERO)
 298:                      {
 299:                          S = ZERO;
 300:                          C = ONE;
 301:                          SESTPR = ZERO;
 302:                      }
 303:                      else
 304:                      {
 305:                          S = ALPHA / S1;
 306:                          C = GAMMA / S1;
 307:                          TMP = Math.Sqrt(S * S + C * C);
 308:                          S = S / TMP;
 309:                          C = C / TMP;
 310:                          SESTPR = S1 * TMP;
 311:                      }
 312:                      return;
 313:                  }
 314:                  else
 315:                  {
 316:                      if (ABSGAM <= EPS * ABSEST)
 317:                      {
 318:                          S = ONE;
 319:                          C = ZERO;
 320:                          TMP = Math.Max(ABSEST, ABSALP);
 321:                          S1 = ABSEST / TMP;
 322:                          S2 = ABSALP / TMP;
 323:                          SESTPR = TMP * Math.Sqrt(S1 * S1 + S2 * S2);
 324:                          return;
 325:                      }
 326:                      else
 327:                      {
 328:                          if (ABSALP <= EPS * ABSEST)
 329:                          {
 330:                              S1 = ABSGAM;
 331:                              S2 = ABSEST;
 332:                              if (S1 <= S2)
 333:                              {
 334:                                  S = ONE;
 335:                                  C = ZERO;
 336:                                  SESTPR = S2;
 337:                              }
 338:                              else
 339:                              {
 340:                                  S = ZERO;
 341:                                  C = ONE;
 342:                                  SESTPR = S1;
 343:                              }
 344:                              return;
 345:                          }
 346:                          else
 347:                          {
 348:                              if (ABSEST <= EPS * ABSALP || ABSEST <= EPS * ABSGAM)
 349:                              {
 350:                                  S1 = ABSGAM;
 351:                                  S2 = ABSALP;
 352:                                  if (S1 <= S2)
 353:                                  {
 354:                                      TMP = S1 / S2;
 355:                                      S = Math.Sqrt(ONE + TMP * TMP);
 356:                                      SESTPR = S2 * S;
 357:                                      C = (GAMMA / S2) / S;
 358:                                      S = FortranLib.Sign(ONE,ALPHA) / S;
 359:                                  }
 360:                                  else
 361:                                  {
 362:                                      TMP = S2 / S1;
 363:                                      C = Math.Sqrt(ONE + TMP * TMP);
 364:                                      SESTPR = S1 * C;
 365:                                      S = (ALPHA / S1) / C;
 366:                                      C = FortranLib.Sign(ONE,GAMMA) / C;
 367:                                  }
 368:                                  return;
 369:                              }
 370:                              else
 371:                              {
 372:                                  // *
 373:                                  // *           normal case
 374:                                  // *
 375:                                  ZETA1 = ALPHA / ABSEST;
 376:                                  ZETA2 = GAMMA / ABSEST;
 377:                                  // *
 378:                                  B = (ONE - ZETA1 * ZETA1 - ZETA2 * ZETA2) * HALF;
 379:                                  C = ZETA1 * ZETA1;
 380:                                  if (B > ZERO)
 381:                                  {
 382:                                      T = C / (B + Math.Sqrt(B * B + C));
 383:                                  }
 384:                                  else
 385:                                  {
 386:                                      T = Math.Sqrt(B * B + C) - B;
 387:                                  }
 388:                                  // *
 389:                                  SINE =  - ZETA1 / T;
 390:                                  COSINE =  - ZETA2 / (ONE + T);
 391:                                  TMP = Math.Sqrt(SINE * SINE + COSINE * COSINE);
 392:                                  S = SINE / TMP;
 393:                                  C = COSINE / TMP;
 394:                                  SESTPR = Math.Sqrt(T + ONE) * ABSEST;
 395:                                  return;
 396:                              }
 397:                          }
 398:                      }
 399:                  }
 400:                  // *
 401:              }
 402:              else
 403:              {
 404:                  if (JOB == 2)
 405:                  {
 406:                      // *
 407:                      // *        Estimating smallest singular value
 408:                      // *
 409:                      // *        special cases
 410:                      // *
 411:                      if (SEST == ZERO)
 412:                      {
 413:                          SESTPR = ZERO;
 414:                          if (Math.Max(ABSGAM, ABSALP) == ZERO)
 415:                          {
 416:                              SINE = ONE;
 417:                              COSINE = ZERO;
 418:                          }
 419:                          else
 420:                          {
 421:                              SINE =  - GAMMA;
 422:                              COSINE = ALPHA;
 423:                          }
 424:                          S1 = Math.Max(Math.Abs(SINE), Math.Abs(COSINE));
 425:                          S = SINE / S1;
 426:                          C = COSINE / S1;
 427:                          TMP = Math.Sqrt(S * S + C * C);
 428:                          S = S / TMP;
 429:                          C = C / TMP;
 430:                          return;
 431:                      }
 432:                      else
 433:                      {
 434:                          if (ABSGAM <= EPS * ABSEST)
 435:                          {
 436:                              S = ZERO;
 437:                              C = ONE;
 438:                              SESTPR = ABSGAM;
 439:                              return;
 440:                          }
 441:                          else
 442:                          {
 443:                              if (ABSALP <= EPS * ABSEST)
 444:                              {
 445:                                  S1 = ABSGAM;
 446:                                  S2 = ABSEST;
 447:                                  if (S1 <= S2)
 448:                                  {
 449:                                      S = ZERO;
 450:                                      C = ONE;
 451:                                      SESTPR = S1;
 452:                                  }
 453:                                  else
 454:                                  {
 455:                                      S = ONE;
 456:                                      C = ZERO;
 457:                                      SESTPR = S2;
 458:                                  }
 459:                                  return;
 460:                              }
 461:                              else
 462:                              {
 463:                                  if (ABSEST <= EPS * ABSALP || ABSEST <= EPS * ABSGAM)
 464:                                  {
 465:                                      S1 = ABSGAM;
 466:                                      S2 = ABSALP;
 467:                                      if (S1 <= S2)
 468:                                      {
 469:                                          TMP = S1 / S2;
 470:                                          C = Math.Sqrt(ONE + TMP * TMP);
 471:                                          SESTPR = ABSEST * (TMP / C);
 472:                                          S =  - (GAMMA / S2) / C;
 473:                                          C = FortranLib.Sign(ONE,ALPHA) / C;
 474:                                      }
 475:                                      else
 476:                                      {
 477:                                          TMP = S2 / S1;
 478:                                          S = Math.Sqrt(ONE + TMP * TMP);
 479:                                          SESTPR = ABSEST / S;
 480:                                          C = (ALPHA / S1) / S;
 481:                                          S =  - FortranLib.Sign(ONE,GAMMA) / S;
 482:                                      }
 483:                                      return;
 484:                                  }
 485:                                  else
 486:                                  {
 487:                                      // *
 488:                                      // *           normal case
 489:                                      // *
 490:                                      ZETA1 = ALPHA / ABSEST;
 491:                                      ZETA2 = GAMMA / ABSEST;
 492:                                      // *
 493:                                      NORMA = Math.Max(ONE + ZETA1 * ZETA1 + Math.Abs(ZETA1 * ZETA2), Math.Abs(ZETA1 * ZETA2) + ZETA2 * ZETA2);
 494:                                      // *
 495:                                      // *           See if root is closer to zero or to ONE
 496:                                      // *
 497:                                      TEST = ONE + TWO * (ZETA1 - ZETA2) * (ZETA1 + ZETA2);
 498:                                      if (TEST >= ZERO)
 499:                                      {
 500:                                          // *
 501:                                          // *              root is close to zero, compute directly
 502:                                          // *
 503:                                          B = (ZETA1 * ZETA1 + ZETA2 * ZETA2 + ONE) * HALF;
 504:                                          C = ZETA2 * ZETA2;
 505:                                          T = C / (B + Math.Sqrt(Math.Abs(B * B - C)));
 506:                                          SINE = ZETA1 / (ONE - T);
 507:                                          COSINE =  - ZETA2 / T;
 508:                                          SESTPR = Math.Sqrt(T + FOUR * EPS * EPS * NORMA) * ABSEST;
 509:                                      }
 510:                                      else
 511:                                      {
 512:                                          // *
 513:                                          // *              root is closer to ONE, shift by that amount
 514:                                          // *
 515:                                          B = (ZETA2 * ZETA2 + ZETA1 * ZETA1 - ONE) * HALF;
 516:                                          C = ZETA1 * ZETA1;
 517:                                          if (B >= ZERO)
 518:                                          {
 519:                                              T =  - C / (B + Math.Sqrt(B * B + C));
 520:                                          }
 521:                                          else
 522:                                          {
 523:                                              T = B - Math.Sqrt(B * B + C);
 524:                                          }
 525:                                          SINE =  - ZETA1 / T;
 526:                                          COSINE =  - ZETA2 / (ONE + T);
 527:                                          SESTPR = Math.Sqrt(ONE + T + FOUR * EPS * EPS * NORMA) * ABSEST;
 528:                                      }
 529:                                      TMP = Math.Sqrt(SINE * SINE + COSINE * COSINE);
 530:                                      S = SINE / TMP;
 531:                                      C = COSINE / TMP;
 532:                                      return;
 533:                                      // *
 534:                                  }
 535:                              }
 536:                          }
 537:                      }
 538:                  }
 539:              }
 540:              return;
 541:              // *
 542:              // *     End of DLAIC1
 543:              // *
 544:   
 545:              #endregion
 546:   
 547:          }
 548:      }
 549:  }