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:      /// DORMLQ 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 defined as the product of k
  33:      /// elementary reflectors
  34:      /// 
  35:      /// Q = H(k) . . . H(2) H(1)
  36:      /// 
  37:      /// as returned by DGELQF. Q is of order M if SIDE = 'L' and of order N
  38:      /// if SIDE = 'R'.
  39:      /// 
  40:      ///</summary>
  41:      public class DORMLQ
  42:      {
  43:      
  44:   
  45:          #region Dependencies
  46:          
  47:          LSAME _lsame; ILAENV _ilaenv; DLARFB _dlarfb; DLARFT _dlarft; DORML2 _dorml2; XERBLA _xerbla; 
  48:   
  49:          #endregion
  50:   
  51:   
  52:          #region Fields
  53:          
  54:          const int NBMAX = 64; const int LDT = NBMAX + 1; bool LEFT = false; bool LQUERY = false; bool NOTRAN = false; 
  55:          string TRANST = new string(' ', 1);int I = 0; int I1 = 0; int I2 = 0; int I3 = 0; int IB = 0; int IC = 0; int IINFO = 0; 
  56:          int IWS = 0;int JC = 0; int LDWORK = 0; int LWKOPT = 0; int MI = 0; int NB = 0; int NBMIN = 0; int NI = 0; int NQ = 0; 
  57:          int NW = 0;double[] T = new double[LDT * NBMAX]; int offset_t = 0; 
  58:   
  59:          #endregion
  60:   
  61:          public DORMLQ(LSAME lsame, ILAENV ilaenv, DLARFB dlarfb, DLARFT dlarft, DORML2 dorml2, XERBLA xerbla)
  62:          {
  63:      
  64:   
  65:              #region Set Dependencies
  66:              
  67:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlarfb = dlarfb; this._dlarft = dlarft; this._dorml2 = dorml2; 
  68:              this._xerbla = xerbla;
  69:   
  70:              #endregion
  71:   
  72:          }
  73:      
  74:          public DORMLQ()
  75:          {
  76:      
  77:   
  78:              #region Dependencies (Initialization)
  79:              
  80:              LSAME lsame = new LSAME();
  81:              IEEECK ieeeck = new IEEECK();
  82:              IPARMQ iparmq = new IPARMQ();
  83:              DCOPY dcopy = new DCOPY();
  84:              XERBLA xerbla = new XERBLA();
  85:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  86:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  87:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  88:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
  89:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  90:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
  91:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
  92:              DGER dger = new DGER(xerbla);
  93:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
  94:              DORML2 dorml2 = new DORML2(lsame, dlarf, xerbla);
  95:   
  96:              #endregion
  97:   
  98:   
  99:              #region Set Dependencies
 100:              
 101:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlarfb = dlarfb; this._dlarft = dlarft; this._dorml2 = dorml2; 
 102:              this._xerbla = xerbla;
 103:   
 104:              #endregion
 105:   
 106:          }
 107:          /// <summary>
 108:          /// Purpose
 109:          /// =======
 110:          /// 
 111:          /// DORMLQ overwrites the general real M-by-N matrix C with
 112:          /// 
 113:          /// SIDE = 'L'     SIDE = 'R'
 114:          /// TRANS = 'N':      Q * C          C * Q
 115:          /// TRANS = 'T':      Q**T * C       C * Q**T
 116:          /// 
 117:          /// where Q is a real orthogonal matrix defined as the product of k
 118:          /// elementary reflectors
 119:          /// 
 120:          /// Q = H(k) . . . H(2) H(1)
 121:          /// 
 122:          /// as returned by DGELQF. Q is of order M if SIDE = 'L' and of order N
 123:          /// if SIDE = 'R'.
 124:          /// 
 125:          ///</summary>
 126:          /// <param name="SIDE">
 127:          /// = 'L'     SIDE = 'R'
 128:          ///</param>
 129:          /// <param name="TRANS">
 130:          /// (input) CHARACTER*1
 131:          /// = 'N':  No transpose, apply Q;
 132:          /// = 'T':  Transpose, apply Q**T.
 133:          ///</param>
 134:          /// <param name="M">
 135:          /// (input) INTEGER
 136:          /// The number of rows of the matrix C. M .GE. 0.
 137:          ///</param>
 138:          /// <param name="N">
 139:          /// (input) INTEGER
 140:          /// The number of columns of the matrix C. N .GE. 0.
 141:          ///</param>
 142:          /// <param name="K">
 143:          /// (input) INTEGER
 144:          /// The number of elementary reflectors whose product defines
 145:          /// the matrix Q.
 146:          /// If SIDE = 'L', M .GE. K .GE. 0;
 147:          /// if SIDE = 'R', N .GE. K .GE. 0.
 148:          ///</param>
 149:          /// <param name="A">
 150:          /// (input) DOUBLE PRECISION array, dimension
 151:          /// (LDA,M) if SIDE = 'L',
 152:          /// (LDA,N) if SIDE = 'R'
 153:          /// The i-th row must contain the vector which defines the
 154:          /// elementary reflector H(i), for i = 1,2,...,k, as returned by
 155:          /// DGELQF in the first k rows of its array argument A.
 156:          /// A is modified by the routine but restored on exit.
 157:          ///</param>
 158:          /// <param name="LDA">
 159:          /// (input) INTEGER
 160:          /// The leading dimension of the array A. LDA .GE. max(1,K).
 161:          ///</param>
 162:          /// <param name="TAU">
 163:          /// (input) DOUBLE PRECISION array, dimension (K)
 164:          /// TAU(i) must contain the scalar factor of the elementary
 165:          /// reflector H(i), as returned by DGELQF.
 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 TRANS, int M, int N, int K, 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);  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:              // *  DORMLQ 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 defined as the product of k
 241:              // *  elementary reflectors
 242:              // *
 243:              // *        Q = H(k) . . . H(2) H(1)
 244:              // *
 245:              // *  as returned by DGELQF. Q is of order M if SIDE = 'L' and of order N
 246:              // *  if SIDE = 'R'.
 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:              // *  TRANS   (input) CHARACTER*1
 256:              // *          = 'N':  No transpose, apply Q;
 257:              // *          = 'T':  Transpose, apply Q**T.
 258:              // *
 259:              // *  M       (input) INTEGER
 260:              // *          The number of rows of the matrix C. M >= 0.
 261:              // *
 262:              // *  N       (input) INTEGER
 263:              // *          The number of columns of the matrix C. N >= 0.
 264:              // *
 265:              // *  K       (input) INTEGER
 266:              // *          The number of elementary reflectors whose product defines
 267:              // *          the matrix Q.
 268:              // *          If SIDE = 'L', M >= K >= 0;
 269:              // *          if SIDE = 'R', N >= K >= 0.
 270:              // *
 271:              // *  A       (input) DOUBLE PRECISION array, dimension
 272:              // *                               (LDA,M) if SIDE = 'L',
 273:              // *                               (LDA,N) if SIDE = 'R'
 274:              // *          The i-th row must contain the vector which defines the
 275:              // *          elementary reflector H(i), for i = 1,2,...,k, as returned by
 276:              // *          DGELQF in the first k rows of its array argument A.
 277:              // *          A is modified by the routine but restored on exit.
 278:              // *
 279:              // *  LDA     (input) INTEGER
 280:              // *          The leading dimension of the array A. LDA >= max(1,K).
 281:              // *
 282:              // *  TAU     (input) DOUBLE PRECISION array, dimension (K)
 283:              // *          TAU(i) must contain the scalar factor of the elementary
 284:              // *          reflector H(i), as returned by DGELQF.
 285:              // *
 286:              // *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 287:              // *          On entry, the M-by-N matrix C.
 288:              // *          On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q.
 289:              // *
 290:              // *  LDC     (input) INTEGER
 291:              // *          The leading dimension of the array C. LDC >= max(1,M).
 292:              // *
 293:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 294:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 295:              // *
 296:              // *  LWORK   (input) INTEGER
 297:              // *          The dimension of the array WORK.
 298:              // *          If SIDE = 'L', LWORK >= max(1,N);
 299:              // *          if SIDE = 'R', LWORK >= max(1,M).
 300:              // *          For optimum performance LWORK >= N*NB if SIDE = 'L', and
 301:              // *          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
 302:              // *          blocksize.
 303:              // *
 304:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 305:              // *          only calculates the optimal size of the WORK array, returns
 306:              // *          this value as the first entry of the WORK array, and no error
 307:              // *          message related to LWORK is issued by XERBLA.
 308:              // *
 309:              // *  INFO    (output) INTEGER
 310:              // *          = 0:  successful exit
 311:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 312:              // *
 313:              // *  =====================================================================
 314:              // *
 315:              // *     .. Parameters ..
 316:              // *     ..
 317:              // *     .. Local Scalars ..
 318:              // *     ..
 319:              // *     .. Local Arrays ..
 320:              // *     ..
 321:              // *     .. External Functions ..
 322:              // *     ..
 323:              // *     .. External Subroutines ..
 324:              // *     ..
 325:              // *     .. Intrinsic Functions ..
 326:              //      INTRINSIC          MAX, MIN;
 327:              // *     ..
 328:              // *     .. Executable Statements ..
 329:              // *
 330:              // *     Test the input arguments
 331:              // *
 332:   
 333:              #endregion
 334:   
 335:   
 336:              #region Body
 337:              
 338:              INFO = 0;
 339:              LEFT = this._lsame.Run(SIDE, "L");
 340:              NOTRAN = this._lsame.Run(TRANS, "N");
 341:              LQUERY = (LWORK ==  - 1);
 342:              // *
 343:              // *     NQ is the order of Q and NW is the minimum dimension of WORK
 344:              // *
 345:              if (LEFT)
 346:              {
 347:                  NQ = M;
 348:                  NW = N;
 349:              }
 350:              else
 351:              {
 352:                  NQ = N;
 353:                  NW = M;
 354:              }
 355:              if (!LEFT && !this._lsame.Run(SIDE, "R"))
 356:              {
 357:                  INFO =  - 1;
 358:              }
 359:              else
 360:              {
 361:                  if (!NOTRAN && !this._lsame.Run(TRANS, "T"))
 362:                  {
 363:                      INFO =  - 2;
 364:                  }
 365:                  else
 366:                  {
 367:                      if (M < 0)
 368:                      {
 369:                          INFO =  - 3;
 370:                      }
 371:                      else
 372:                      {
 373:                          if (N < 0)
 374:                          {
 375:                              INFO =  - 4;
 376:                          }
 377:                          else
 378:                          {
 379:                              if (K < 0 || K > NQ)
 380:                              {
 381:                                  INFO =  - 5;
 382:                              }
 383:                              else
 384:                              {
 385:                                  if (LDA < Math.Max(1, K))
 386:                                  {
 387:                                      INFO =  - 7;
 388:                                  }
 389:                                  else
 390:                                  {
 391:                                      if (LDC < Math.Max(1, M))
 392:                                      {
 393:                                          INFO =  - 10;
 394:                                      }
 395:                                      else
 396:                                      {
 397:                                          if (LWORK < Math.Max(1, NW) && !LQUERY)
 398:                                          {
 399:                                              INFO =  - 12;
 400:                                          }
 401:                                      }
 402:                                  }
 403:                              }
 404:                          }
 405:                      }
 406:                  }
 407:              }
 408:              // *
 409:              if (INFO == 0)
 410:              {
 411:                  // *
 412:                  // *        Determine the block size.  NB may be at most NBMAX, where NBMAX
 413:                  // *        is used to define the local array T.
 414:                  // *
 415:                  NB = Math.Min(NBMAX, this._ilaenv.Run(1, "DORMLQ", SIDE + TRANS, M, N, K,  - 1));
 416:                  LWKOPT = Math.Max(1, NW) * NB;
 417:                  WORK[1 + o_work] = LWKOPT;
 418:              }
 419:              // *
 420:              if (INFO != 0)
 421:              {
 422:                  this._xerbla.Run("DORMLQ",  - INFO);
 423:                  return;
 424:              }
 425:              else
 426:              {
 427:                  if (LQUERY)
 428:                  {
 429:                      return;
 430:                  }
 431:              }
 432:              // *
 433:              // *     Quick return if possible
 434:              // *
 435:              if (M == 0 || N == 0 || K == 0)
 436:              {
 437:                  WORK[1 + o_work] = 1;
 438:                  return;
 439:              }
 440:              // *
 441:              NBMIN = 2;
 442:              LDWORK = NW;
 443:              if (NB > 1 && NB < K)
 444:              {
 445:                  IWS = NW * NB;
 446:                  if (LWORK < IWS)
 447:                  {
 448:                      NB = LWORK / LDWORK;
 449:                      NBMIN = Math.Max(2, this._ilaenv.Run(2, "DORMLQ", SIDE + TRANS, M, N, K,  - 1));
 450:                  }
 451:              }
 452:              else
 453:              {
 454:                  IWS = NW;
 455:              }
 456:              // *
 457:              if (NB < NBMIN || NB >= K)
 458:              {
 459:                  // *
 460:                  // *        Use unblocked code
 461:                  // *
 462:                  this._dorml2.Run(SIDE, TRANS, M, N, K, ref A, offset_a
 463:                                   , LDA, TAU, offset_tau, ref C, offset_c, LDC, ref WORK, offset_work, ref IINFO);
 464:              }
 465:              else
 466:              {
 467:                  // *
 468:                  // *        Use blocked code
 469:                  // *
 470:                  if ((LEFT && NOTRAN) || (!LEFT && !NOTRAN))
 471:                  {
 472:                      I1 = 1;
 473:                      I2 = K;
 474:                      I3 = NB;
 475:                  }
 476:                  else
 477:                  {
 478:                      I1 = ((K - 1) / NB) * NB + 1;
 479:                      I2 = 1;
 480:                      I3 =  - NB;
 481:                  }
 482:                  // *
 483:                  if (LEFT)
 484:                  {
 485:                      NI = N;
 486:                      JC = 1;
 487:                  }
 488:                  else
 489:                  {
 490:                      MI = M;
 491:                      IC = 1;
 492:                  }
 493:                  // *
 494:                  if (NOTRAN)
 495:                  {
 496:                      FortranLib.Copy(ref TRANST , "T");
 497:                  }
 498:                  else
 499:                  {
 500:                      FortranLib.Copy(ref TRANST , "N");
 501:                  }
 502:                  // *
 503:                  for (I = I1; (I3 >= 0) ? (I <= I2) : (I >= I2); I += I3)
 504:                  {
 505:                      IB = Math.Min(NB, K - I + 1);
 506:                      // *
 507:                      // *           Form the triangular factor of the block reflector
 508:                      // *           H = H(i) H(i+1) . . . H(i+ib-1)
 509:                      // *
 510:                      this._dlarft.Run("Forward", "Rowwise", NQ - I + 1, IB, ref A, I+I * LDA + o_a, LDA
 511:                                       , TAU, I + o_tau, ref T, offset_t, LDT);
 512:                      if (LEFT)
 513:                      {
 514:                          // *
 515:                          // *              H or H' is applied to C(i:m,1:n)
 516:                          // *
 517:                          MI = M - I + 1;
 518:                          IC = I;
 519:                      }
 520:                      else
 521:                      {
 522:                          // *
 523:                          // *              H or H' is applied to C(1:m,i:n)
 524:                          // *
 525:                          NI = N - I + 1;
 526:                          JC = I;
 527:                      }
 528:                      // *
 529:                      // *           Apply H or H'
 530:                      // *
 531:                      this._dlarfb.Run(SIDE, TRANST, "Forward", "Rowwise", MI, NI
 532:                                       , IB, A, I+I * LDA + o_a, LDA, T, offset_t, LDT, ref C, IC+JC * LDC + o_c
 533:                                       , LDC, ref WORK, offset_work, LDWORK);
 534:                  }
 535:              }
 536:              WORK[1 + o_work] = LWKOPT;
 537:              return;
 538:              // *
 539:              // *     End of DORMLQ
 540:              // *
 541:   
 542:              #endregion
 543:   
 544:          }
 545:      }
 546:  }