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