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) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DGEHD2 reduces a real general matrix A to upper Hessenberg form H by
  27:      /// an orthogonal similarity transformation:  Q' * A * Q = H .
  28:      /// 
  29:      ///</summary>
  30:      public class DGEHD2
  31:      {
  32:      
  33:   
  34:          #region Dependencies
  35:          
  36:          DLARF _dlarf; DLARFG _dlarfg; XERBLA _xerbla; 
  37:   
  38:          #endregion
  39:   
  40:   
  41:          #region Fields
  42:          
  43:          const double ONE = 1.0E+0; int I = 0; double AII = 0; 
  44:   
  45:          #endregion
  46:   
  47:          public DGEHD2(DLARF dlarf, DLARFG dlarfg, XERBLA xerbla)
  48:          {
  49:      
  50:   
  51:              #region Set Dependencies
  52:              
  53:              this._dlarf = dlarf; this._dlarfg = dlarfg; this._xerbla = xerbla; 
  54:   
  55:              #endregion
  56:   
  57:          }
  58:      
  59:          public DGEHD2()
  60:          {
  61:      
  62:   
  63:              #region Dependencies (Initialization)
  64:              
  65:              LSAME lsame = new LSAME();
  66:              XERBLA xerbla = new XERBLA();
  67:              DLAMC3 dlamc3 = new DLAMC3();
  68:              DLAPY2 dlapy2 = new DLAPY2();
  69:              DNRM2 dnrm2 = new DNRM2();
  70:              DSCAL dscal = new DSCAL();
  71:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  72:              DGER dger = new DGER(xerbla);
  73:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  74:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  75:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  76:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  77:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  78:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  79:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
  80:   
  81:              #endregion
  82:   
  83:   
  84:              #region Set Dependencies
  85:              
  86:              this._dlarf = dlarf; this._dlarfg = dlarfg; this._xerbla = xerbla; 
  87:   
  88:              #endregion
  89:   
  90:          }
  91:          /// <summary>
  92:          /// Purpose
  93:          /// =======
  94:          /// 
  95:          /// DGEHD2 reduces a real general matrix A to upper Hessenberg form H by
  96:          /// an orthogonal similarity transformation:  Q' * A * Q = H .
  97:          /// 
  98:          ///</summary>
  99:          /// <param name="N">
 100:          /// (input) INTEGER
 101:          /// The order of the matrix A.  N .GE. 0.
 102:          ///</param>
 103:          /// <param name="ILO">
 104:          /// (input) INTEGER
 105:          ///</param>
 106:          /// <param name="IHI">
 107:          /// (input) INTEGER
 108:          /// It is assumed that A is already upper triangular in rows
 109:          /// and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
 110:          /// set by a previous call to DGEBAL; otherwise they should be
 111:          /// set to 1 and N respectively. See Further Details.
 112:          /// 1 .LE. ILO .LE. IHI .LE. max(1,N).
 113:          ///</param>
 114:          /// <param name="A">
 115:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 116:          /// On entry, the n by n general matrix to be reduced.
 117:          /// On exit, the upper triangle and the first subdiagonal of A
 118:          /// are overwritten with the upper Hessenberg matrix H, and the
 119:          /// elements below the first subdiagonal, with the array TAU,
 120:          /// represent the orthogonal matrix Q as a product of elementary
 121:          /// reflectors. See Further Details.
 122:          ///</param>
 123:          /// <param name="LDA">
 124:          /// (input) INTEGER
 125:          /// The leading dimension of the array A.  LDA .GE. max(1,N).
 126:          ///</param>
 127:          /// <param name="TAU">
 128:          /// (output) DOUBLE PRECISION array, dimension (N-1)
 129:          /// The scalar factors of the elementary reflectors (see Further
 130:          /// Details).
 131:          ///</param>
 132:          /// <param name="WORK">
 133:          /// (workspace) DOUBLE PRECISION array, dimension (N)
 134:          ///</param>
 135:          /// <param name="INFO">
 136:          /// (output) INTEGER
 137:          /// = 0:  successful exit.
 138:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 139:          ///</param>
 140:          public void Run(int N, int ILO, int IHI, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau
 141:                           , ref double[] WORK, int offset_work, ref int INFO)
 142:          {
 143:   
 144:              #region Array Index Correction
 145:              
 146:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_work = -1 + offset_work; 
 147:   
 148:              #endregion
 149:   
 150:   
 151:              #region Prolog
 152:              
 153:              // *
 154:              // *  -- LAPACK routine (version 3.1) --
 155:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 156:              // *     November 2006
 157:              // *
 158:              // *     .. Scalar Arguments ..
 159:              // *     ..
 160:              // *     .. Array Arguments ..
 161:              // *     ..
 162:              // *
 163:              // *  Purpose
 164:              // *  =======
 165:              // *
 166:              // *  DGEHD2 reduces a real general matrix A to upper Hessenberg form H by
 167:              // *  an orthogonal similarity transformation:  Q' * A * Q = H .
 168:              // *
 169:              // *  Arguments
 170:              // *  =========
 171:              // *
 172:              // *  N       (input) INTEGER
 173:              // *          The order of the matrix A.  N >= 0.
 174:              // *
 175:              // *  ILO     (input) INTEGER
 176:              // *  IHI     (input) INTEGER
 177:              // *          It is assumed that A is already upper triangular in rows
 178:              // *          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
 179:              // *          set by a previous call to DGEBAL; otherwise they should be
 180:              // *          set to 1 and N respectively. See Further Details.
 181:              // *          1 <= ILO <= IHI <= max(1,N).
 182:              // *
 183:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 184:              // *          On entry, the n by n general matrix to be reduced.
 185:              // *          On exit, the upper triangle and the first subdiagonal of A
 186:              // *          are overwritten with the upper Hessenberg matrix H, and the
 187:              // *          elements below the first subdiagonal, with the array TAU,
 188:              // *          represent the orthogonal matrix Q as a product of elementary
 189:              // *          reflectors. See Further Details.
 190:              // *
 191:              // *  LDA     (input) INTEGER
 192:              // *          The leading dimension of the array A.  LDA >= max(1,N).
 193:              // *
 194:              // *  TAU     (output) DOUBLE PRECISION array, dimension (N-1)
 195:              // *          The scalar factors of the elementary reflectors (see Further
 196:              // *          Details).
 197:              // *
 198:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (N)
 199:              // *
 200:              // *  INFO    (output) INTEGER
 201:              // *          = 0:  successful exit.
 202:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 203:              // *
 204:              // *  Further Details
 205:              // *  ===============
 206:              // *
 207:              // *  The matrix Q is represented as a product of (ihi-ilo) elementary
 208:              // *  reflectors
 209:              // *
 210:              // *     Q = H(ilo) H(ilo+1) . . . H(ihi-1).
 211:              // *
 212:              // *  Each H(i) has the form
 213:              // *
 214:              // *     H(i) = I - tau * v * v'
 215:              // *
 216:              // *  where tau is a real scalar, and v is a real vector with
 217:              // *  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
 218:              // *  exit in A(i+2:ihi,i), and tau in TAU(i).
 219:              // *
 220:              // *  The contents of A are illustrated by the following example, with
 221:              // *  n = 7, ilo = 2 and ihi = 6:
 222:              // *
 223:              // *  on entry,                        on exit,
 224:              // *
 225:              // *  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
 226:              // *  (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
 227:              // *  (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
 228:              // *  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
 229:              // *  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
 230:              // *  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
 231:              // *  (                         a )    (                          a )
 232:              // *
 233:              // *  where a denotes an element of the original matrix A, h denotes a
 234:              // *  modified element of the upper Hessenberg matrix H, and vi denotes an
 235:              // *  element of the vector defining H(i).
 236:              // *
 237:              // *  =====================================================================
 238:              // *
 239:              // *     .. Parameters ..
 240:              // *     ..
 241:              // *     .. Local Scalars ..
 242:              // *     ..
 243:              // *     .. External Subroutines ..
 244:              // *     ..
 245:              // *     .. Intrinsic Functions ..
 246:              //      INTRINSIC          MAX, MIN;
 247:              // *     ..
 248:              // *     .. Executable Statements ..
 249:              // *
 250:              // *     Test the input parameters
 251:              // *
 252:   
 253:              #endregion
 254:   
 255:   
 256:              #region Body
 257:              
 258:              INFO = 0;
 259:              if (N < 0)
 260:              {
 261:                  INFO =  - 1;
 262:              }
 263:              else
 264:              {
 265:                  if (ILO < 1 || ILO > Math.Max(1, N))
 266:                  {
 267:                      INFO =  - 2;
 268:                  }
 269:                  else
 270:                  {
 271:                      if (IHI < Math.Min(ILO, N) || IHI > N)
 272:                      {
 273:                          INFO =  - 3;
 274:                      }
 275:                      else
 276:                      {
 277:                          if (LDA < Math.Max(1, N))
 278:                          {
 279:                              INFO =  - 5;
 280:                          }
 281:                      }
 282:                  }
 283:              }
 284:              if (INFO != 0)
 285:              {
 286:                  this._xerbla.Run("DGEHD2",  - INFO);
 287:                  return;
 288:              }
 289:              // *
 290:              for (I = ILO; I <= IHI - 1; I++)
 291:              {
 292:                  // *
 293:                  // *        Compute elementary reflector H(i) to annihilate A(i+2:ihi,i)
 294:                  // *
 295:                  this._dlarfg.Run(IHI - I, ref A[I + 1+I * LDA + o_a], ref A, Math.Min(I + 2, N)+I * LDA + o_a, 1, ref TAU[I + o_tau]);
 296:                  AII = A[I + 1+I * LDA + o_a];
 297:                  A[I + 1+I * LDA + o_a] = ONE;
 298:                  // *
 299:                  // *        Apply H(i) to A(1:ihi,i+1:ihi) from the right
 300:                  // *
 301:                  this._dlarf.Run("Right", IHI, IHI - I, A, I + 1+I * LDA + o_a, 1, TAU[I + o_tau]
 302:                                  , ref A, 1+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
 303:                  // *
 304:                  // *        Apply H(i) to A(i+1:ihi,i+1:n) from the left
 305:                  // *
 306:                  this._dlarf.Run("Left", IHI - I, N - I, A, I + 1+I * LDA + o_a, 1, TAU[I + o_tau]
 307:                                  , ref A, I + 1+(I + 1) * LDA + o_a, LDA, ref WORK, offset_work);
 308:                  // *
 309:                  A[I + 1+I * LDA + o_a] = AII;
 310:              }
 311:              // *
 312:              return;
 313:              // *
 314:              // *     End of DGEHD2
 315:              // *
 316:   
 317:              #endregion
 318:   
 319:          }
 320:      }
 321:  }