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:      /// DORMTR overwrites the general real M-by-N matrix C with
  27:      /// 
  28:      /// SIDE = 'L'     SIDE = 'R'
  29:      /// TRANS = 'N':      Q * C          C * Q
  30:      /// TRANS = 'T':      Q**T * C       C * Q**T
  31:      /// 
  32:      /// where Q is a real orthogonal matrix of order nq, with nq = m if
  33:      /// SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
  34:      /// nq-1 elementary reflectors, as returned by DSYTRD:
  35:      /// 
  36:      /// if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
  37:      /// 
  38:      /// if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
  39:      /// 
  40:      ///</summary>
  41:      public class DORMTR
  42:      {
  43:      
  44:   
  45:          #region Dependencies
  46:          
  47:          LSAME _lsame; ILAENV _ilaenv; DORMQL _dormql; DORMQR _dormqr; XERBLA _xerbla; 
  48:   
  49:          #endregion
  50:   
  51:   
  52:          #region Fields
  53:          
  54:          bool LEFT = false; bool LQUERY = false; bool UPPER = false; int I1 = 0; int I2 = 0; int IINFO = 0; int LWKOPT = 0; 
  55:          int MI = 0;int NB = 0; int NI = 0; int NQ = 0; int NW = 0; 
  56:   
  57:          #endregion
  58:   
  59:          public DORMTR(LSAME lsame, ILAENV ilaenv, DORMQL dormql, DORMQR dormqr, XERBLA xerbla)
  60:          {
  61:      
  62:   
  63:              #region Set Dependencies
  64:              
  65:              this._lsame = lsame; this._ilaenv = ilaenv; this._dormql = dormql; this._dormqr = dormqr; this._xerbla = xerbla; 
  66:   
  67:              #endregion
  68:   
  69:          }
  70:      
  71:          public DORMTR()
  72:          {
  73:      
  74:   
  75:              #region Dependencies (Initialization)
  76:              
  77:              LSAME lsame = new LSAME();
  78:              IEEECK ieeeck = new IEEECK();
  79:              IPARMQ iparmq = new IPARMQ();
  80:              DCOPY dcopy = new DCOPY();
  81:              XERBLA xerbla = new XERBLA();
  82:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  83:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  84:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  85:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
  86:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  87:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
  88:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
  89:              DGER dger = new DGER(xerbla);
  90:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  91:              DORM2L dorm2l = new DORM2L(lsame, dlarf, xerbla);
  92:              DORMQL dormql = new DORMQL(lsame, ilaenv, dlarfb, dlarft, dorm2l, xerbla);
  93:              DORM2R dorm2r = new DORM2R(lsame, dlarf, xerbla);
  94:              DORMQR dormqr = new DORMQR(lsame, ilaenv, dlarfb, dlarft, dorm2r, xerbla);
  95:   
  96:              #endregion
  97:   
  98:   
  99:              #region Set Dependencies
 100:              
 101:              this._lsame = lsame; this._ilaenv = ilaenv; this._dormql = dormql; this._dormqr = dormqr; this._xerbla = xerbla; 
 102:   
 103:              #endregion
 104:   
 105:          }
 106:          /// <summary>
 107:          /// Purpose
 108:          /// =======
 109:          /// 
 110:          /// DORMTR overwrites the general real M-by-N matrix C with
 111:          /// 
 112:          /// SIDE = 'L'     SIDE = 'R'
 113:          /// TRANS = 'N':      Q * C          C * Q
 114:          /// TRANS = 'T':      Q**T * C       C * Q**T
 115:          /// 
 116:          /// where Q is a real orthogonal matrix of order nq, with nq = m if
 117:          /// SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
 118:          /// nq-1 elementary reflectors, as returned by DSYTRD:
 119:          /// 
 120:          /// if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
 121:          /// 
 122:          /// if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
 123:          /// 
 124:          ///</summary>
 125:          /// <param name="SIDE">
 126:          /// = 'L'     SIDE = 'R'
 127:          ///</param>
 128:          /// <param name="UPLO">
 129:          /// (input) CHARACTER*1
 130:          /// = 'U': Upper triangle of A contains elementary reflectors
 131:          /// from DSYTRD;
 132:          /// = 'L': Lower triangle of A contains elementary reflectors
 133:          /// from DSYTRD.
 134:          ///</param>
 135:          /// <param name="TRANS">
 136:          /// (input) CHARACTER*1
 137:          /// = 'N':  No transpose, apply Q;
 138:          /// = 'T':  Transpose, apply Q**T.
 139:          ///</param>
 140:          /// <param name="M">
 141:          /// (input) INTEGER
 142:          /// The number of rows of the matrix C. M .GE. 0.
 143:          ///</param>
 144:          /// <param name="N">
 145:          /// (input) INTEGER
 146:          /// The number of columns of the matrix C. N .GE. 0.
 147:          ///</param>
 148:          /// <param name="A">
 149:          /// (input) DOUBLE PRECISION array, dimension
 150:          /// (LDA,M) if SIDE = 'L'
 151:          /// (LDA,N) if SIDE = 'R'
 152:          /// The vectors which define the elementary reflectors, as
 153:          /// returned by DSYTRD.
 154:          ///</param>
 155:          /// <param name="LDA">
 156:          /// (input) INTEGER
 157:          /// The leading dimension of the array A.
 158:          /// LDA .GE. max(1,M) if SIDE = 'L'; LDA .GE. max(1,N) if SIDE = 'R'.
 159:          ///</param>
 160:          /// <param name="TAU">
 161:          /// (input) DOUBLE PRECISION array, dimension
 162:          /// (M-1) if SIDE = 'L'
 163:          /// (N-1) if SIDE = 'R'
 164:          /// TAU(i) must contain the scalar factor of the elementary
 165:          /// reflector H(i), as returned by DSYTRD.
 166:          ///</param>
 167:          /// <param name="C">
 168:          /// (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 169:          /// On entry, the M-by-N matrix C.
 170:          /// On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
 171:          ///</param>
 172:          /// <param name="LDC">
 173:          /// (input) INTEGER
 174:          /// The leading dimension of the array C. LDC .GE. max(1,M).
 175:          ///</param>
 176:          /// <param name="WORK">
 177:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 178:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 179:          ///</param>
 180:          /// <param name="LWORK">
 181:          /// (input) INTEGER
 182:          /// The dimension of the array WORK.
 183:          /// If SIDE = 'L', LWORK .GE. max(1,N);
 184:          /// if SIDE = 'R', LWORK .GE. max(1,M).
 185:          /// For optimum performance LWORK .GE. N*NB if SIDE = 'L', and
 186:          /// LWORK .GE. M*NB if SIDE = 'R', where NB is the optimal
 187:          /// blocksize.
 188:          /// 
 189:          /// If LWORK = -1, then a workspace query is assumed; the routine
 190:          /// only calculates the optimal size of the WORK array, returns
 191:          /// this value as the first entry of the WORK array, and no error
 192:          /// message related to LWORK is issued by XERBLA.
 193:          ///</param>
 194:          /// <param name="INFO">
 195:          /// (output) INTEGER
 196:          /// = 0:  successful exit
 197:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 198:          ///</param>
 199:          public void Run(string SIDE, string UPLO, string TRANS, int M, int N, ref double[] A, int offset_a
 200:                           , int LDA, double[] TAU, int offset_tau, ref double[] C, int offset_c, int LDC, ref double[] WORK, int offset_work, int LWORK
 201:                           , ref int INFO)
 202:          {
 203:   
 204:              #region Array Index Correction
 205:              
 206:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_c = -1 - LDC + offset_c; 
 207:               int o_work = -1 + offset_work;
 208:   
 209:              #endregion
 210:   
 211:   
 212:              #region Strings
 213:              
 214:              SIDE = SIDE.Substring(0, 1);  UPLO = UPLO.Substring(0, 1);  TRANS = TRANS.Substring(0, 1);  
 215:   
 216:              #endregion
 217:   
 218:   
 219:              #region Prolog
 220:              
 221:              // *
 222:              // *  -- LAPACK routine (version 3.1) --
 223:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 224:              // *     November 2006
 225:              // *
 226:              // *     .. Scalar Arguments ..
 227:              // *     ..
 228:              // *     .. Array Arguments ..
 229:              // *     ..
 230:              // *
 231:              // *  Purpose
 232:              // *  =======
 233:              // *
 234:              // *  DORMTR overwrites the general real M-by-N matrix C with
 235:              // *
 236:              // *                  SIDE = 'L'     SIDE = 'R'
 237:              // *  TRANS = 'N':      Q * C          C * Q
 238:              // *  TRANS = 'T':      Q**T * C       C * Q**T
 239:              // *
 240:              // *  where Q is a real orthogonal matrix of order nq, with nq = m if
 241:              // *  SIDE = 'L' and nq = n if SIDE = 'R'. Q is defined as the product of
 242:              // *  nq-1 elementary reflectors, as returned by DSYTRD:
 243:              // *
 244:              // *  if UPLO = 'U', Q = H(nq-1) . . . H(2) H(1);
 245:              // *
 246:              // *  if UPLO = 'L', Q = H(1) H(2) . . . H(nq-1).
 247:              // *
 248:              // *  Arguments
 249:              // *  =========
 250:              // *
 251:              // *  SIDE    (input) CHARACTER*1
 252:              // *          = 'L': apply Q or Q**T from the Left;
 253:              // *          = 'R': apply Q or Q**T from the Right.
 254:              // *
 255:              // *  UPLO    (input) CHARACTER*1
 256:              // *          = 'U': Upper triangle of A contains elementary reflectors
 257:              // *                 from DSYTRD;
 258:              // *          = 'L': Lower triangle of A contains elementary reflectors
 259:              // *                 from DSYTRD.
 260:              // *
 261:              // *  TRANS   (input) CHARACTER*1
 262:              // *          = 'N':  No transpose, apply Q;
 263:              // *          = 'T':  Transpose, apply Q**T.
 264:              // *
 265:              // *  M       (input) INTEGER
 266:              // *          The number of rows of the matrix C. M >= 0.
 267:              // *
 268:              // *  N       (input) INTEGER
 269:              // *          The number of columns of the matrix C. N >= 0.
 270:              // *
 271:              // *  A       (input) DOUBLE PRECISION array, dimension
 272:              // *                               (LDA,M) if SIDE = 'L'
 273:              // *                               (LDA,N) if SIDE = 'R'
 274:              // *          The vectors which define the elementary reflectors, as
 275:              // *          returned by DSYTRD.
 276:              // *
 277:              // *  LDA     (input) INTEGER
 278:              // *          The leading dimension of the array A.
 279:              // *          LDA >= max(1,M) if SIDE = 'L'; LDA >= max(1,N) if SIDE = 'R'.
 280:              // *
 281:              // *  TAU     (input) DOUBLE PRECISION array, dimension
 282:              // *                               (M-1) if SIDE = 'L'
 283:              // *                               (N-1) if SIDE = 'R'
 284:              // *          TAU(i) must contain the scalar factor of the elementary
 285:              // *          reflector H(i), as returned by DSYTRD.
 286:              // *
 287:              // *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 288:              // *          On entry, the M-by-N matrix C.
 289:              // *          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
 290:              // *
 291:              // *  LDC     (input) INTEGER
 292:              // *          The leading dimension of the array C. LDC >= max(1,M).
 293:              // *
 294:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 295:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 296:              // *
 297:              // *  LWORK   (input) INTEGER
 298:              // *          The dimension of the array WORK.
 299:              // *          If SIDE = 'L', LWORK >= max(1,N);
 300:              // *          if SIDE = 'R', LWORK >= max(1,M).
 301:              // *          For optimum performance LWORK >= N*NB if SIDE = 'L', and
 302:              // *          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
 303:              // *          blocksize.
 304:              // *
 305:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 306:              // *          only calculates the optimal size of the WORK array, returns
 307:              // *          this value as the first entry of the WORK array, and no error
 308:              // *          message related to LWORK is issued by XERBLA.
 309:              // *
 310:              // *  INFO    (output) INTEGER
 311:              // *          = 0:  successful exit
 312:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 313:              // *
 314:              // *  =====================================================================
 315:              // *
 316:              // *     .. Local Scalars ..
 317:              // *     ..
 318:              // *     .. External Functions ..
 319:              // *     ..
 320:              // *     .. External Subroutines ..
 321:              // *     ..
 322:              // *     .. Intrinsic Functions ..
 323:              //      INTRINSIC          MAX;
 324:              // *     ..
 325:              // *     .. Executable Statements ..
 326:              // *
 327:              // *     Test the input arguments
 328:              // *
 329:   
 330:              #endregion
 331:   
 332:   
 333:              #region Body
 334:              
 335:              INFO = 0;
 336:              LEFT = this._lsame.Run(SIDE, "L");
 337:              UPPER = this._lsame.Run(UPLO, "U");
 338:              LQUERY = (LWORK ==  - 1);
 339:              // *
 340:              // *     NQ is the order of Q and NW is the minimum dimension of WORK
 341:              // *
 342:              if (LEFT)
 343:              {
 344:                  NQ = M;
 345:                  NW = N;
 346:              }
 347:              else
 348:              {
 349:                  NQ = N;
 350:                  NW = M;
 351:              }
 352:              if (!LEFT && !this._lsame.Run(SIDE, "R"))
 353:              {
 354:                  INFO =  - 1;
 355:              }
 356:              else
 357:              {
 358:                  if (!UPPER && !this._lsame.Run(UPLO, "L"))
 359:                  {
 360:                      INFO =  - 2;
 361:                  }
 362:                  else
 363:                  {
 364:                      if (!this._lsame.Run(TRANS, "N") && !this._lsame.Run(TRANS, "T"))
 365:                      {
 366:                          INFO =  - 3;
 367:                      }
 368:                      else
 369:                      {
 370:                          if (M < 0)
 371:                          {
 372:                              INFO =  - 4;
 373:                          }
 374:                          else
 375:                          {
 376:                              if (N < 0)
 377:                              {
 378:                                  INFO =  - 5;
 379:                              }
 380:                              else
 381:                              {
 382:                                  if (LDA < Math.Max(1, NQ))
 383:                                  {
 384:                                      INFO =  - 7;
 385:                                  }
 386:                                  else
 387:                                  {
 388:                                      if (LDC < Math.Max(1, M))
 389:                                      {
 390:                                          INFO =  - 10;
 391:                                      }
 392:                                      else
 393:                                      {
 394:                                          if (LWORK < Math.Max(1, NW) && !LQUERY)
 395:                                          {
 396:                                              INFO =  - 12;
 397:                                          }
 398:                                      }
 399:                                  }
 400:                              }
 401:                          }
 402:                      }
 403:                  }
 404:              }
 405:              // *
 406:              if (INFO == 0)
 407:              {
 408:                  if (UPPER)
 409:                  {
 410:                      if (LEFT)
 411:                      {
 412:                          NB = this._ilaenv.Run(1, "DORMQL", SIDE + TRANS, M - 1, N, M - 1,  - 1);
 413:                      }
 414:                      else
 415:                      {
 416:                          NB = this._ilaenv.Run(1, "DORMQL", SIDE + TRANS, M, N - 1, N - 1,  - 1);
 417:                      }
 418:                  }
 419:                  else
 420:                  {
 421:                      if (LEFT)
 422:                      {
 423:                          NB = this._ilaenv.Run(1, "DORMQR", SIDE + TRANS, M - 1, N, M - 1,  - 1);
 424:                      }
 425:                      else
 426:                      {
 427:                          NB = this._ilaenv.Run(1, "DORMQR", SIDE + TRANS, M, N - 1, N - 1,  - 1);
 428:                      }
 429:                  }
 430:                  LWKOPT = Math.Max(1, NW) * NB;
 431:                  WORK[1 + o_work] = LWKOPT;
 432:              }
 433:              // *
 434:              if (INFO != 0)
 435:              {
 436:                  this._xerbla.Run("DORMTR",  - INFO);
 437:                  return;
 438:              }
 439:              else
 440:              {
 441:                  if (LQUERY)
 442:                  {
 443:                      return;
 444:                  }
 445:              }
 446:              // *
 447:              // *     Quick return if possible
 448:              // *
 449:              if (M == 0 || N == 0 || NQ == 1)
 450:              {
 451:                  WORK[1 + o_work] = 1;
 452:                  return;
 453:              }
 454:              // *
 455:              if (LEFT)
 456:              {
 457:                  MI = M - 1;
 458:                  NI = N;
 459:              }
 460:              else
 461:              {
 462:                  MI = M;
 463:                  NI = N - 1;
 464:              }
 465:              // *
 466:              if (UPPER)
 467:              {
 468:                  // *
 469:                  // *        Q was determined by a call to DSYTRD with UPLO = 'U'
 470:                  // *
 471:                  this._dormql.Run(SIDE, TRANS, MI, NI, NQ - 1, ref A, 1+2 * LDA + o_a
 472:                                   , LDA, TAU, offset_tau, ref C, offset_c, LDC, ref WORK, offset_work, LWORK
 473:                                   , ref IINFO);
 474:              }
 475:              else
 476:              {
 477:                  // *
 478:                  // *        Q was determined by a call to DSYTRD with UPLO = 'L'
 479:                  // *
 480:                  if (LEFT)
 481:                  {
 482:                      I1 = 2;
 483:                      I2 = 1;
 484:                  }
 485:                  else
 486:                  {
 487:                      I1 = 1;
 488:                      I2 = 2;
 489:                  }
 490:                  this._dormqr.Run(SIDE, TRANS, MI, NI, NQ - 1, ref A, 2+1 * LDA + o_a
 491:                                   , LDA, TAU, offset_tau, ref C, I1+I2 * LDC + o_c, LDC, ref WORK, offset_work, LWORK
 492:                                   , ref IINFO);
 493:              }
 494:              WORK[1 + o_work] = LWKOPT;
 495:              return;
 496:              // *
 497:              // *     End of DORMTR
 498:              // *
 499:   
 500:              #endregion
 501:   
 502:          }
 503:      }
 504:  }