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.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// January 2007
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DORMRZ 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(1) H(2) . . . H(k)
  36:      /// 
  37:      /// as returned by DTZRZF. Q is of order M if SIDE = 'L' and of order N
  38:      /// if SIDE = 'R'.
  39:      /// 
  40:      ///</summary>
  41:      public class DORMRZ
  42:      {
  43:      
  44:   
  45:          #region Dependencies
  46:          
  47:          LSAME _lsame; ILAENV _ilaenv; DLARZB _dlarzb; DLARZT _dlarzt; DORMR3 _dormr3; 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 JA = 0; int JC = 0; int LDWORK = 0; int LWKOPT = 0; int MI = 0; int NB = 0; int NBMIN = 0; int NI = 0; 
  57:          int NQ = 0;int NW = 0; double[] T = new double[LDT * NBMAX]; int offset_t = 0;
  58:   
  59:          #endregion
  60:   
  61:          public DORMRZ(LSAME lsame, ILAENV ilaenv, DLARZB dlarzb, DLARZT dlarzt, DORMR3 dormr3, XERBLA xerbla)
  62:          {
  63:      
  64:   
  65:              #region Set Dependencies
  66:              
  67:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlarzb = dlarzb; this._dlarzt = dlarzt; this._dormr3 = dormr3; 
  68:              this._xerbla = xerbla;
  69:   
  70:              #endregion
  71:   
  72:          }
  73:      
  74:          public DORMRZ()
  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:              DAXPY daxpy = new DAXPY();
  86:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  87:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  88:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  89:              DLARZB dlarzb = new DLARZB(lsame, dcopy, dgemm, dtrmm, xerbla);
  90:              DGEMV dgemv = new DGEMV(lsame, xerbla);
  91:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
  92:              DLARZT dlarzt = new DLARZT(dgemv, dtrmv, xerbla, lsame);
  93:              DGER dger = new DGER(xerbla);
  94:              DLARZ dlarz = new DLARZ(daxpy, dcopy, dgemv, dger, lsame);
  95:              DORMR3 dormr3 = new DORMR3(lsame, dlarz, xerbla);
  96:   
  97:              #endregion
  98:   
  99:   
 100:              #region Set Dependencies
 101:              
 102:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlarzb = dlarzb; this._dlarzt = dlarzt; this._dormr3 = dormr3; 
 103:              this._xerbla = xerbla;
 104:   
 105:              #endregion
 106:   
 107:          }
 108:          /// <summary>
 109:          /// Purpose
 110:          /// =======
 111:          /// 
 112:          /// DORMRZ overwrites the general real M-by-N matrix C with
 113:          /// 
 114:          /// SIDE = 'L'     SIDE = 'R'
 115:          /// TRANS = 'N':      Q * C          C * Q
 116:          /// TRANS = 'T':      Q**T * C       C * Q**T
 117:          /// 
 118:          /// where Q is a real orthogonal matrix defined as the product of k
 119:          /// elementary reflectors
 120:          /// 
 121:          /// Q = H(1) H(2) . . . H(k)
 122:          /// 
 123:          /// as returned by DTZRZF. Q is of order M if SIDE = 'L' and of order N
 124:          /// if SIDE = 'R'.
 125:          /// 
 126:          ///</summary>
 127:          /// <param name="SIDE">
 128:          /// = 'L'     SIDE = 'R'
 129:          ///</param>
 130:          /// <param name="TRANS">
 131:          /// (input) CHARACTER*1
 132:          /// = 'N':  No transpose, apply Q;
 133:          /// = 'T':  Transpose, apply Q**T.
 134:          ///</param>
 135:          /// <param name="M">
 136:          /// (input) INTEGER
 137:          /// The number of rows of the matrix C. M .GE. 0.
 138:          ///</param>
 139:          /// <param name="N">
 140:          /// (input) INTEGER
 141:          /// The number of columns of the matrix C. N .GE. 0.
 142:          ///</param>
 143:          /// <param name="K">
 144:          /// (input) INTEGER
 145:          /// The number of elementary reflectors whose product defines
 146:          /// the matrix Q.
 147:          /// If SIDE = 'L', M .GE. K .GE. 0;
 148:          /// if SIDE = 'R', N .GE. K .GE. 0.
 149:          ///</param>
 150:          /// <param name="L">
 151:          /// (input) INTEGER
 152:          /// The number of columns of the matrix A containing
 153:          /// the meaningful part of the Householder reflectors.
 154:          /// If SIDE = 'L', M .GE. L .GE. 0, if SIDE = 'R', N .GE. L .GE. 0.
 155:          ///</param>
 156:          /// <param name="A">
 157:          /// (input) DOUBLE PRECISION array, dimension
 158:          /// (LDA,M) if SIDE = 'L',
 159:          /// (LDA,N) if SIDE = 'R'
 160:          /// The i-th row must contain the vector which defines the
 161:          /// elementary reflector H(i), for i = 1,2,...,k, as returned by
 162:          /// DTZRZF in the last k rows of its array argument A.
 163:          /// A is modified by the routine but restored on exit.
 164:          ///</param>
 165:          /// <param name="LDA">
 166:          /// (input) INTEGER
 167:          /// The leading dimension of the array A. LDA .GE. max(1,K).
 168:          ///</param>
 169:          /// <param name="TAU">
 170:          /// (input) DOUBLE PRECISION array, dimension (K)
 171:          /// TAU(i) must contain the scalar factor of the elementary
 172:          /// reflector H(i), as returned by DTZRZF.
 173:          ///</param>
 174:          /// <param name="C">
 175:          /// (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 176:          /// On entry, the M-by-N matrix C.
 177:          /// On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
 178:          ///</param>
 179:          /// <param name="LDC">
 180:          /// (input) INTEGER
 181:          /// The leading dimension of the array C. LDC .GE. max(1,M).
 182:          ///</param>
 183:          /// <param name="WORK">
 184:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 185:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 186:          ///</param>
 187:          /// <param name="LWORK">
 188:          /// (input) INTEGER
 189:          /// The dimension of the array WORK.
 190:          /// If SIDE = 'L', LWORK .GE. max(1,N);
 191:          /// if SIDE = 'R', LWORK .GE. max(1,M).
 192:          /// For optimum performance LWORK .GE. N*NB if SIDE = 'L', and
 193:          /// LWORK .GE. M*NB if SIDE = 'R', where NB is the optimal
 194:          /// blocksize.
 195:          /// 
 196:          /// If LWORK = -1, then a workspace query is assumed; the routine
 197:          /// only calculates the optimal size of the WORK array, returns
 198:          /// this value as the first entry of the WORK array, and no error
 199:          /// message related to LWORK is issued by XERBLA.
 200:          ///</param>
 201:          /// <param name="INFO">
 202:          /// (output) INTEGER
 203:          /// = 0:  successful exit
 204:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 205:          ///</param>
 206:          public void Run(string SIDE, string TRANS, int M, int N, int K, int L
 207:                           , double[] A, int offset_a, int LDA, double[] TAU, int offset_tau, ref double[] C, int offset_c, int LDC, ref double[] WORK, int offset_work
 208:                           , int LWORK, ref int INFO)
 209:          {
 210:   
 211:              #region Array Index Correction
 212:              
 213:               int o_a = -1 - LDA + offset_a;  int o_tau = -1 + offset_tau;  int o_c = -1 - LDC + offset_c; 
 214:               int o_work = -1 + offset_work;
 215:   
 216:              #endregion
 217:   
 218:   
 219:              #region Strings
 220:              
 221:              SIDE = SIDE.Substring(0, 1);  TRANS = TRANS.Substring(0, 1);  
 222:   
 223:              #endregion
 224:   
 225:   
 226:              #region Prolog
 227:              
 228:              // *
 229:              // *  -- LAPACK routine (version 3.1.1) --
 230:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 231:              // *     January 2007
 232:              // *
 233:              // *     .. Scalar Arguments ..
 234:              // *     ..
 235:              // *     .. Array Arguments ..
 236:              // *     ..
 237:              // *
 238:              // *  Purpose
 239:              // *  =======
 240:              // *
 241:              // *  DORMRZ overwrites the general real M-by-N matrix C with
 242:              // *
 243:              // *                  SIDE = 'L'     SIDE = 'R'
 244:              // *  TRANS = 'N':      Q * C          C * Q
 245:              // *  TRANS = 'T':      Q**T * C       C * Q**T
 246:              // *
 247:              // *  where Q is a real orthogonal matrix defined as the product of k
 248:              // *  elementary reflectors
 249:              // *
 250:              // *        Q = H(1) H(2) . . . H(k)
 251:              // *
 252:              // *  as returned by DTZRZF. Q is of order M if SIDE = 'L' and of order N
 253:              // *  if SIDE = 'R'.
 254:              // *
 255:              // *  Arguments
 256:              // *  =========
 257:              // *
 258:              // *  SIDE    (input) CHARACTER*1
 259:              // *          = 'L': apply Q or Q**T from the Left;
 260:              // *          = 'R': apply Q or Q**T from the Right.
 261:              // *
 262:              // *  TRANS   (input) CHARACTER*1
 263:              // *          = 'N':  No transpose, apply Q;
 264:              // *          = 'T':  Transpose, apply Q**T.
 265:              // *
 266:              // *  M       (input) INTEGER
 267:              // *          The number of rows of the matrix C. M >= 0.
 268:              // *
 269:              // *  N       (input) INTEGER
 270:              // *          The number of columns of the matrix C. N >= 0.
 271:              // *
 272:              // *  K       (input) INTEGER
 273:              // *          The number of elementary reflectors whose product defines
 274:              // *          the matrix Q.
 275:              // *          If SIDE = 'L', M >= K >= 0;
 276:              // *          if SIDE = 'R', N >= K >= 0.
 277:              // *
 278:              // *  L       (input) INTEGER
 279:              // *          The number of columns of the matrix A containing
 280:              // *          the meaningful part of the Householder reflectors.
 281:              // *          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
 282:              // *
 283:              // *  A       (input) DOUBLE PRECISION array, dimension
 284:              // *                               (LDA,M) if SIDE = 'L',
 285:              // *                               (LDA,N) if SIDE = 'R'
 286:              // *          The i-th row must contain the vector which defines the
 287:              // *          elementary reflector H(i), for i = 1,2,...,k, as returned by
 288:              // *          DTZRZF in the last k rows of its array argument A.
 289:              // *          A is modified by the routine but restored on exit.
 290:              // *
 291:              // *  LDA     (input) INTEGER
 292:              // *          The leading dimension of the array A. LDA >= max(1,K).
 293:              // *
 294:              // *  TAU     (input) DOUBLE PRECISION array, dimension (K)
 295:              // *          TAU(i) must contain the scalar factor of the elementary
 296:              // *          reflector H(i), as returned by DTZRZF.
 297:              // *
 298:              // *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 299:              // *          On entry, the M-by-N matrix C.
 300:              // *          On exit, C is overwritten by Q*C or Q**H*C or C*Q**H or C*Q.
 301:              // *
 302:              // *  LDC     (input) INTEGER
 303:              // *          The leading dimension of the array C. LDC >= max(1,M).
 304:              // *
 305:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 306:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 307:              // *
 308:              // *  LWORK   (input) INTEGER
 309:              // *          The dimension of the array WORK.
 310:              // *          If SIDE = 'L', LWORK >= max(1,N);
 311:              // *          if SIDE = 'R', LWORK >= max(1,M).
 312:              // *          For optimum performance LWORK >= N*NB if SIDE = 'L', and
 313:              // *          LWORK >= M*NB if SIDE = 'R', where NB is the optimal
 314:              // *          blocksize.
 315:              // *
 316:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 317:              // *          only calculates the optimal size of the WORK array, returns
 318:              // *          this value as the first entry of the WORK array, and no error
 319:              // *          message related to LWORK is issued by XERBLA.
 320:              // *
 321:              // *  INFO    (output) INTEGER
 322:              // *          = 0:  successful exit
 323:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 324:              // *
 325:              // *  Further Details
 326:              // *  ===============
 327:              // *
 328:              // *  Based on contributions by
 329:              // *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 330:              // *
 331:              // *  =====================================================================
 332:              // *
 333:              // *     .. Parameters ..
 334:              // *     ..
 335:              // *     .. Local Scalars ..
 336:              // *     ..
 337:              // *     .. Local Arrays ..
 338:              // *     ..
 339:              // *     .. External Functions ..
 340:              // *     ..
 341:              // *     .. External Subroutines ..
 342:              // *     ..
 343:              // *     .. Intrinsic Functions ..
 344:              //      INTRINSIC          MAX, MIN;
 345:              // *     ..
 346:              // *     .. Executable Statements ..
 347:              // *
 348:              // *     Test the input arguments
 349:              // *
 350:   
 351:              #endregion
 352:   
 353:   
 354:              #region Body
 355:              
 356:              INFO = 0;
 357:              LEFT = this._lsame.Run(SIDE, "L");
 358:              NOTRAN = this._lsame.Run(TRANS, "N");
 359:              LQUERY = (LWORK ==  - 1);
 360:              // *
 361:              // *     NQ is the order of Q and NW is the minimum dimension of WORK
 362:              // *
 363:              if (LEFT)
 364:              {
 365:                  NQ = M;
 366:                  NW = Math.Max(1, N);
 367:              }
 368:              else
 369:              {
 370:                  NQ = N;
 371:                  NW = Math.Max(1, M);
 372:              }
 373:              if (!LEFT && !this._lsame.Run(SIDE, "R"))
 374:              {
 375:                  INFO =  - 1;
 376:              }
 377:              else
 378:              {
 379:                  if (!NOTRAN && !this._lsame.Run(TRANS, "T"))
 380:                  {
 381:                      INFO =  - 2;
 382:                  }
 383:                  else
 384:                  {
 385:                      if (M < 0)
 386:                      {
 387:                          INFO =  - 3;
 388:                      }
 389:                      else
 390:                      {
 391:                          if (N < 0)
 392:                          {
 393:                              INFO =  - 4;
 394:                          }
 395:                          else
 396:                          {
 397:                              if (K < 0 || K > NQ)
 398:                              {
 399:                                  INFO =  - 5;
 400:                              }
 401:                              else
 402:                              {
 403:                                  if (L < 0 || (LEFT && (L > M)) || (!LEFT && (L > N)))
 404:                                  {
 405:                                      INFO =  - 6;
 406:                                  }
 407:                                  else
 408:                                  {
 409:                                      if (LDA < Math.Max(1, K))
 410:                                      {
 411:                                          INFO =  - 8;
 412:                                      }
 413:                                      else
 414:                                      {
 415:                                          if (LDC < Math.Max(1, M))
 416:                                          {
 417:                                              INFO =  - 11;
 418:                                          }
 419:                                      }
 420:                                  }
 421:                              }
 422:                          }
 423:                      }
 424:                  }
 425:              }
 426:              // *
 427:              if (INFO == 0)
 428:              {
 429:                  if (M == 0 || N == 0)
 430:                  {
 431:                      LWKOPT = 1;
 432:                  }
 433:                  else
 434:                  {
 435:                      // *
 436:                      // *           Determine the block size.  NB may be at most NBMAX, where
 437:                      // *           NBMAX is used to define the local array T.
 438:                      // *
 439:                      NB = Math.Min(NBMAX, this._ilaenv.Run(1, "DORMRQ", SIDE + TRANS, M, N, K,  - 1));
 440:                      LWKOPT = NW * NB;
 441:                  }
 442:                  WORK[1 + o_work] = LWKOPT;
 443:                  // *
 444:                  if (LWORK < Math.Max(1, NW) && !LQUERY)
 445:                  {
 446:                      INFO =  - 13;
 447:                  }
 448:              }
 449:              // *
 450:              if (INFO != 0)
 451:              {
 452:                  this._xerbla.Run("DORMRZ",  - INFO);
 453:                  return;
 454:              }
 455:              else
 456:              {
 457:                  if (LQUERY)
 458:                  {
 459:                      return;
 460:                  }
 461:              }
 462:              // *
 463:              // *     Quick return if possible
 464:              // *
 465:              if (M == 0 || N == 0)
 466:              {
 467:                  WORK[1 + o_work] = 1;
 468:                  return;
 469:              }
 470:              // *
 471:              NBMIN = 2;
 472:              LDWORK = NW;
 473:              if (NB > 1 && NB < K)
 474:              {
 475:                  IWS = NW * NB;
 476:                  if (LWORK < IWS)
 477:                  {
 478:                      NB = LWORK / LDWORK;
 479:                      NBMIN = Math.Max(2, this._ilaenv.Run(2, "DORMRQ", SIDE + TRANS, M, N, K,  - 1));
 480:                  }
 481:              }
 482:              else
 483:              {
 484:                  IWS = NW;
 485:              }
 486:              // *
 487:              if (NB < NBMIN || NB >= K)
 488:              {
 489:                  // *
 490:                  // *        Use unblocked code
 491:                  // *
 492:                  this._dormr3.Run(SIDE, TRANS, M, N, K, L
 493:                                   , A, offset_a, LDA, TAU, offset_tau, ref C, offset_c, LDC, ref WORK, offset_work
 494:                                   , ref IINFO);
 495:              }
 496:              else
 497:              {
 498:                  // *
 499:                  // *        Use blocked code
 500:                  // *
 501:                  if ((LEFT && !NOTRAN) || (!LEFT && NOTRAN))
 502:                  {
 503:                      I1 = 1;
 504:                      I2 = K;
 505:                      I3 = NB;
 506:                  }
 507:                  else
 508:                  {
 509:                      I1 = ((K - 1) / NB) * NB + 1;
 510:                      I2 = 1;
 511:                      I3 =  - NB;
 512:                  }
 513:                  // *
 514:                  if (LEFT)
 515:                  {
 516:                      NI = N;
 517:                      JC = 1;
 518:                      JA = M - L + 1;
 519:                  }
 520:                  else
 521:                  {
 522:                      MI = M;
 523:                      IC = 1;
 524:                      JA = N - L + 1;
 525:                  }
 526:                  // *
 527:                  if (NOTRAN)
 528:                  {
 529:                      FortranLib.Copy(ref TRANST , "T");
 530:                  }
 531:                  else
 532:                  {
 533:                      FortranLib.Copy(ref TRANST , "N");
 534:                  }
 535:                  // *
 536:                  for (I = I1; (I3 >= 0) ? (I <= I2) : (I >= I2); I += I3)
 537:                  {
 538:                      IB = Math.Min(NB, K - I + 1);
 539:                      // *
 540:                      // *           Form the triangular factor of the block reflector
 541:                      // *           H = H(i+ib-1) . . . H(i+1) H(i)
 542:                      // *
 543:                      this._dlarzt.Run("Backward", "Rowwise", L, IB, A, I+JA * LDA + o_a, LDA
 544:                                       , TAU, I + o_tau, ref T, offset_t, LDT);
 545:                      // *
 546:                      if (LEFT)
 547:                      {
 548:                          // *
 549:                          // *              H or H' is applied to C(i:m,1:n)
 550:                          // *
 551:                          MI = M - I + 1;
 552:                          IC = I;
 553:                      }
 554:                      else
 555:                      {
 556:                          // *
 557:                          // *              H or H' is applied to C(1:m,i:n)
 558:                          // *
 559:                          NI = N - I + 1;
 560:                          JC = I;
 561:                      }
 562:                      // *
 563:                      // *           Apply H or H'
 564:                      // *
 565:                      this._dlarzb.Run(SIDE, TRANST, "Backward", "Rowwise", MI, NI
 566:                                       , IB, L, A, I+JA * LDA + o_a, LDA, T, offset_t, LDT
 567:                                       , ref C, IC+JC * LDC + o_c, LDC, ref WORK, offset_work, LDWORK);
 568:                  }
 569:                  // *
 570:              }
 571:              // *
 572:              WORK[1 + o_work] = LWKOPT;
 573:              // *
 574:              return;
 575:              // *
 576:              // *     End of DORMRZ
 577:              // *
 578:   
 579:              #endregion
 580:   
 581:          }
 582:      }
 583:  }