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:      ///</summary>
  24:      public class DLAQR1
  25:      {
  26:      
  27:   
  28:          #region Fields
  29:          
  30:          const double ZERO = 0.0E0; double H21S = 0; double H31S = 0; double S = 0; 
  31:   
  32:          #endregion
  33:   
  34:          public DLAQR1()
  35:          {
  36:      
  37:          }
  38:      
  39:          /// <param name="N">
  40:          /// (input) integer
  41:          /// Order of the matrix H. N must be either 2 or 3.
  42:          ///</param>
  43:          /// <param name="H">
  44:          /// (input) DOUBLE PRECISION array of dimension (LDH,N)
  45:          /// The 2-by-2 or 3-by-3 matrix H in (*).
  46:          ///</param>
  47:          /// <param name="LDH">
  48:          /// (input) integer
  49:          /// The leading dimension of H as declared in
  50:          /// the calling procedure.  LDH.GE.N
  51:          ///</param>
  52:          /// <param name="SR1">
  53:          /// (input) DOUBLE PRECISION
  54:          ///</param>
  55:          /// <param name="SI1">
  56:          /// The shifts in (*).
  57:          ///</param>
  58:          /// <param name="V">
  59:          /// (output) DOUBLE PRECISION array of dimension N
  60:          /// A scalar multiple of the first column of the
  61:          /// matrix K in (*).
  62:          ///</param>
  63:          public void Run(int N, double[] H, int offset_h, int LDH, double SR1, double SI1, double SR2
  64:                           , double SI2, ref double[] V, int offset_v)
  65:          {
  66:   
  67:              #region Array Index Correction
  68:              
  69:               int o_h = -1 - LDH + offset_h;  int o_v = -1 + offset_v; 
  70:   
  71:              #endregion
  72:   
  73:   
  74:              #region Prolog
  75:              
  76:              // *
  77:              // *  -- LAPACK auxiliary routine (version 3.1) --
  78:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  79:              // *     November 2006
  80:              // *
  81:              // *     .. Scalar Arguments ..
  82:              // *     ..
  83:              // *     .. Array Arguments ..
  84:              // *     ..
  85:              // *
  86:              // *       Given a 2-by-2 or 3-by-3 matrix H, DLAQR1 sets v to a
  87:              // *       scalar multiple of the first column of the product
  88:              // *
  89:              // *       (*)  K = (H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)
  90:              // *
  91:              // *       scaling to avoid overflows and most underflows. It
  92:              // *       is assumed that either
  93:              // *
  94:              // *               1) sr1 = sr2 and si1 = -si2
  95:              // *           or
  96:              // *               2) si1 = si2 = 0.
  97:              // *
  98:              // *       This is useful for starting double implicit shift bulges
  99:              // *       in the QR algorithm.
 100:              // *
 101:              // *
 102:              // *       N      (input) integer
 103:              // *              Order of the matrix H. N must be either 2 or 3.
 104:              // *
 105:              // *       H      (input) DOUBLE PRECISION array of dimension (LDH,N)
 106:              // *              The 2-by-2 or 3-by-3 matrix H in (*).
 107:              // *
 108:              // *       LDH    (input) integer
 109:              // *              The leading dimension of H as declared in
 110:              // *              the calling procedure.  LDH.GE.N
 111:              // *
 112:              // *       SR1    (input) DOUBLE PRECISION
 113:              // *       SI1    The shifts in (*).
 114:              // *       SR2
 115:              // *       SI2
 116:              // *
 117:              // *       V      (output) DOUBLE PRECISION array of dimension N
 118:              // *              A scalar multiple of the first column of the
 119:              // *              matrix K in (*).
 120:              // *
 121:              // *     ================================================================
 122:              // *     Based on contributions by
 123:              // *        Karen Braman and Ralph Byers, Department of Mathematics,
 124:              // *        University of Kansas, USA
 125:              // *
 126:              // *     ================================================================
 127:              // *
 128:              // *     .. Parameters ..
 129:              // *     ..
 130:              // *     .. Local Scalars ..
 131:              // *     ..
 132:              // *     .. Intrinsic Functions ..
 133:              //      INTRINSIC          ABS;
 134:              // *     ..
 135:              // *     .. Executable Statements ..
 136:   
 137:              #endregion
 138:   
 139:   
 140:              #region Body
 141:              
 142:              if (N == 2)
 143:              {
 144:                  S = Math.Abs(H[1+1 * LDH + o_h] - SR2) + Math.Abs(SI2) + Math.Abs(H[2+1 * LDH + o_h]);
 145:                  if (S == ZERO)
 146:                  {
 147:                      V[1 + o_v] = ZERO;
 148:                      V[2 + o_v] = ZERO;
 149:                  }
 150:                  else
 151:                  {
 152:                      H21S = H[2+1 * LDH + o_h] / S;
 153:                      V[1 + o_v] = H21S * H[1+2 * LDH + o_h] + (H[1+1 * LDH + o_h] - SR1) * ((H[1+1 * LDH + o_h] - SR2) / S) - SI1 * (SI2 / S);
 154:                      V[2 + o_v] = H21S * (H[1+1 * LDH + o_h] + H[2+2 * LDH + o_h] - SR1 - SR2);
 155:                  }
 156:              }
 157:              else
 158:              {
 159:                  S = Math.Abs(H[1+1 * LDH + o_h] - SR2) + Math.Abs(SI2) + Math.Abs(H[2+1 * LDH + o_h]) + Math.Abs(H[3+1 * LDH + o_h]);
 160:                  if (S == ZERO)
 161:                  {
 162:                      V[1 + o_v] = ZERO;
 163:                      V[2 + o_v] = ZERO;
 164:                      V[3 + o_v] = ZERO;
 165:                  }
 166:                  else
 167:                  {
 168:                      H21S = H[2+1 * LDH + o_h] / S;
 169:                      H31S = H[3+1 * LDH + o_h] / S;
 170:                      V[1 + o_v] = (H[1+1 * LDH + o_h] - SR1) * ((H[1+1 * LDH + o_h] - SR2) / S) - SI1 * (SI2 / S) + H[1+2 * LDH + o_h] * H21S + H[1+3 * LDH + o_h] * H31S;
 171:                      V[2 + o_v] = H21S * (H[1+1 * LDH + o_h] + H[2+2 * LDH + o_h] - SR1 - SR2) + H[2+3 * LDH + o_h] * H31S;
 172:                      V[3 + o_v] = H31S * (H[1+1 * LDH + o_h] + H[3+3 * LDH + o_h] - SR1 - SR2) + H21S * H[3+2 * LDH + o_h];
 173:                  }
 174:              }
 175:   
 176:              #endregion
 177:   
 178:          }
 179:      }
 180:  }