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:      /// DLASQ6 computes one dqd (shift equal to zero) transform in
  27:      /// ping-pong form, with protection against underflow and overflow.
  28:      /// 
  29:      ///</summary>
  30:      public class DLASQ6
  31:      {
  32:      
  33:   
  34:          #region Dependencies
  35:          
  36:          DLAMCH _dlamch; 
  37:   
  38:          #endregion
  39:   
  40:   
  41:          #region Fields
  42:          
  43:          const double ZERO = 0.0E0; int J4 = 0; int J4P2 = 0; double D = 0; double EMIN = 0; double SAFMIN = 0; double TEMP = 0; 
  44:   
  45:          #endregion
  46:   
  47:          public DLASQ6(DLAMCH dlamch)
  48:          {
  49:      
  50:   
  51:              #region Set Dependencies
  52:              
  53:              this._dlamch = dlamch; 
  54:   
  55:              #endregion
  56:   
  57:          }
  58:      
  59:          public DLASQ6()
  60:          {
  61:      
  62:   
  63:              #region Dependencies (Initialization)
  64:              
  65:              LSAME lsame = new LSAME();
  66:              DLAMC3 dlamc3 = new DLAMC3();
  67:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  68:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  69:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  70:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  71:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  72:   
  73:              #endregion
  74:   
  75:   
  76:              #region Set Dependencies
  77:              
  78:              this._dlamch = dlamch; 
  79:   
  80:              #endregion
  81:   
  82:          }
  83:          /// <summary>
  84:          /// Purpose
  85:          /// =======
  86:          /// 
  87:          /// DLASQ6 computes one dqd (shift equal to zero) transform in
  88:          /// ping-pong form, with protection against underflow and overflow.
  89:          /// 
  90:          ///</summary>
  91:          /// <param name="I0">
  92:          /// (input) INTEGER
  93:          /// First index.
  94:          ///</param>
  95:          /// <param name="N0">
  96:          /// (input) INTEGER
  97:          /// Last index.
  98:          ///</param>
  99:          /// <param name="Z">
 100:          /// (input) DOUBLE PRECISION array, dimension ( 4*N )
 101:          /// Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
 102:          /// an extra argument.
 103:          ///</param>
 104:          /// <param name="PP">
 105:          /// (input) INTEGER
 106:          /// PP=0 for ping, PP=1 for pong.
 107:          ///</param>
 108:          /// <param name="DMIN">
 109:          /// (output) DOUBLE PRECISION
 110:          /// Minimum value of d.
 111:          ///</param>
 112:          /// <param name="DMIN1">
 113:          /// (output) DOUBLE PRECISION
 114:          /// Minimum value of d, excluding D( N0 ).
 115:          ///</param>
 116:          /// <param name="DMIN2">
 117:          /// (output) DOUBLE PRECISION
 118:          /// Minimum value of d, excluding D( N0 ) and D( N0-1 ).
 119:          ///</param>
 120:          /// <param name="DN">
 121:          /// (output) DOUBLE PRECISION
 122:          /// d(N0), the last value of d.
 123:          ///</param>
 124:          /// <param name="DNM1">
 125:          /// (output) DOUBLE PRECISION
 126:          /// d(N0-1).
 127:          ///</param>
 128:          /// <param name="DNM2">
 129:          /// (output) DOUBLE PRECISION
 130:          /// d(N0-2).
 131:          ///</param>
 132:          public void Run(int I0, int N0, ref double[] Z, int offset_z, int PP, ref double DMIN, ref double DMIN1
 133:                           , ref double DMIN2, ref double DN, ref double DNM1, ref double DNM2)
 134:          {
 135:   
 136:              #region Array Index Correction
 137:              
 138:               int o_z = -1 + offset_z; 
 139:   
 140:              #endregion
 141:   
 142:   
 143:              #region Prolog
 144:              
 145:              // *
 146:              // *  -- LAPACK auxiliary routine (version 3.1) --
 147:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 148:              // *     November 2006
 149:              // *
 150:              // *     .. Scalar Arguments ..
 151:              // *     ..
 152:              // *     .. Array Arguments ..
 153:              // *     ..
 154:              // *
 155:              // *  Purpose
 156:              // *  =======
 157:              // *
 158:              // *  DLASQ6 computes one dqd (shift equal to zero) transform in
 159:              // *  ping-pong form, with protection against underflow and overflow.
 160:              // *
 161:              // *  Arguments
 162:              // *  =========
 163:              // *
 164:              // *  I0    (input) INTEGER
 165:              // *        First index.
 166:              // *
 167:              // *  N0    (input) INTEGER
 168:              // *        Last index.
 169:              // *
 170:              // *  Z     (input) DOUBLE PRECISION array, dimension ( 4*N )
 171:              // *        Z holds the qd array. EMIN is stored in Z(4*N0) to avoid
 172:              // *        an extra argument.
 173:              // *
 174:              // *  PP    (input) INTEGER
 175:              // *        PP=0 for ping, PP=1 for pong.
 176:              // *
 177:              // *  DMIN  (output) DOUBLE PRECISION
 178:              // *        Minimum value of d.
 179:              // *
 180:              // *  DMIN1 (output) DOUBLE PRECISION
 181:              // *        Minimum value of d, excluding D( N0 ).
 182:              // *
 183:              // *  DMIN2 (output) DOUBLE PRECISION
 184:              // *        Minimum value of d, excluding D( N0 ) and D( N0-1 ).
 185:              // *
 186:              // *  DN    (output) DOUBLE PRECISION
 187:              // *        d(N0), the last value of d.
 188:              // *
 189:              // *  DNM1  (output) DOUBLE PRECISION
 190:              // *        d(N0-1).
 191:              // *
 192:              // *  DNM2  (output) DOUBLE PRECISION
 193:              // *        d(N0-2).
 194:              // *
 195:              // *  =====================================================================
 196:              // *
 197:              // *     .. Parameter ..
 198:              // *     ..
 199:              // *     .. Local Scalars ..
 200:              // *     ..
 201:              // *     .. External Function ..
 202:              // *     ..
 203:              // *     .. Intrinsic Functions ..
 204:              //      INTRINSIC          MIN;
 205:              // *     ..
 206:              // *     .. Executable Statements ..
 207:              // *
 208:   
 209:              #endregion
 210:   
 211:   
 212:              #region Body
 213:              
 214:              if ((N0 - I0 - 1) <= 0) return;
 215:              // *
 216:              SAFMIN = this._dlamch.Run("Safe minimum");
 217:              J4 = 4 * I0 + PP - 3;
 218:              EMIN = Z[J4 + 4 + o_z];
 219:              D = Z[J4 + o_z];
 220:              DMIN = D;
 221:              // *
 222:              if (PP == 0)
 223:              {
 224:                  for (J4 = 4 * I0; J4 <= 4 * (N0 - 3); J4 += 4)
 225:                  {
 226:                      Z[J4 - 2 + o_z] = D + Z[J4 - 1 + o_z];
 227:                      if (Z[J4 - 2 + o_z] == ZERO)
 228:                      {
 229:                          Z[J4 + o_z] = ZERO;
 230:                          D = Z[J4 + 1 + o_z];
 231:                          DMIN = D;
 232:                          EMIN = ZERO;
 233:                      }
 234:                      else
 235:                      {
 236:                          if (SAFMIN * Z[J4 + 1 + o_z] < Z[J4 - 2 + o_z] && SAFMIN * Z[J4 - 2 + o_z] < Z[J4 + 1 + o_z])
 237:                          {
 238:                              TEMP = Z[J4 + 1 + o_z] / Z[J4 - 2 + o_z];
 239:                              Z[J4 + o_z] = Z[J4 - 1 + o_z] * TEMP;
 240:                              D = D * TEMP;
 241:                          }
 242:                          else
 243:                          {
 244:                              Z[J4 + o_z] = Z[J4 + 1 + o_z] * (Z[J4 - 1 + o_z] / Z[J4 - 2 + o_z]);
 245:                              D = Z[J4 + 1 + o_z] * (D / Z[J4 - 2 + o_z]);
 246:                          }
 247:                      }
 248:                      DMIN = Math.Min(DMIN, D);
 249:                      EMIN = Math.Min(EMIN, Z[J4 + o_z]);
 250:                  }
 251:              }
 252:              else
 253:              {
 254:                  for (J4 = 4 * I0; J4 <= 4 * (N0 - 3); J4 += 4)
 255:                  {
 256:                      Z[J4 - 3 + o_z] = D + Z[J4 + o_z];
 257:                      if (Z[J4 - 3 + o_z] == ZERO)
 258:                      {
 259:                          Z[J4 - 1 + o_z] = ZERO;
 260:                          D = Z[J4 + 2 + o_z];
 261:                          DMIN = D;
 262:                          EMIN = ZERO;
 263:                      }
 264:                      else
 265:                      {
 266:                          if (SAFMIN * Z[J4 + 2 + o_z] < Z[J4 - 3 + o_z] && SAFMIN * Z[J4 - 3 + o_z] < Z[J4 + 2 + o_z])
 267:                          {
 268:                              TEMP = Z[J4 + 2 + o_z] / Z[J4 - 3 + o_z];
 269:                              Z[J4 - 1 + o_z] = Z[J4 + o_z] * TEMP;
 270:                              D = D * TEMP;
 271:                          }
 272:                          else
 273:                          {
 274:                              Z[J4 - 1 + o_z] = Z[J4 + 2 + o_z] * (Z[J4 + o_z] / Z[J4 - 3 + o_z]);
 275:                              D = Z[J4 + 2 + o_z] * (D / Z[J4 - 3 + o_z]);
 276:                          }
 277:                      }
 278:                      DMIN = Math.Min(DMIN, D);
 279:                      EMIN = Math.Min(EMIN, Z[J4 - 1 + o_z]);
 280:                  }
 281:              }
 282:              // *
 283:              // *     Unroll last two steps. 
 284:              // *
 285:              DNM2 = D;
 286:              DMIN2 = DMIN;
 287:              J4 = 4 * (N0 - 2) - PP;
 288:              J4P2 = J4 + 2 * PP - 1;
 289:              Z[J4 - 2 + o_z] = DNM2 + Z[J4P2 + o_z];
 290:              if (Z[J4 - 2 + o_z] == ZERO)
 291:              {
 292:                  Z[J4 + o_z] = ZERO;
 293:                  DNM1 = Z[J4P2 + 2 + o_z];
 294:                  DMIN = DNM1;
 295:                  EMIN = ZERO;
 296:              }
 297:              else
 298:              {
 299:                  if (SAFMIN * Z[J4P2 + 2 + o_z] < Z[J4 - 2 + o_z] && SAFMIN * Z[J4 - 2 + o_z] < Z[J4P2 + 2 + o_z])
 300:                  {
 301:                      TEMP = Z[J4P2 + 2 + o_z] / Z[J4 - 2 + o_z];
 302:                      Z[J4 + o_z] = Z[J4P2 + o_z] * TEMP;
 303:                      DNM1 = DNM2 * TEMP;
 304:                  }
 305:                  else
 306:                  {
 307:                      Z[J4 + o_z] = Z[J4P2 + 2 + o_z] * (Z[J4P2 + o_z] / Z[J4 - 2 + o_z]);
 308:                      DNM1 = Z[J4P2 + 2 + o_z] * (DNM2 / Z[J4 - 2 + o_z]);
 309:                  }
 310:              }
 311:              DMIN = Math.Min(DMIN, DNM1);
 312:              // *
 313:              DMIN1 = DMIN;
 314:              J4 = J4 + 4;
 315:              J4P2 = J4 + 2 * PP - 1;
 316:              Z[J4 - 2 + o_z] = DNM1 + Z[J4P2 + o_z];
 317:              if (Z[J4 - 2 + o_z] == ZERO)
 318:              {
 319:                  Z[J4 + o_z] = ZERO;
 320:                  DN = Z[J4P2 + 2 + o_z];
 321:                  DMIN = DN;
 322:                  EMIN = ZERO;
 323:              }
 324:              else
 325:              {
 326:                  if (SAFMIN * Z[J4P2 + 2 + o_z] < Z[J4 - 2 + o_z] && SAFMIN * Z[J4 - 2 + o_z] < Z[J4P2 + 2 + o_z])
 327:                  {
 328:                      TEMP = Z[J4P2 + 2 + o_z] / Z[J4 - 2 + o_z];
 329:                      Z[J4 + o_z] = Z[J4P2 + o_z] * TEMP;
 330:                      DN = DNM1 * TEMP;
 331:                  }
 332:                  else
 333:                  {
 334:                      Z[J4 + o_z] = Z[J4P2 + 2 + o_z] * (Z[J4P2 + o_z] / Z[J4 - 2 + o_z]);
 335:                      DN = Z[J4P2 + 2 + o_z] * (DNM1 / Z[J4 - 2 + o_z]);
 336:                  }
 337:              }
 338:              DMIN = Math.Min(DMIN, DN);
 339:              // *
 340:              Z[J4 + 2 + o_z] = DN;
 341:              Z[4 * N0 - PP + o_z] = EMIN;
 342:              return;
 343:              // *
 344:              // *     End of DLASQ6
 345:              // *
 346:   
 347:              #endregion
 348:   
 349:          }
 350:      }
 351:  }