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:      /// DLAZQ3 checks for deflation, computes a shift (TAU) and calls dqds.
  27:      /// In case of failure it changes shifts, and tries again until output
  28:      /// is positive.
  29:      /// 
  30:      ///</summary>
  31:      public class DLAZQ3
  32:      {
  33:      
  34:   
  35:          #region Dependencies
  36:          
  37:          DLASQ5 _dlasq5; DLASQ6 _dlasq6; DLAZQ4 _dlazq4; DLAMCH _dlamch; 
  38:   
  39:          #endregion
  40:   
  41:   
  42:          #region Fields
  43:          
  44:          const double CBIAS = 1.50E0; const double ZERO = 0.0E0; const double QURTR = 0.250E0; const double HALF = 0.5E0; 
  45:          const double ONE = 1.0E0;const double TWO = 2.0E0; const double HUNDRD = 100.0E0; int IPN4 = 0; int J4 = 0; int N0IN = 0; 
  46:          int NN = 0;double EPS = 0; double G = 0; double S = 0; double SAFMIN = 0; double T = 0; double TEMP = 0; double TOL = 0; 
  47:          double TOL2 = 0;
  48:   
  49:          #endregion
  50:   
  51:          public DLAZQ3(DLASQ5 dlasq5, DLASQ6 dlasq6, DLAZQ4 dlazq4, DLAMCH dlamch)
  52:          {
  53:      
  54:   
  55:              #region Set Dependencies
  56:              
  57:              this._dlasq5 = dlasq5; this._dlasq6 = dlasq6; this._dlazq4 = dlazq4; this._dlamch = dlamch; 
  58:   
  59:              #endregion
  60:   
  61:          }
  62:      
  63:          public DLAZQ3()
  64:          {
  65:      
  66:   
  67:              #region Dependencies (Initialization)
  68:              
  69:              DLASQ5 dlasq5 = new DLASQ5();
  70:              LSAME lsame = new LSAME();
  71:              DLAMC3 dlamc3 = new DLAMC3();
  72:              DLAZQ4 dlazq4 = new DLAZQ4();
  73:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  74:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  75:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  76:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  77:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  78:              DLASQ6 dlasq6 = new DLASQ6(dlamch);
  79:   
  80:              #endregion
  81:   
  82:   
  83:              #region Set Dependencies
  84:              
  85:              this._dlasq5 = dlasq5; this._dlasq6 = dlasq6; this._dlazq4 = dlazq4; this._dlamch = dlamch; 
  86:   
  87:              #endregion
  88:   
  89:          }
  90:          /// <summary>
  91:          /// Purpose
  92:          /// =======
  93:          /// 
  94:          /// DLAZQ3 checks for deflation, computes a shift (TAU) and calls dqds.
  95:          /// In case of failure it changes shifts, and tries again until output
  96:          /// is positive.
  97:          /// 
  98:          ///</summary>
  99:          /// <param name="I0">
 100:          /// (input) INTEGER
 101:          /// First index.
 102:          ///</param>
 103:          /// <param name="N0">
 104:          /// (input) INTEGER
 105:          /// Last index.
 106:          ///</param>
 107:          /// <param name="Z">
 108:          /// (input) DOUBLE PRECISION array, dimension ( 4*N )
 109:          /// Z holds the qd array.
 110:          ///</param>
 111:          /// <param name="PP">
 112:          /// (input) INTEGER
 113:          /// PP=0 for ping, PP=1 for pong.
 114:          ///</param>
 115:          /// <param name="DMIN">
 116:          /// (output) DOUBLE PRECISION
 117:          /// Minimum value of d.
 118:          ///</param>
 119:          /// <param name="SIGMA">
 120:          /// (output) DOUBLE PRECISION
 121:          /// Sum of shifts used in current segment.
 122:          ///</param>
 123:          /// <param name="DESIG">
 124:          /// (input/output) DOUBLE PRECISION
 125:          /// Lower order part of SIGMA
 126:          ///</param>
 127:          /// <param name="QMAX">
 128:          /// (input) DOUBLE PRECISION
 129:          /// Maximum value of q.
 130:          ///</param>
 131:          /// <param name="NFAIL">
 132:          /// (output) INTEGER
 133:          /// Number of times shift was too big.
 134:          ///</param>
 135:          /// <param name="ITER">
 136:          /// (output) INTEGER
 137:          /// Number of iterations.
 138:          ///</param>
 139:          /// <param name="NDIV">
 140:          /// (output) INTEGER
 141:          /// Number of divisions.
 142:          ///</param>
 143:          /// <param name="IEEE">
 144:          /// (input) LOGICAL
 145:          /// Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
 146:          ///</param>
 147:          /// <param name="TTYPE">
 148:          /// (input/output) INTEGER
 149:          /// Shift type.  TTYPE is passed as an argument in order to save
 150:          /// its value between calls to DLAZQ3
 151:          ///</param>
 152:          /// <param name="DMIN1">
 153:          /// (input/output) REAL
 154:          ///</param>
 155:          /// <param name="DMIN2">
 156:          /// (input/output) REAL
 157:          ///</param>
 158:          /// <param name="DN">
 159:          /// (input/output) REAL
 160:          ///</param>
 161:          /// <param name="DN1">
 162:          /// (input/output) REAL
 163:          ///</param>
 164:          /// <param name="DN2">
 165:          /// (input/output) REAL
 166:          ///</param>
 167:          /// <param name="TAU">
 168:          /// (input/output) REAL
 169:          /// These are passed as arguments in order to save their values
 170:          /// between calls to DLAZQ3
 171:          ///</param>
 172:          public void Run(int I0, ref int N0, ref double[] Z, int offset_z, int PP, ref double DMIN, ref double SIGMA
 173:                           , ref double DESIG, ref double QMAX, ref int NFAIL, ref int ITER, ref int NDIV, bool IEEE
 174:                           , ref int TTYPE, ref double DMIN1, ref double DMIN2, ref double DN, ref double DN1, ref double DN2
 175:                           , ref double TAU)
 176:          {
 177:   
 178:              #region Array Index Correction
 179:              
 180:               int o_z = -1 + offset_z; 
 181:   
 182:              #endregion
 183:   
 184:   
 185:              #region Prolog
 186:              
 187:              // *
 188:              // *  -- LAPACK auxiliary routine (version 3.1) --
 189:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 190:              // *     November 2006
 191:              // *
 192:              // *     .. Scalar Arguments ..
 193:              // *     ..
 194:              // *     .. Array Arguments ..
 195:              // *     ..
 196:              // *
 197:              // *  Purpose
 198:              // *  =======
 199:              // *
 200:              // *  DLAZQ3 checks for deflation, computes a shift (TAU) and calls dqds.
 201:              // *  In case of failure it changes shifts, and tries again until output
 202:              // *  is positive.
 203:              // *
 204:              // *  Arguments
 205:              // *  =========
 206:              // *
 207:              // *  I0     (input) INTEGER
 208:              // *         First index.
 209:              // *
 210:              // *  N0     (input) INTEGER
 211:              // *         Last index.
 212:              // *
 213:              // *  Z      (input) DOUBLE PRECISION array, dimension ( 4*N )
 214:              // *         Z holds the qd array.
 215:              // *
 216:              // *  PP     (input) INTEGER
 217:              // *         PP=0 for ping, PP=1 for pong.
 218:              // *
 219:              // *  DMIN   (output) DOUBLE PRECISION
 220:              // *         Minimum value of d.
 221:              // *
 222:              // *  SIGMA  (output) DOUBLE PRECISION
 223:              // *         Sum of shifts used in current segment.
 224:              // *
 225:              // *  DESIG  (input/output) DOUBLE PRECISION
 226:              // *         Lower order part of SIGMA
 227:              // *
 228:              // *  QMAX   (input) DOUBLE PRECISION
 229:              // *         Maximum value of q.
 230:              // *
 231:              // *  NFAIL  (output) INTEGER
 232:              // *         Number of times shift was too big.
 233:              // *
 234:              // *  ITER   (output) INTEGER
 235:              // *         Number of iterations.
 236:              // *
 237:              // *  NDIV   (output) INTEGER
 238:              // *         Number of divisions.
 239:              // *
 240:              // *  IEEE   (input) LOGICAL
 241:              // *         Flag for IEEE or non IEEE arithmetic (passed to DLASQ5).
 242:              // *
 243:              // *  TTYPE  (input/output) INTEGER
 244:              // *         Shift type.  TTYPE is passed as an argument in order to save
 245:              // *         its value between calls to DLAZQ3
 246:              // *
 247:              // *  DMIN1  (input/output) REAL
 248:              // *  DMIN2  (input/output) REAL
 249:              // *  DN     (input/output) REAL
 250:              // *  DN1    (input/output) REAL
 251:              // *  DN2    (input/output) REAL
 252:              // *  TAU    (input/output) REAL
 253:              // *         These are passed as arguments in order to save their values
 254:              // *         between calls to DLAZQ3
 255:              // *
 256:              // *  This is a thread safe version of DLASQ3, which passes TTYPE, DMIN1,
 257:              // *  DMIN2, DN, DN1. DN2 and TAU through the argument list in place of
 258:              // *  declaring them in a SAVE statment.
 259:              // *
 260:              // *  =====================================================================
 261:              // *
 262:              // *     .. Parameters ..
 263:              // *     ..
 264:              // *     .. Local Scalars ..
 265:              // *     ..
 266:              // *     .. External Subroutines ..
 267:              // *     ..
 268:              // *     .. External Function ..
 269:              // *     ..
 270:              // *     .. Intrinsic Functions ..
 271:              //      INTRINSIC          ABS, MIN, SQRT;
 272:              // *     ..
 273:              // *     .. Executable Statements ..
 274:              // *
 275:   
 276:              #endregion
 277:   
 278:   
 279:              #region Body
 280:              
 281:              N0IN = N0;
 282:              EPS = this._dlamch.Run("Precision");
 283:              SAFMIN = this._dlamch.Run("Safe minimum");
 284:              TOL = EPS * HUNDRD;
 285:              TOL2 = Math.Pow(TOL,2);
 286:              G = ZERO;
 287:              // *
 288:              // *     Check for deflation.
 289:              // *
 290:          LABEL10:;
 291:              // *
 292:              if (N0 < I0) return;
 293:              if (N0 == I0) goto LABEL20;
 294:              NN = 4 * N0 + PP;
 295:              if (N0 == (I0 + 1)) goto LABEL40;
 296:              // *
 297:              // *     Check whether E(N0-1) is negligible, 1 eigenvalue.
 298:              // *
 299:              if (Z[NN - 5 + o_z] > TOL2 * (SIGMA + Z[NN - 3 + o_z]) && Z[NN - 2 * PP - 4 + o_z] > TOL2 * Z[NN - 7 + o_z]) goto LABEL30;
 300:              // *
 301:          LABEL20:;
 302:              // *
 303:              Z[4 * N0 - 3 + o_z] = Z[4 * N0 + PP - 3 + o_z] + SIGMA;
 304:              N0 = N0 - 1;
 305:              goto LABEL10;
 306:              // *
 307:              // *     Check  whether E(N0-2) is negligible, 2 eigenvalues.
 308:              // *
 309:          LABEL30:;
 310:              // *
 311:              if (Z[NN - 9 + o_z] > TOL2 * SIGMA && Z[NN - 2 * PP - 8 + o_z] > TOL2 * Z[NN - 11 + o_z]) goto LABEL50;
 312:              // *
 313:          LABEL40:;
 314:              // *
 315:              if (Z[NN - 3 + o_z] > Z[NN - 7 + o_z])
 316:              {
 317:                  S = Z[NN - 3 + o_z];
 318:                  Z[NN - 3 + o_z] = Z[NN - 7 + o_z];
 319:                  Z[NN - 7 + o_z] = S;
 320:              }
 321:              if (Z[NN - 5 + o_z] > Z[NN - 3 + o_z] * TOL2)
 322:              {
 323:                  T = HALF * ((Z[NN - 7 + o_z] - Z[NN - 3 + o_z]) + Z[NN - 5 + o_z]);
 324:                  S = Z[NN - 3 + o_z] * (Z[NN - 5 + o_z] / T);
 325:                  if (S <= T)
 326:                  {
 327:                      S = Z[NN - 3 + o_z] * (Z[NN - 5 + o_z] / (T * (ONE + Math.Sqrt(ONE + S / T))));
 328:                  }
 329:                  else
 330:                  {
 331:                      S = Z[NN - 3 + o_z] * (Z[NN - 5 + o_z] / (T + Math.Sqrt(T) * Math.Sqrt(T + S)));
 332:                  }
 333:                  T = Z[NN - 7 + o_z] + (S + Z[NN - 5 + o_z]);
 334:                  Z[NN - 3 + o_z] = Z[NN - 3 + o_z] * (Z[NN - 7 + o_z] / T);
 335:                  Z[NN - 7 + o_z] = T;
 336:              }
 337:              Z[4 * N0 - 7 + o_z] = Z[NN - 7 + o_z] + SIGMA;
 338:              Z[4 * N0 - 3 + o_z] = Z[NN - 3 + o_z] + SIGMA;
 339:              N0 = N0 - 2;
 340:              goto LABEL10;
 341:              // *
 342:          LABEL50:;
 343:              // *
 344:              // *     Reverse the qd-array, if warranted.
 345:              // *
 346:              if (DMIN <= ZERO || N0 < N0IN)
 347:              {
 348:                  if (CBIAS * Z[4 * I0 + PP - 3 + o_z] < Z[4 * N0 + PP - 3 + o_z])
 349:                  {
 350:                      IPN4 = 4 * (I0 + N0);
 351:                      for (J4 = 4 * I0; J4 <= 2 * (I0 + N0 - 1); J4 += 4)
 352:                      {
 353:                          TEMP = Z[J4 - 3 + o_z];
 354:                          Z[J4 - 3 + o_z] = Z[IPN4 - J4 - 3 + o_z];
 355:                          Z[IPN4 - J4 - 3 + o_z] = TEMP;
 356:                          TEMP = Z[J4 - 2 + o_z];
 357:                          Z[J4 - 2 + o_z] = Z[IPN4 - J4 - 2 + o_z];
 358:                          Z[IPN4 - J4 - 2 + o_z] = TEMP;
 359:                          TEMP = Z[J4 - 1 + o_z];
 360:                          Z[J4 - 1 + o_z] = Z[IPN4 - J4 - 5 + o_z];
 361:                          Z[IPN4 - J4 - 5 + o_z] = TEMP;
 362:                          TEMP = Z[J4 + o_z];
 363:                          Z[J4 + o_z] = Z[IPN4 - J4 - 4 + o_z];
 364:                          Z[IPN4 - J4 - 4 + o_z] = TEMP;
 365:                      }
 366:                      if (N0 - I0 <= 4)
 367:                      {
 368:                          Z[4 * N0 + PP - 1 + o_z] = Z[4 * I0 + PP - 1 + o_z];
 369:                          Z[4 * N0 - PP + o_z] = Z[4 * I0 - PP + o_z];
 370:                      }
 371:                      DMIN2 = Math.Min(DMIN2, Z[4 * N0 + PP - 1 + o_z]);
 372:                      Z[4 * N0 + PP - 1 + o_z] = Math.Min(Z[4 * N0 + PP - 1 + o_z], Math.Min(Z[4 * I0 + PP - 1 + o_z], Z[4 * I0 + PP + 3 + o_z]));
 373:                      Z[4 * N0 - PP + o_z] = Math.Min(Z[4 * N0 - PP + o_z], Math.Min(Z[4 * I0 - PP + o_z], Z[4 * I0 - PP + 4 + o_z]));
 374:                      QMAX = Math.Max(QMAX, Math.Max(Z[4 * I0 + PP - 3 + o_z], Z[4 * I0 + PP + 1 + o_z]));
 375:                      DMIN =  - ZERO;
 376:                  }
 377:              }
 378:              // *
 379:              if (DMIN < ZERO || SAFMIN * QMAX < Math.Min(Z[4 * N0 + PP - 1 + o_z], Math.Min(Z[4 * N0 + PP - 9 + o_z], DMIN2 + Z[4 * N0 - PP + o_z])))
 380:              {
 381:                  // *
 382:                  // *        Choose a shift.
 383:                  // *
 384:                  this._dlazq4.Run(I0, N0, Z, offset_z, PP, N0IN, DMIN
 385:                                   , DMIN1, DMIN2, DN, DN1, DN2, ref TAU
 386:                                   , ref TTYPE, ref G);
 387:                  // *
 388:                  // *        Call dqds until DMIN > 0.
 389:                  // *
 390:              LABEL80:;
 391:                  // *
 392:                  this._dlasq5.Run(I0, N0, ref Z, offset_z, PP, TAU, ref DMIN
 393:                                   , ref DMIN1, ref DMIN2, ref DN, ref DN1, ref DN2, IEEE);
 394:                  // *
 395:                  NDIV = NDIV + (N0 - I0 + 2);
 396:                  ITER = ITER + 1;
 397:                  // *
 398:                  // *        Check status.
 399:                  // *
 400:                  if (DMIN >= ZERO && DMIN1 > ZERO)
 401:                  {
 402:                      // *
 403:                      // *           Success.
 404:                      // *
 405:                      goto LABEL100;
 406:                      // *
 407:                  }
 408:                  else
 409:                  {
 410:                      if (DMIN < ZERO && DMIN1 > ZERO && Z[4 * (N0 - 1) - PP + o_z] < TOL * (SIGMA + DN1) && Math.Abs(DN) < TOL * SIGMA)
 411:                      {
 412:                          // *
 413:                          // *           Convergence hidden by negative DN.
 414:                          // *
 415:                          Z[4 * (N0 - 1) - PP + 2 + o_z] = ZERO;
 416:                          DMIN = ZERO;
 417:                          goto LABEL100;
 418:                      }
 419:                      else
 420:                      {
 421:                          if (DMIN < ZERO)
 422:                          {
 423:                              // *
 424:                              // *           TAU too big. Select new TAU and try again.
 425:                              // *
 426:                              NFAIL = NFAIL + 1;
 427:                              if (TTYPE <  - 22)
 428:                              {
 429:                                  // *
 430:                                  // *              Failed twice. Play it safe.
 431:                                  // *
 432:                                  TAU = ZERO;
 433:                              }
 434:                              else
 435:                              {
 436:                                  if (DMIN1 > ZERO)
 437:                                  {
 438:                                      // *
 439:                                      // *              Late failure. Gives excellent shift.
 440:                                      // *
 441:                                      TAU = (TAU + DMIN) * (ONE - TWO * EPS);
 442:                                      TTYPE = TTYPE - 11;
 443:                                  }
 444:                                  else
 445:                                  {
 446:                                      // *
 447:                                      // *              Early failure. Divide by 4.
 448:                                      // *
 449:                                      TAU = QURTR * TAU;
 450:                                      TTYPE = TTYPE - 12;
 451:                                  }
 452:                              }
 453:                              goto LABEL80;
 454:                          }
 455:                          else
 456:                          {
 457:                              if (DMIN != DMIN)
 458:                              {
 459:                                  // *
 460:                                  // *           NaN.
 461:                                  // *
 462:                                  TAU = ZERO;
 463:                                  goto LABEL80;
 464:                              }
 465:                              else
 466:                              {
 467:                                  // *
 468:                                  // *           Possible underflow. Play it safe.
 469:                                  // *
 470:                                  goto LABEL90;
 471:                              }
 472:                          }
 473:                      }
 474:                  }
 475:              }
 476:              // *
 477:              // *     Risk of underflow.
 478:              // *
 479:          LABEL90:;
 480:              this._dlasq6.Run(I0, N0, ref Z, offset_z, PP, ref DMIN, ref DMIN1
 481:                               , ref DMIN2, ref DN, ref DN1, ref DN2);
 482:              NDIV = NDIV + (N0 - I0 + 2);
 483:              ITER = ITER + 1;
 484:              TAU = ZERO;
 485:              // *
 486:          LABEL100:;
 487:              if (TAU < SIGMA)
 488:              {
 489:                  DESIG = DESIG + TAU;
 490:                  T = SIGMA + DESIG;
 491:                  DESIG = DESIG - (T - SIGMA);
 492:              }
 493:              else
 494:              {
 495:                  T = SIGMA + TAU;
 496:                  DESIG = SIGMA - (T - TAU) + DESIG;
 497:              }
 498:              SIGMA = T;
 499:              // *
 500:              return;
 501:              // *
 502:              // *     End of DLAZQ3
 503:              // *
 504:   
 505:              #endregion
 506:   
 507:          }
 508:      }
 509:  }