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