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:      /// DGEHRD 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 DGEHRD
  31:      {
  32:      
  33:   
  34:          #region Dependencies
  35:          
  36:          DAXPY _daxpy; DGEHD2 _dgehd2; DGEMM _dgemm; DLAHR2 _dlahr2; DLARFB _dlarfb; DTRMM _dtrmm; XERBLA _xerbla; ILAENV _ilaenv; 
  37:   
  38:          #endregion
  39:   
  40:   
  41:          #region Fields
  42:          
  43:          const int NBMAX = 64; const int LDT = NBMAX + 1; const double ZERO = 0.0E+0; const double ONE = 1.0E+0; 
  44:          bool LQUERY = false;int I = 0; int IB = 0; int IINFO = 0; int IWS = 0; int J = 0; int LDWORK = 0; int LWKOPT = 0; 
  45:          int NB = 0;int NBMIN = 0; int NH = 0; int NX = 0; double EI = 0; double[] T = new double[LDT * NBMAX]; int offset_t = 0;
  46:   
  47:          #endregion
  48:   
  49:          public DGEHRD(DAXPY daxpy, DGEHD2 dgehd2, DGEMM dgemm, DLAHR2 dlahr2, DLARFB dlarfb, DTRMM dtrmm, XERBLA xerbla, ILAENV ilaenv)
  50:          {
  51:      
  52:   
  53:              #region Set Dependencies
  54:              
  55:              this._daxpy = daxpy; this._dgehd2 = dgehd2; this._dgemm = dgemm; this._dlahr2 = dlahr2; this._dlarfb = dlarfb; 
  56:              this._dtrmm = dtrmm;this._xerbla = xerbla; this._ilaenv = ilaenv; 
  57:   
  58:              #endregion
  59:   
  60:          }
  61:      
  62:          public DGEHRD()
  63:          {
  64:      
  65:   
  66:              #region Dependencies (Initialization)
  67:              
  68:              DAXPY daxpy = new DAXPY();
  69:              LSAME lsame = new LSAME();
  70:              XERBLA xerbla = new XERBLA();
  71:              DLAMC3 dlamc3 = new DLAMC3();
  72:              DLAPY2 dlapy2 = new DLAPY2();
  73:              DNRM2 dnrm2 = new DNRM2();
  74:              DSCAL dscal = new DSCAL();
  75:              DCOPY dcopy = new DCOPY();
  76:              IEEECK ieeeck = new IEEECK();
  77:              IPARMQ iparmq = new IPARMQ();
  78:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  79:              DGER dger = new DGER(xerbla);
  80:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  81:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  82:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  83:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  84:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  85:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  86:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
  87:              DGEHD2 dgehd2 = new DGEHD2(dlarf, dlarfg, xerbla);
  88:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  89:              DLACPY dlacpy = new DLACPY(lsame);
  90:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  91:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
  92:              DLAHR2 dlahr2 = new DLAHR2(daxpy, dcopy, dgemm, dgemv, dlacpy, dlarfg, dscal, dtrmm, dtrmv);
  93:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
  94:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  95:   
  96:              #endregion
  97:   
  98:   
  99:              #region Set Dependencies
 100:              
 101:              this._daxpy = daxpy; this._dgehd2 = dgehd2; this._dgemm = dgemm; this._dlahr2 = dlahr2; this._dlarfb = dlarfb; 
 102:              this._dtrmm = dtrmm;this._xerbla = xerbla; this._ilaenv = ilaenv; 
 103:   
 104:              #endregion
 105:   
 106:          }
 107:          /// <summary>
 108:          /// Purpose
 109:          /// =======
 110:          /// 
 111:          /// DGEHRD reduces a real general matrix A to upper Hessenberg form H by
 112:          /// an orthogonal similarity transformation:  Q' * A * Q = H .
 113:          /// 
 114:          ///</summary>
 115:          /// <param name="N">
 116:          /// (input) INTEGER
 117:          /// The order of the matrix A.  N .GE. 0.
 118:          ///</param>
 119:          /// <param name="ILO">
 120:          /// (input) INTEGER
 121:          ///</param>
 122:          /// <param name="IHI">
 123:          /// (input) INTEGER
 124:          /// It is assumed that A is already upper triangular in rows
 125:          /// and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
 126:          /// set by a previous call to DGEBAL; otherwise they should be
 127:          /// set to 1 and N respectively. See Further Details.
 128:          /// 1 .LE. ILO .LE. IHI .LE. N, if N .GT. 0; ILO=1 and IHI=0, if N=0.
 129:          ///</param>
 130:          /// <param name="A">
 131:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 132:          /// On entry, the N-by-N general matrix to be reduced.
 133:          /// On exit, the upper triangle and the first subdiagonal of A
 134:          /// are overwritten with the upper Hessenberg matrix H, and the
 135:          /// elements below the first subdiagonal, with the array TAU,
 136:          /// represent the orthogonal matrix Q as a product of elementary
 137:          /// reflectors. See Further Details.
 138:          ///</param>
 139:          /// <param name="LDA">
 140:          /// (input) INTEGER
 141:          /// The leading dimension of the array A.  LDA .GE. max(1,N).
 142:          ///</param>
 143:          /// <param name="TAU">
 144:          /// (output) DOUBLE PRECISION array, dimension (N-1)
 145:          /// The scalar factors of the elementary reflectors (see Further
 146:          /// Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
 147:          /// zero.
 148:          ///</param>
 149:          /// <param name="WORK">
 150:          /// (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
 151:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 152:          ///</param>
 153:          /// <param name="LWORK">
 154:          /// (input) INTEGER
 155:          /// The length of the array WORK.  LWORK .GE. max(1,N).
 156:          /// For optimum performance LWORK .GE. N*NB, where NB is the
 157:          /// optimal blocksize.
 158:          /// 
 159:          /// If LWORK = -1, then a workspace query is assumed; the routine
 160:          /// only calculates the optimal size of the WORK array, returns
 161:          /// this value as the first entry of the WORK array, and no error
 162:          /// message related to LWORK is issued by XERBLA.
 163:          ///</param>
 164:          /// <param name="INFO">
 165:          /// (output) INTEGER
 166:          /// = 0:  successful exit
 167:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 168:          ///</param>
 169:          public void Run(int N, int ILO, int IHI, ref double[] A, int offset_a, int LDA, ref double[] TAU, int offset_tau
 170:                           , ref double[] WORK, int offset_work, int LWORK, ref int INFO)
 171:          {
 172:   
 173:              #region Array Index Correction
 174:              
 175:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_work = -1 + offset_work; 
 176:   
 177:              #endregion
 178:   
 179:   
 180:              #region Prolog
 181:              
 182:              // *
 183:              // *  -- LAPACK routine (version 3.1) --
 184:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 185:              // *     November 2006
 186:              // *
 187:              // *     .. Scalar Arguments ..
 188:              // *     ..
 189:              // *     .. Array Arguments ..
 190:              // *     ..
 191:              // *
 192:              // *  Purpose
 193:              // *  =======
 194:              // *
 195:              // *  DGEHRD reduces a real general matrix A to upper Hessenberg form H by
 196:              // *  an orthogonal similarity transformation:  Q' * A * Q = H .
 197:              // *
 198:              // *  Arguments
 199:              // *  =========
 200:              // *
 201:              // *  N       (input) INTEGER
 202:              // *          The order of the matrix A.  N >= 0.
 203:              // *
 204:              // *  ILO     (input) INTEGER
 205:              // *  IHI     (input) INTEGER
 206:              // *          It is assumed that A is already upper triangular in rows
 207:              // *          and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally
 208:              // *          set by a previous call to DGEBAL; otherwise they should be
 209:              // *          set to 1 and N respectively. See Further Details.
 210:              // *          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
 211:              // *
 212:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 213:              // *          On entry, the N-by-N general matrix to be reduced.
 214:              // *          On exit, the upper triangle and the first subdiagonal of A
 215:              // *          are overwritten with the upper Hessenberg matrix H, and the
 216:              // *          elements below the first subdiagonal, with the array TAU,
 217:              // *          represent the orthogonal matrix Q as a product of elementary
 218:              // *          reflectors. See Further Details.
 219:              // *
 220:              // *  LDA     (input) INTEGER
 221:              // *          The leading dimension of the array A.  LDA >= max(1,N).
 222:              // *
 223:              // *  TAU     (output) DOUBLE PRECISION array, dimension (N-1)
 224:              // *          The scalar factors of the elementary reflectors (see Further
 225:              // *          Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to
 226:              // *          zero.
 227:              // *
 228:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
 229:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 230:              // *
 231:              // *  LWORK   (input) INTEGER
 232:              // *          The length of the array WORK.  LWORK >= max(1,N).
 233:              // *          For optimum performance LWORK >= N*NB, where NB is the
 234:              // *          optimal blocksize.
 235:              // *
 236:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 237:              // *          only calculates the optimal size of the WORK array, returns
 238:              // *          this value as the first entry of the WORK array, and no error
 239:              // *          message related to LWORK is issued by XERBLA.
 240:              // *
 241:              // *  INFO    (output) INTEGER
 242:              // *          = 0:  successful exit
 243:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 244:              // *
 245:              // *  Further Details
 246:              // *  ===============
 247:              // *
 248:              // *  The matrix Q is represented as a product of (ihi-ilo) elementary
 249:              // *  reflectors
 250:              // *
 251:              // *     Q = H(ilo) H(ilo+1) . . . H(ihi-1).
 252:              // *
 253:              // *  Each H(i) has the form
 254:              // *
 255:              // *     H(i) = I - tau * v * v'
 256:              // *
 257:              // *  where tau is a real scalar, and v is a real vector with
 258:              // *  v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on
 259:              // *  exit in A(i+2:ihi,i), and tau in TAU(i).
 260:              // *
 261:              // *  The contents of A are illustrated by the following example, with
 262:              // *  n = 7, ilo = 2 and ihi = 6:
 263:              // *
 264:              // *  on entry,                        on exit,
 265:              // *
 266:              // *  ( a   a   a   a   a   a   a )    (  a   a   h   h   h   h   a )
 267:              // *  (     a   a   a   a   a   a )    (      a   h   h   h   h   a )
 268:              // *  (     a   a   a   a   a   a )    (      h   h   h   h   h   h )
 269:              // *  (     a   a   a   a   a   a )    (      v2  h   h   h   h   h )
 270:              // *  (     a   a   a   a   a   a )    (      v2  v3  h   h   h   h )
 271:              // *  (     a   a   a   a   a   a )    (      v2  v3  v4  h   h   h )
 272:              // *  (                         a )    (                          a )
 273:              // *
 274:              // *  where a denotes an element of the original matrix A, h denotes a
 275:              // *  modified element of the upper Hessenberg matrix H, and vi denotes an
 276:              // *  element of the vector defining H(i).
 277:              // *
 278:              // *  This file is a slight modification of LAPACK-3.0's DGEHRD
 279:              // *  subroutine incorporating improvements proposed by Quintana-Orti and
 280:              // *  Van de Geijn (2005). 
 281:              // *
 282:              // *  =====================================================================
 283:              // *
 284:              // *     .. Parameters ..
 285:              // *     ..
 286:              // *     .. Local Scalars ..
 287:              // *     ..
 288:              // *     .. Local Arrays ..
 289:              // *     ..
 290:              // *     .. External Subroutines ..
 291:              // *     ..
 292:              // *     .. Intrinsic Functions ..
 293:              //      INTRINSIC          MAX, MIN;
 294:              // *     ..
 295:              // *     .. External Functions ..
 296:              // *     ..
 297:              // *     .. Executable Statements ..
 298:              // *
 299:              // *     Test the input parameters
 300:              // *
 301:   
 302:              #endregion
 303:   
 304:   
 305:              #region Body
 306:              
 307:              INFO = 0;
 308:              NB = Math.Min(NBMAX, this._ilaenv.Run(1, "DGEHRD", " ", N, ILO, IHI,  - 1));
 309:              LWKOPT = N * NB;
 310:              WORK[1 + o_work] = LWKOPT;
 311:              LQUERY = (LWORK ==  - 1);
 312:              if (N < 0)
 313:              {
 314:                  INFO =  - 1;
 315:              }
 316:              else
 317:              {
 318:                  if (ILO < 1 || ILO > Math.Max(1, N))
 319:                  {
 320:                      INFO =  - 2;
 321:                  }
 322:                  else
 323:                  {
 324:                      if (IHI < Math.Min(ILO, N) || IHI > N)
 325:                      {
 326:                          INFO =  - 3;
 327:                      }
 328:                      else
 329:                      {
 330:                          if (LDA < Math.Max(1, N))
 331:                          {
 332:                              INFO =  - 5;
 333:                          }
 334:                          else
 335:                          {
 336:                              if (LWORK < Math.Max(1, N) && !LQUERY)
 337:                              {
 338:                                  INFO =  - 8;
 339:                              }
 340:                          }
 341:                      }
 342:                  }
 343:              }
 344:              if (INFO != 0)
 345:              {
 346:                  this._xerbla.Run("DGEHRD",  - INFO);
 347:                  return;
 348:              }
 349:              else
 350:              {
 351:                  if (LQUERY)
 352:                  {
 353:                      return;
 354:                  }
 355:              }
 356:              // *
 357:              // *     Set elements 1:ILO-1 and IHI:N-1 of TAU to zero
 358:              // *
 359:              for (I = 1; I <= ILO - 1; I++)
 360:              {
 361:                  TAU[I + o_tau] = ZERO;
 362:              }
 363:              for (I = Math.Max(1, IHI); I <= N - 1; I++)
 364:              {
 365:                  TAU[I + o_tau] = ZERO;
 366:              }
 367:              // *
 368:              // *     Quick return if possible
 369:              // *
 370:              NH = IHI - ILO + 1;
 371:              if (NH <= 1)
 372:              {
 373:                  WORK[1 + o_work] = 1;
 374:                  return;
 375:              }
 376:              // *
 377:              // *     Determine the block size
 378:              // *
 379:              NB = Math.Min(NBMAX, this._ilaenv.Run(1, "DGEHRD", " ", N, ILO, IHI,  - 1));
 380:              NBMIN = 2;
 381:              IWS = 1;
 382:              if (NB > 1 && NB < NH)
 383:              {
 384:                  // *
 385:                  // *        Determine when to cross over from blocked to unblocked code
 386:                  // *        (last block is always handled by unblocked code)
 387:                  // *
 388:                  NX = Math.Max(NB, this._ilaenv.Run(3, "DGEHRD", " ", N, ILO, IHI,  - 1));
 389:                  if (NX < NH)
 390:                  {
 391:                      // *
 392:                      // *           Determine if workspace is large enough for blocked code
 393:                      // *
 394:                      IWS = N * NB;
 395:                      if (LWORK < IWS)
 396:                      {
 397:                          // *
 398:                          // *              Not enough workspace to use optimal NB:  determine the
 399:                          // *              minimum value of NB, and reduce NB or force use of
 400:                          // *              unblocked code
 401:                          // *
 402:                          NBMIN = Math.Max(2, this._ilaenv.Run(2, "DGEHRD", " ", N, ILO, IHI,  - 1));
 403:                          if (LWORK >= N * NBMIN)
 404:                          {
 405:                              NB = LWORK / N;
 406:                          }
 407:                          else
 408:                          {
 409:                              NB = 1;
 410:                          }
 411:                      }
 412:                  }
 413:              }
 414:              LDWORK = N;
 415:              // *
 416:              if (NB < NBMIN || NB >= NH)
 417:              {
 418:                  // *
 419:                  // *        Use unblocked code below
 420:                  // *
 421:                  I = ILO;
 422:                  // *
 423:              }
 424:              else
 425:              {
 426:                  // *
 427:                  // *        Use blocked code
 428:                  // *
 429:                  for (I = ILO; (NB >= 0) ? (I <= IHI - 1 - NX) : (I >= IHI - 1 - NX); I += NB)
 430:                  {
 431:                      IB = Math.Min(NB, IHI - I);
 432:                      // *
 433:                      // *           Reduce columns i:i+ib-1 to Hessenberg form, returning the
 434:                      // *           matrices V and T of the block reflector H = I - V*T*V'
 435:                      // *           which performs the reduction, and also the matrix Y = A*V*T
 436:                      // *
 437:                      this._dlahr2.Run(IHI, I, IB, ref A, 1+I * LDA + o_a, LDA, ref TAU, I + o_tau
 438:                                       , ref T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 439:                      // *
 440:                      // *           Apply the block reflector H to A(1:ihi,i+ib:ihi) from the
 441:                      // *           right, computing  A := A - Y * V'. V(i+ib,ib-1) must be set
 442:                      // *           to 1
 443:                      // *
 444:                      EI = A[I + IB+(I + IB - 1) * LDA + o_a];
 445:                      A[I + IB+(I + IB - 1) * LDA + o_a] = ONE;
 446:                      this._dgemm.Run("No transpose", "Transpose", IHI, IHI - I - IB + 1, IB,  - ONE
 447:                                      , WORK, offset_work, LDWORK, A, I + IB+I * LDA + o_a, LDA, ONE, ref A, 1+(I + IB) * LDA + o_a
 448:                                      , LDA);
 449:                      A[I + IB+(I + IB - 1) * LDA + o_a] = EI;
 450:                      // *
 451:                      // *           Apply the block reflector H to A(1:i,i+1:i+ib-1) from the
 452:                      // *           right
 453:                      // *
 454:                      this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", I, IB - 1
 455:                                      , ONE, A, I + 1+I * LDA + o_a, LDA, ref WORK, offset_work, LDWORK);
 456:                      for (J = 0; J <= IB - 2; J++)
 457:                      {
 458:                          this._daxpy.Run(I,  - ONE, WORK, LDWORK * J + 1 + o_work, 1, ref A, 1+(I + J + 1) * LDA + o_a, 1);
 459:                      }
 460:                      // *
 461:                      // *           Apply the block reflector H to A(i+1:ihi,i+ib:n) from the
 462:                      // *           left
 463:                      // *
 464:                      this._dlarfb.Run("Left", "Transpose", "Forward", "Columnwise", IHI - I, N - I - IB + 1
 465:                                       , IB, A, I + 1+I * LDA + o_a, LDA, T, offset_t, LDT, ref A, I + 1+(I + IB) * LDA + o_a
 466:                                       , LDA, ref WORK, offset_work, LDWORK);
 467:                  }
 468:              }
 469:              // *
 470:              // *     Use unblocked code to reduce the rest of the matrix
 471:              // *
 472:              this._dgehd2.Run(N, I, IHI, ref A, offset_a, LDA, ref TAU, offset_tau
 473:                               , ref WORK, offset_work, ref IINFO);
 474:              WORK[1 + o_work] = IWS;
 475:              // *
 476:              return;
 477:              // *
 478:              // *     End of DGEHRD
 479:              // *
 480:   
 481:              #endregion
 482:   
 483:          }
 484:      }
 485:  }