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:      /// DORGHR generates a real orthogonal matrix Q which is defined as the
  27:      /// product of IHI-ILO elementary reflectors of order N, as returned by
  28:      /// DGEHRD:
  29:      /// 
  30:      /// Q = H(ilo) H(ilo+1) . . . H(ihi-1).
  31:      /// 
  32:      ///</summary>
  33:      public class DORGHR
  34:      {
  35:      
  36:   
  37:          #region Dependencies
  38:          
  39:          DORGQR _dorgqr; XERBLA _xerbla; ILAENV _ilaenv; 
  40:   
  41:          #endregion
  42:   
  43:   
  44:          #region Fields
  45:          
  46:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool LQUERY = false; int I = 0; int IINFO = 0; int J = 0; 
  47:          int LWKOPT = 0;int NB = 0; int NH = 0; 
  48:   
  49:          #endregion
  50:   
  51:          public DORGHR(DORGQR dorgqr, XERBLA xerbla, ILAENV ilaenv)
  52:          {
  53:      
  54:   
  55:              #region Set Dependencies
  56:              
  57:              this._dorgqr = dorgqr; this._xerbla = xerbla; this._ilaenv = ilaenv; 
  58:   
  59:              #endregion
  60:   
  61:          }
  62:      
  63:          public DORGHR()
  64:          {
  65:      
  66:   
  67:              #region Dependencies (Initialization)
  68:              
  69:              LSAME lsame = new LSAME();
  70:              DCOPY dcopy = new DCOPY();
  71:              XERBLA xerbla = new XERBLA();
  72:              DSCAL dscal = new DSCAL();
  73:              IEEECK ieeeck = new IEEECK();
  74:              IPARMQ iparmq = new IPARMQ();
  75:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  76:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  77:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
  78:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  79:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
  80:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
  81:              DGER dger = new DGER(xerbla);
  82:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  83:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
  84:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  85:              DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
  86:   
  87:              #endregion
  88:   
  89:   
  90:              #region Set Dependencies
  91:              
  92:              this._dorgqr = dorgqr; this._xerbla = xerbla; this._ilaenv = ilaenv; 
  93:   
  94:              #endregion
  95:   
  96:          }
  97:          /// <summary>
  98:          /// Purpose
  99:          /// =======
 100:          /// 
 101:          /// DORGHR generates a real orthogonal matrix Q which is defined as the
 102:          /// product of IHI-ILO elementary reflectors of order N, as returned by
 103:          /// DGEHRD:
 104:          /// 
 105:          /// Q = H(ilo) H(ilo+1) . . . H(ihi-1).
 106:          /// 
 107:          ///</summary>
 108:          /// <param name="N">
 109:          /// (input) INTEGER
 110:          /// The order of the matrix Q. N .GE. 0.
 111:          ///</param>
 112:          /// <param name="ILO">
 113:          /// (input) INTEGER
 114:          ///</param>
 115:          /// <param name="IHI">
 116:          /// (input) INTEGER
 117:          /// ILO and IHI must have the same values as in the previous call
 118:          /// of DGEHRD. Q is equal to the unit matrix except in the
 119:          /// submatrix Q(ilo+1:ihi,ilo+1:ihi).
 120:          /// 1 .LE. ILO .LE. IHI .LE. N, if N .GT. 0; ILO=1 and IHI=0, if N=0.
 121:          ///</param>
 122:          /// <param name="A">
 123:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 124:          /// On entry, the vectors which define the elementary reflectors,
 125:          /// as returned by DGEHRD.
 126:          /// On exit, the N-by-N orthogonal matrix Q.
 127:          ///</param>
 128:          /// <param name="LDA">
 129:          /// (input) INTEGER
 130:          /// The leading dimension of the array A. LDA .GE. max(1,N).
 131:          ///</param>
 132:          /// <param name="TAU">
 133:          /// (input) DOUBLE PRECISION array, dimension (N-1)
 134:          /// TAU(i) must contain the scalar factor of the elementary
 135:          /// reflector H(i), as returned by DGEHRD.
 136:          ///</param>
 137:          /// <param name="WORK">
 138:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 139:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 140:          ///</param>
 141:          /// <param name="LWORK">
 142:          /// (input) INTEGER
 143:          /// The dimension of the array WORK. LWORK .GE. IHI-ILO.
 144:          /// For optimum performance LWORK .GE. (IHI-ILO)*NB, where NB is
 145:          /// the optimal blocksize.
 146:          /// 
 147:          /// If LWORK = -1, then a workspace query is assumed; the routine
 148:          /// only calculates the optimal size of the WORK array, returns
 149:          /// this value as the first entry of the WORK array, and no error
 150:          /// message related to LWORK is issued by XERBLA.
 151:          ///</param>
 152:          /// <param name="INFO">
 153:          /// (output) INTEGER
 154:          /// = 0:  successful exit
 155:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 156:          ///</param>
 157:          public void Run(int N, int ILO, int IHI, ref double[] A, int offset_a, int LDA, double[] TAU, int offset_tau
 158:                           , ref double[] WORK, int offset_work, int LWORK, ref int INFO)
 159:          {
 160:   
 161:              #region Array Index Correction
 162:              
 163:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_work = -1 + offset_work; 
 164:   
 165:              #endregion
 166:   
 167:   
 168:              #region Prolog
 169:              
 170:              // *
 171:              // *  -- LAPACK routine (version 3.1) --
 172:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 173:              // *     November 2006
 174:              // *
 175:              // *     .. Scalar Arguments ..
 176:              // *     ..
 177:              // *     .. Array Arguments ..
 178:              // *     ..
 179:              // *
 180:              // *  Purpose
 181:              // *  =======
 182:              // *
 183:              // *  DORGHR generates a real orthogonal matrix Q which is defined as the
 184:              // *  product of IHI-ILO elementary reflectors of order N, as returned by
 185:              // *  DGEHRD:
 186:              // *
 187:              // *  Q = H(ilo) H(ilo+1) . . . H(ihi-1).
 188:              // *
 189:              // *  Arguments
 190:              // *  =========
 191:              // *
 192:              // *  N       (input) INTEGER
 193:              // *          The order of the matrix Q. N >= 0.
 194:              // *
 195:              // *  ILO     (input) INTEGER
 196:              // *  IHI     (input) INTEGER
 197:              // *          ILO and IHI must have the same values as in the previous call
 198:              // *          of DGEHRD. Q is equal to the unit matrix except in the
 199:              // *          submatrix Q(ilo+1:ihi,ilo+1:ihi).
 200:              // *          1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
 201:              // *
 202:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 203:              // *          On entry, the vectors which define the elementary reflectors,
 204:              // *          as returned by DGEHRD.
 205:              // *          On exit, the N-by-N orthogonal matrix Q.
 206:              // *
 207:              // *  LDA     (input) INTEGER
 208:              // *          The leading dimension of the array A. LDA >= max(1,N).
 209:              // *
 210:              // *  TAU     (input) DOUBLE PRECISION array, dimension (N-1)
 211:              // *          TAU(i) must contain the scalar factor of the elementary
 212:              // *          reflector H(i), as returned by DGEHRD.
 213:              // *
 214:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 215:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 216:              // *
 217:              // *  LWORK   (input) INTEGER
 218:              // *          The dimension of the array WORK. LWORK >= IHI-ILO.
 219:              // *          For optimum performance LWORK >= (IHI-ILO)*NB, where NB is
 220:              // *          the optimal blocksize.
 221:              // *
 222:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 223:              // *          only calculates the optimal size of the WORK array, returns
 224:              // *          this value as the first entry of the WORK array, and no error
 225:              // *          message related to LWORK is issued by XERBLA.
 226:              // *
 227:              // *  INFO    (output) INTEGER
 228:              // *          = 0:  successful exit
 229:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 230:              // *
 231:              // *  =====================================================================
 232:              // *
 233:              // *     .. Parameters ..
 234:              // *     ..
 235:              // *     .. Local Scalars ..
 236:              // *     ..
 237:              // *     .. External Subroutines ..
 238:              // *     ..
 239:              // *     .. External Functions ..
 240:              // *     ..
 241:              // *     .. Intrinsic Functions ..
 242:              //      INTRINSIC          MAX, MIN;
 243:              // *     ..
 244:              // *     .. Executable Statements ..
 245:              // *
 246:              // *     Test the input arguments
 247:              // *
 248:   
 249:              #endregion
 250:   
 251:   
 252:              #region Body
 253:              
 254:              INFO = 0;
 255:              NH = IHI - ILO;
 256:              LQUERY = (LWORK ==  - 1);
 257:              if (N < 0)
 258:              {
 259:                  INFO =  - 1;
 260:              }
 261:              else
 262:              {
 263:                  if (ILO < 1 || ILO > Math.Max(1, N))
 264:                  {
 265:                      INFO =  - 2;
 266:                  }
 267:                  else
 268:                  {
 269:                      if (IHI < Math.Min(ILO, N) || IHI > N)
 270:                      {
 271:                          INFO =  - 3;
 272:                      }
 273:                      else
 274:                      {
 275:                          if (LDA < Math.Max(1, N))
 276:                          {
 277:                              INFO =  - 5;
 278:                          }
 279:                          else
 280:                          {
 281:                              if (LWORK < Math.Max(1, NH) && !LQUERY)
 282:                              {
 283:                                  INFO =  - 8;
 284:                              }
 285:                          }
 286:                      }
 287:                  }
 288:              }
 289:              // *
 290:              if (INFO == 0)
 291:              {
 292:                  NB = this._ilaenv.Run(1, "DORGQR", " ", NH, NH, NH,  - 1);
 293:                  LWKOPT = Math.Max(1, NH) * NB;
 294:                  WORK[1 + o_work] = LWKOPT;
 295:              }
 296:              // *
 297:              if (INFO != 0)
 298:              {
 299:                  this._xerbla.Run("DORGHR",  - INFO);
 300:                  return;
 301:              }
 302:              else
 303:              {
 304:                  if (LQUERY)
 305:                  {
 306:                      return;
 307:                  }
 308:              }
 309:              // *
 310:              // *     Quick return if possible
 311:              // *
 312:              if (N == 0)
 313:              {
 314:                  WORK[1 + o_work] = 1;
 315:                  return;
 316:              }
 317:              // *
 318:              // *     Shift the vectors which define the elementary reflectors one
 319:              // *     column to the right, and set the first ilo and the last n-ihi
 320:              // *     rows and columns to those of the unit matrix
 321:              // *
 322:              for (J = IHI; J >= ILO + 1; J +=  - 1)
 323:              {
 324:                  for (I = 1; I <= J - 1; I++)
 325:                  {
 326:                      A[I+J * LDA + o_a] = ZERO;
 327:                  }
 328:                  for (I = J + 1; I <= IHI; I++)
 329:                  {
 330:                      A[I+J * LDA + o_a] = A[I+(J - 1) * LDA + o_a];
 331:                  }
 332:                  for (I = IHI + 1; I <= N; I++)
 333:                  {
 334:                      A[I+J * LDA + o_a] = ZERO;
 335:                  }
 336:              }
 337:              for (J = 1; J <= ILO; J++)
 338:              {
 339:                  for (I = 1; I <= N; I++)
 340:                  {
 341:                      A[I+J * LDA + o_a] = ZERO;
 342:                  }
 343:                  A[J+J * LDA + o_a] = ONE;
 344:              }
 345:              for (J = IHI + 1; J <= N; J++)
 346:              {
 347:                  for (I = 1; I <= N; I++)
 348:                  {
 349:                      A[I+J * LDA + o_a] = ZERO;
 350:                  }
 351:                  A[J+J * LDA + o_a] = ONE;
 352:              }
 353:              // *
 354:              if (NH > 0)
 355:              {
 356:                  // *
 357:                  // *        Generate Q(ilo+1:ihi,ilo+1:ihi)
 358:                  // *
 359:                  this._dorgqr.Run(NH, NH, NH, ref A, ILO + 1+(ILO + 1) * LDA + o_a, LDA, TAU, ILO + o_tau
 360:                                   , ref WORK, offset_work, LWORK, ref IINFO);
 361:              }
 362:              WORK[1 + o_work] = LWKOPT;
 363:              return;
 364:              // *
 365:              // *     End of DORGHR
 366:              // *
 367:   
 368:              #endregion
 369:   
 370:          }
 371:      }
 372:  }