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 routine (version 3.1.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// February 2007
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DLAED6 computes the positive or negative root (closest to the origin)
  27:      /// of
  28:      /// z(1)        z(2)        z(3)
  29:      /// f(x) =   rho + --------- + ---------- + ---------
  30:      /// d(1)-x      d(2)-x      d(3)-x
  31:      /// 
  32:      /// It is assumed that
  33:      /// 
  34:      /// if ORGATI = .true. the root is between d(2) and d(3);
  35:      /// otherwise it is between d(1) and d(2)
  36:      /// 
  37:      /// This routine will be called by DLAED4 when necessary. In most cases,
  38:      /// the root sought is the smallest in magnitude, though it might not be
  39:      /// in some extremely rare situations.
  40:      /// 
  41:      ///</summary>
  42:      public class DLAED6
  43:      {
  44:      
  45:   
  46:          #region Dependencies
  47:          
  48:          DLAMCH _dlamch; 
  49:   
  50:          #endregion
  51:   
  52:   
  53:          #region Fields
  54:          
  55:          const int MAXIT = 40; const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; 
  56:          const double THREE = 3.0E0;const double FOUR = 4.0E0; const double EIGHT = 8.0E0; 
  57:          double[] DSCALE = new double[3]; int o_dscale = -1;
  58:          double[] ZSCALE = new double[3]; int o_zscale = -1;bool SCALE = false; int I = 0; int ITER = 0; 
  59:          int NITER = 0;double A = 0; double B = 0; double BASE = 0; double C = 0; double DDF = 0; double DF = 0; double EPS = 0; 
  60:          double ERRETM = 0;double ETA = 0; double F = 0; double FC = 0; double SCLFAC = 0; double SCLINV = 0; double SMALL1 = 0; 
  61:          double SMALL2 = 0;double SMINV1 = 0; double SMINV2 = 0; double TEMP = 0; double TEMP1 = 0; double TEMP2 = 0; 
  62:          double TEMP3 = 0;double TEMP4 = 0; double LBD = 0; double UBD = 0; 
  63:   
  64:          #endregion
  65:   
  66:          public DLAED6(DLAMCH dlamch)
  67:          {
  68:      
  69:   
  70:              #region Set Dependencies
  71:              
  72:              this._dlamch = dlamch; 
  73:   
  74:              #endregion
  75:   
  76:          }
  77:      
  78:          public DLAED6()
  79:          {
  80:      
  81:   
  82:              #region Dependencies (Initialization)
  83:              
  84:              LSAME lsame = new LSAME();
  85:              DLAMC3 dlamc3 = new DLAMC3();
  86:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  87:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  88:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  89:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  90:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  91:   
  92:              #endregion
  93:   
  94:   
  95:              #region Set Dependencies
  96:              
  97:              this._dlamch = dlamch; 
  98:   
  99:              #endregion
 100:   
 101:          }
 102:          /// <summary>
 103:          /// Purpose
 104:          /// =======
 105:          /// 
 106:          /// DLAED6 computes the positive or negative root (closest to the origin)
 107:          /// of
 108:          /// z(1)        z(2)        z(3)
 109:          /// f(x) =   rho + --------- + ---------- + ---------
 110:          /// d(1)-x      d(2)-x      d(3)-x
 111:          /// 
 112:          /// It is assumed that
 113:          /// 
 114:          /// if ORGATI = .true. the root is between d(2) and d(3);
 115:          /// otherwise it is between d(1) and d(2)
 116:          /// 
 117:          /// This routine will be called by DLAED4 when necessary. In most cases,
 118:          /// the root sought is the smallest in magnitude, though it might not be
 119:          /// in some extremely rare situations.
 120:          /// 
 121:          ///</summary>
 122:          /// <param name="KNITER">
 123:          /// (input) INTEGER
 124:          /// Refer to DLAED4 for its significance.
 125:          ///</param>
 126:          /// <param name="ORGATI">
 127:          /// (input) LOGICAL
 128:          /// If ORGATI is true, the needed root is between d(2) and
 129:          /// d(3); otherwise it is between d(1) and d(2).  See
 130:          /// DLAED4 for further details.
 131:          ///</param>
 132:          /// <param name="RHO">
 133:          /// (input) DOUBLE PRECISION
 134:          /// Refer to the equation f(x) above.
 135:          ///</param>
 136:          /// <param name="D">
 137:          /// (input) DOUBLE PRECISION array, dimension (3)
 138:          /// D satisfies d(1) .LT. d(2) .LT. d(3).
 139:          ///</param>
 140:          /// <param name="Z">
 141:          /// (input) DOUBLE PRECISION array, dimension (3)
 142:          /// Each of the elements in z must be positive.
 143:          ///</param>
 144:          /// <param name="FINIT">
 145:          /// (input) DOUBLE PRECISION
 146:          /// The value of f at 0. It is more accurate than the one
 147:          /// evaluated inside this routine (if someone wants to do
 148:          /// so).
 149:          ///</param>
 150:          /// <param name="TAU">
 151:          /// (output) DOUBLE PRECISION
 152:          /// The root of the equation f(x).
 153:          ///</param>
 154:          /// <param name="INFO">
 155:          /// (output) INTEGER
 156:          /// = 0: successful exit
 157:          /// .GT. 0: if INFO = 1, failure to converge
 158:          ///</param>
 159:          public void Run(int KNITER, bool ORGATI, double RHO, double[] D, int offset_d, double[] Z, int offset_z, double FINIT
 160:                           , ref double TAU, ref int INFO)
 161:          {
 162:   
 163:              #region Array Index Correction
 164:              
 165:               int o_d = -1 + offset_d;  int o_z = -1 + offset_z; 
 166:   
 167:              #endregion
 168:   
 169:   
 170:              #region Prolog
 171:              
 172:              // *
 173:              // *  -- LAPACK routine (version 3.1.1) --
 174:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 175:              // *     February 2007
 176:              // *
 177:              // *     .. Scalar Arguments ..
 178:              // *     ..
 179:              // *     .. Array Arguments ..
 180:              // *     ..
 181:              // *
 182:              // *  Purpose
 183:              // *  =======
 184:              // *
 185:              // *  DLAED6 computes the positive or negative root (closest to the origin)
 186:              // *  of
 187:              // *                   z(1)        z(2)        z(3)
 188:              // *  f(x) =   rho + --------- + ---------- + ---------
 189:              // *                  d(1)-x      d(2)-x      d(3)-x
 190:              // *
 191:              // *  It is assumed that
 192:              // *
 193:              // *        if ORGATI = .true. the root is between d(2) and d(3);
 194:              // *        otherwise it is between d(1) and d(2)
 195:              // *
 196:              // *  This routine will be called by DLAED4 when necessary. In most cases,
 197:              // *  the root sought is the smallest in magnitude, though it might not be
 198:              // *  in some extremely rare situations.
 199:              // *
 200:              // *  Arguments
 201:              // *  =========
 202:              // *
 203:              // *  KNITER       (input) INTEGER
 204:              // *               Refer to DLAED4 for its significance.
 205:              // *
 206:              // *  ORGATI       (input) LOGICAL
 207:              // *               If ORGATI is true, the needed root is between d(2) and
 208:              // *               d(3); otherwise it is between d(1) and d(2).  See
 209:              // *               DLAED4 for further details.
 210:              // *
 211:              // *  RHO          (input) DOUBLE PRECISION
 212:              // *               Refer to the equation f(x) above.
 213:              // *
 214:              // *  D            (input) DOUBLE PRECISION array, dimension (3)
 215:              // *               D satisfies d(1) < d(2) < d(3).
 216:              // *
 217:              // *  Z            (input) DOUBLE PRECISION array, dimension (3)
 218:              // *               Each of the elements in z must be positive.
 219:              // *
 220:              // *  FINIT        (input) DOUBLE PRECISION
 221:              // *               The value of f at 0. It is more accurate than the one
 222:              // *               evaluated inside this routine (if someone wants to do
 223:              // *               so).
 224:              // *
 225:              // *  TAU          (output) DOUBLE PRECISION
 226:              // *               The root of the equation f(x).
 227:              // *
 228:              // *  INFO         (output) INTEGER
 229:              // *               = 0: successful exit
 230:              // *               > 0: if INFO = 1, failure to converge
 231:              // *
 232:              // *  Further Details
 233:              // *  ===============
 234:              // *
 235:              // *  30/06/99: Based on contributions by
 236:              // *     Ren-Cang Li, Computer Science Division, University of California
 237:              // *     at Berkeley, USA
 238:              // *
 239:              // *  10/02/03: This version has a few statements commented out for thread
 240:              // *  safety (machine parameters are computed on each entry). SJH.
 241:              // *
 242:              // *  05/10/06: Modified from a new version of Ren-Cang Li, use
 243:              // *     Gragg-Thornton-Warner cubic convergent scheme for better stability.
 244:              // *
 245:              // *  =====================================================================
 246:              // *
 247:              // *     .. Parameters ..
 248:              // *     ..
 249:              // *     .. External Functions ..
 250:              // *     ..
 251:              // *     .. Local Arrays ..
 252:              // *     ..
 253:              // *     .. Local Scalars ..
 254:              // *     ..
 255:              // *     .. Intrinsic Functions ..
 256:              //      INTRINSIC          ABS, INT, LOG, MAX, MIN, SQRT;
 257:              // *     ..
 258:              // *     .. Executable Statements ..
 259:              // *
 260:   
 261:              #endregion
 262:   
 263:   
 264:              #region Body
 265:              
 266:              INFO = 0;
 267:              // *
 268:              if (ORGATI)
 269:              {
 270:                  LBD = D[2 + o_d];
 271:                  UBD = D[3 + o_d];
 272:              }
 273:              else
 274:              {
 275:                  LBD = D[1 + o_d];
 276:                  UBD = D[2 + o_d];
 277:              }
 278:              if (FINIT < ZERO)
 279:              {
 280:                  LBD = ZERO;
 281:              }
 282:              else
 283:              {
 284:                  UBD = ZERO;
 285:              }
 286:              // *
 287:              NITER = 1;
 288:              TAU = ZERO;
 289:              if (KNITER == 2)
 290:              {
 291:                  if (ORGATI)
 292:                  {
 293:                      TEMP = (D[3 + o_d] - D[2 + o_d]) / TWO;
 294:                      C = RHO + Z[1 + o_z] / ((D[1 + o_d] - D[2 + o_d]) - TEMP);
 295:                      A = C * (D[2 + o_d] + D[3 + o_d]) + Z[2 + o_z] + Z[3 + o_z];
 296:                      B = C * D[2 + o_d] * D[3 + o_d] + Z[2 + o_z] * D[3 + o_d] + Z[3 + o_z] * D[2 + o_d];
 297:                  }
 298:                  else
 299:                  {
 300:                      TEMP = (D[1 + o_d] - D[2 + o_d]) / TWO;
 301:                      C = RHO + Z[3 + o_z] / ((D[3 + o_d] - D[2 + o_d]) - TEMP);
 302:                      A = C * (D[1 + o_d] + D[2 + o_d]) + Z[1 + o_z] + Z[2 + o_z];
 303:                      B = C * D[1 + o_d] * D[2 + o_d] + Z[1 + o_z] * D[2 + o_d] + Z[2 + o_z] * D[1 + o_d];
 304:                  }
 305:                  TEMP = Math.Max(Math.Abs(A), Math.Max(Math.Abs(B), Math.Abs(C)));
 306:                  A = A / TEMP;
 307:                  B = B / TEMP;
 308:                  C = C / TEMP;
 309:                  if (C == ZERO)
 310:                  {
 311:                      TAU = B / A;
 312:                  }
 313:                  else
 314:                  {
 315:                      if (A <= ZERO)
 316:                      {
 317:                          TAU = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
 318:                      }
 319:                      else
 320:                      {
 321:                          TAU = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
 322:                      }
 323:                  }
 324:                  if (TAU < LBD || TAU > UBD) TAU = (LBD + UBD) / TWO;
 325:                  if (D[1 + o_d] == TAU || D[2 + o_d] == TAU || D[3 + o_d] == TAU)
 326:                  {
 327:                      TAU = ZERO;
 328:                  }
 329:                  else
 330:                  {
 331:                      TEMP = FINIT + TAU * Z[1 + o_z] / (D[1 + o_d] * (D[1 + o_d] - TAU)) + TAU * Z[2 + o_z] / (D[2 + o_d] * (D[2 + o_d] - TAU)) + TAU * Z[3 + o_z] / (D[3 + o_d] * (D[3 + o_d] - TAU));
 332:                      if (TEMP <= ZERO)
 333:                      {
 334:                          LBD = TAU;
 335:                      }
 336:                      else
 337:                      {
 338:                          UBD = TAU;
 339:                      }
 340:                      if (Math.Abs(FINIT) <= Math.Abs(TEMP)) TAU = ZERO;
 341:                  }
 342:              }
 343:              // *
 344:              // *     get machine parameters for possible scaling to avoid overflow
 345:              // *
 346:              // *     modified by Sven: parameters SMALL1, SMINV1, SMALL2,
 347:              // *     SMINV2, EPS are not SAVEd anymore between one call to the
 348:              // *     others but recomputed at each call
 349:              // *
 350:              EPS = this._dlamch.Run("Epsilon");
 351:              BASE = this._dlamch.Run("Base");
 352:              SMALL1 = Math.Pow(BASE,Convert.ToInt32(Math.Truncate(Math.Log(this._dlamch.Run("SafMin")) / Math.Log(BASE) / THREE)));
 353:              SMINV1 = ONE / SMALL1;
 354:              SMALL2 = SMALL1 * SMALL1;
 355:              SMINV2 = SMINV1 * SMINV1;
 356:              // *
 357:              // *     Determine if scaling of inputs necessary to avoid overflow
 358:              // *     when computing 1/TEMP**3
 359:              // *
 360:              if (ORGATI)
 361:              {
 362:                  TEMP = Math.Min(Math.Abs(D[2 + o_d] - TAU), Math.Abs(D[3 + o_d] - TAU));
 363:              }
 364:              else
 365:              {
 366:                  TEMP = Math.Min(Math.Abs(D[1 + o_d] - TAU), Math.Abs(D[2 + o_d] - TAU));
 367:              }
 368:              SCALE = false;
 369:              if (TEMP <= SMALL1)
 370:              {
 371:                  SCALE = true;
 372:                  if (TEMP <= SMALL2)
 373:                  {
 374:                      // *
 375:                      // *        Scale up by power of radix nearest 1/SAFMIN**(2/3)
 376:                      // *
 377:                      SCLFAC = SMINV2;
 378:                      SCLINV = SMALL2;
 379:                  }
 380:                  else
 381:                  {
 382:                      // *
 383:                      // *        Scale up by power of radix nearest 1/SAFMIN**(1/3)
 384:                      // *
 385:                      SCLFAC = SMINV1;
 386:                      SCLINV = SMALL1;
 387:                  }
 388:                  // *
 389:                  // *        Scaling up safe because D, Z, TAU scaled elsewhere to be O(1)
 390:                  // *
 391:                  for (I = 1; I <= 3; I++)
 392:                  {
 393:                      DSCALE[I + o_dscale] = D[I + o_d] * SCLFAC;
 394:                      ZSCALE[I + o_zscale] = Z[I + o_z] * SCLFAC;
 395:                  }
 396:                  TAU = TAU * SCLFAC;
 397:                  LBD = LBD * SCLFAC;
 398:                  UBD = UBD * SCLFAC;
 399:              }
 400:              else
 401:              {
 402:                  // *
 403:                  // *        Copy D and Z to DSCALE and ZSCALE
 404:                  // *
 405:                  for (I = 1; I <= 3; I++)
 406:                  {
 407:                      DSCALE[I + o_dscale] = D[I + o_d];
 408:                      ZSCALE[I + o_zscale] = Z[I + o_z];
 409:                  }
 410:              }
 411:              // *
 412:              FC = ZERO;
 413:              DF = ZERO;
 414:              DDF = ZERO;
 415:              for (I = 1; I <= 3; I++)
 416:              {
 417:                  TEMP = ONE / (DSCALE[I + o_dscale] - TAU);
 418:                  TEMP1 = ZSCALE[I + o_zscale] * TEMP;
 419:                  TEMP2 = TEMP1 * TEMP;
 420:                  TEMP3 = TEMP2 * TEMP;
 421:                  FC = FC + TEMP1 / DSCALE[I + o_dscale];
 422:                  DF = DF + TEMP2;
 423:                  DDF = DDF + TEMP3;
 424:              }
 425:              F = FINIT + TAU * FC;
 426:              // *
 427:              if (Math.Abs(F) <= ZERO) goto LABEL60;
 428:              if (F <= ZERO)
 429:              {
 430:                  LBD = TAU;
 431:              }
 432:              else
 433:              {
 434:                  UBD = TAU;
 435:              }
 436:              // *
 437:              // *        Iteration begins -- Use Gragg-Thornton-Warner cubic convergent
 438:              // *                            scheme
 439:              // *
 440:              // *     It is not hard to see that
 441:              // *
 442:              // *           1) Iterations will go up monotonically
 443:              // *              if FINIT < 0;
 444:              // *
 445:              // *           2) Iterations will go down monotonically
 446:              // *              if FINIT > 0.
 447:              // *
 448:              ITER = NITER + 1;
 449:              // *
 450:              for (NITER = ITER; NITER <= MAXIT; NITER++)
 451:              {
 452:                  // *
 453:                  if (ORGATI)
 454:                  {
 455:                      TEMP1 = DSCALE[2 + o_dscale] - TAU;
 456:                      TEMP2 = DSCALE[3 + o_dscale] - TAU;
 457:                  }
 458:                  else
 459:                  {
 460:                      TEMP1 = DSCALE[1 + o_dscale] - TAU;
 461:                      TEMP2 = DSCALE[2 + o_dscale] - TAU;
 462:                  }
 463:                  A = (TEMP1 + TEMP2) * F - TEMP1 * TEMP2 * DF;
 464:                  B = TEMP1 * TEMP2 * F;
 465:                  C = F - (TEMP1 + TEMP2) * DF + TEMP1 * TEMP2 * DDF;
 466:                  TEMP = Math.Max(Math.Abs(A), Math.Max(Math.Abs(B), Math.Abs(C)));
 467:                  A = A / TEMP;
 468:                  B = B / TEMP;
 469:                  C = C / TEMP;
 470:                  if (C == ZERO)
 471:                  {
 472:                      ETA = B / A;
 473:                  }
 474:                  else
 475:                  {
 476:                      if (A <= ZERO)
 477:                      {
 478:                          ETA = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);
 479:                      }
 480:                      else
 481:                      {
 482:                          ETA = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));
 483:                      }
 484:                  }
 485:                  if (F * ETA >= ZERO)
 486:                  {
 487:                      ETA =  - F / DF;
 488:                  }
 489:                  // *
 490:                  TAU = TAU + ETA;
 491:                  if (TAU < LBD || TAU > UBD) TAU = (LBD + UBD) / TWO;
 492:                  // *
 493:                  FC = ZERO;
 494:                  ERRETM = ZERO;
 495:                  DF = ZERO;
 496:                  DDF = ZERO;
 497:                  for (I = 1; I <= 3; I++)
 498:                  {
 499:                      TEMP = ONE / (DSCALE[I + o_dscale] - TAU);
 500:                      TEMP1 = ZSCALE[I + o_zscale] * TEMP;
 501:                      TEMP2 = TEMP1 * TEMP;
 502:                      TEMP3 = TEMP2 * TEMP;
 503:                      TEMP4 = TEMP1 / DSCALE[I + o_dscale];
 504:                      FC = FC + TEMP4;
 505:                      ERRETM = ERRETM + Math.Abs(TEMP4);
 506:                      DF = DF + TEMP2;
 507:                      DDF = DDF + TEMP3;
 508:                  }
 509:                  F = FINIT + TAU * FC;
 510:                  ERRETM = EIGHT * (Math.Abs(FINIT) + Math.Abs(TAU) * ERRETM) + Math.Abs(TAU) * DF;
 511:                  if (Math.Abs(F) <= EPS * ERRETM) goto LABEL60;
 512:                  if (F <= ZERO)
 513:                  {
 514:                      LBD = TAU;
 515:                  }
 516:                  else
 517:                  {
 518:                      UBD = TAU;
 519:                  }
 520:              }
 521:              INFO = 1;
 522:          LABEL60:;
 523:              // *
 524:              // *     Undo scaling
 525:              // *
 526:              if (SCALE) TAU = TAU * SCLINV;
 527:              return;
 528:              // *
 529:              // *     End of DLAED6
 530:              // *
 531:   
 532:              #endregion
 533:   
 534:          }
 535:      }
 536:  }