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:      /// DLARZB applies a real block reflector H or its transpose H**T to
  27:      /// a real distributed M-by-N  C from the left or the right.
  28:      /// 
  29:      /// Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
  30:      /// 
  31:      ///</summary>
  32:      public class DLARZB
  33:      {
  34:      
  35:   
  36:          #region Dependencies
  37:          
  38:          LSAME _lsame; DCOPY _dcopy; DGEMM _dgemm; DTRMM _dtrmm; XERBLA _xerbla; 
  39:   
  40:          #endregion
  41:   
  42:   
  43:          #region Fields
  44:          
  45:          const double ONE = 1.0E+0; string TRANST = new string(' ', 1); int I = 0; int INFO = 0; int J = 0; 
  46:   
  47:          #endregion
  48:   
  49:          public DLARZB(LSAME lsame, DCOPY dcopy, DGEMM dgemm, DTRMM dtrmm, XERBLA xerbla)
  50:          {
  51:      
  52:   
  53:              #region Set Dependencies
  54:              
  55:              this._lsame = lsame; this._dcopy = dcopy; this._dgemm = dgemm; this._dtrmm = dtrmm; this._xerbla = xerbla; 
  56:   
  57:              #endregion
  58:   
  59:          }
  60:      
  61:          public DLARZB()
  62:          {
  63:      
  64:   
  65:              #region Dependencies (Initialization)
  66:              
  67:              LSAME lsame = new LSAME();
  68:              DCOPY dcopy = new DCOPY();
  69:              XERBLA xerbla = new XERBLA();
  70:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  71:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  72:   
  73:              #endregion
  74:   
  75:   
  76:              #region Set Dependencies
  77:              
  78:              this._lsame = lsame; this._dcopy = dcopy; this._dgemm = dgemm; this._dtrmm = dtrmm; this._xerbla = xerbla; 
  79:   
  80:              #endregion
  81:   
  82:          }
  83:          /// <summary>
  84:          /// Purpose
  85:          /// =======
  86:          /// 
  87:          /// DLARZB applies a real block reflector H or its transpose H**T to
  88:          /// a real distributed M-by-N  C from the left or the right.
  89:          /// 
  90:          /// Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
  91:          /// 
  92:          ///</summary>
  93:          /// <param name="SIDE">
  94:          /// (input) CHARACTER*1
  95:          /// = 'L': apply H or H' from the Left
  96:          /// = 'R': apply H or H' from the Right
  97:          ///</param>
  98:          /// <param name="TRANS">
  99:          /// (input) CHARACTER*1
 100:          /// = 'N': apply H (No transpose)
 101:          /// = 'C': apply H' (Transpose)
 102:          ///</param>
 103:          /// <param name="DIRECT">
 104:          /// (input) CHARACTER*1
 105:          /// Indicates how H is formed from a product of elementary
 106:          /// reflectors
 107:          /// = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
 108:          /// = 'B': H = H(k) . . . H(2) H(1) (Backward)
 109:          ///</param>
 110:          /// <param name="STOREV">
 111:          /// (input) CHARACTER*1
 112:          /// Indicates how the vectors which define the elementary
 113:          /// reflectors are stored:
 114:          /// = 'C': Columnwise                        (not supported yet)
 115:          /// = 'R': Rowwise
 116:          ///</param>
 117:          /// <param name="M">
 118:          /// (input) INTEGER
 119:          /// The number of rows of the matrix C.
 120:          ///</param>
 121:          /// <param name="N">
 122:          /// (input) INTEGER
 123:          /// The number of columns of the matrix C.
 124:          ///</param>
 125:          /// <param name="K">
 126:          /// (input) INTEGER
 127:          /// The order of the matrix T (= the number of elementary
 128:          /// reflectors whose product defines the block reflector).
 129:          ///</param>
 130:          /// <param name="L">
 131:          /// (input) INTEGER
 132:          /// The number of columns of the matrix V containing the
 133:          /// meaningful part of the Householder reflectors.
 134:          /// If SIDE = 'L', M .GE. L .GE. 0, if SIDE = 'R', N .GE. L .GE. 0.
 135:          ///</param>
 136:          /// <param name="V">
 137:          /// (input) DOUBLE PRECISION array, dimension (LDV,NV).
 138:          /// If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
 139:          ///</param>
 140:          /// <param name="LDV">
 141:          /// (input) INTEGER
 142:          /// The leading dimension of the array V.
 143:          /// If STOREV = 'C', LDV .GE. L; if STOREV = 'R', LDV .GE. K.
 144:          ///</param>
 145:          /// <param name="T">
 146:          /// (input) DOUBLE PRECISION array, dimension (LDT,K)
 147:          /// The triangular K-by-K matrix T in the representation of the
 148:          /// block reflector.
 149:          ///</param>
 150:          /// <param name="LDT">
 151:          /// (input) INTEGER
 152:          /// The leading dimension of the array T. LDT .GE. K.
 153:          ///</param>
 154:          /// <param name="C">
 155:          /// (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 156:          /// On entry, the M-by-N matrix C.
 157:          /// On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
 158:          ///</param>
 159:          /// <param name="LDC">
 160:          /// (input) INTEGER
 161:          /// The leading dimension of the array C. LDC .GE. max(1,M).
 162:          ///</param>
 163:          /// <param name="WORK">
 164:          /// (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
 165:          ///</param>
 166:          /// <param name="LDWORK">
 167:          /// (input) INTEGER
 168:          /// The leading dimension of the array WORK.
 169:          /// If SIDE = 'L', LDWORK .GE. max(1,N);
 170:          /// if SIDE = 'R', LDWORK .GE. max(1,M).
 171:          ///</param>
 172:          public void Run(string SIDE, string TRANS, string DIRECT, string STOREV, int M, int N
 173:                           , int K, int L, double[] V, int offset_v, int LDV, double[] T, int offset_t, int LDT
 174:                           , ref double[] C, int offset_c, int LDC, ref double[] WORK, int offset_work, int LDWORK)
 175:          {
 176:   
 177:              #region Array Index Correction
 178:              
 179:               int o_v = -1 - LDV + offset_v;  int o_t = -1 - LDT + offset_t;  int o_c = -1 - LDC + offset_c; 
 180:               int o_work = -1 - LDWORK + offset_work;
 181:   
 182:              #endregion
 183:   
 184:   
 185:              #region Strings
 186:              
 187:              SIDE = SIDE.Substring(0, 1);  TRANS = TRANS.Substring(0, 1);  DIRECT = DIRECT.Substring(0, 1);  
 188:              STOREV = STOREV.Substring(0, 1); 
 189:   
 190:              #endregion
 191:   
 192:   
 193:              #region Prolog
 194:              
 195:              // *
 196:              // *  -- LAPACK routine (version 3.1) --
 197:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 198:              // *     November 2006
 199:              // *
 200:              // *     .. Scalar Arguments ..
 201:              // *     ..
 202:              // *     .. Array Arguments ..
 203:              // *     ..
 204:              // *
 205:              // *  Purpose
 206:              // *  =======
 207:              // *
 208:              // *  DLARZB applies a real block reflector H or its transpose H**T to
 209:              // *  a real distributed M-by-N  C from the left or the right.
 210:              // *
 211:              // *  Currently, only STOREV = 'R' and DIRECT = 'B' are supported.
 212:              // *
 213:              // *  Arguments
 214:              // *  =========
 215:              // *
 216:              // *  SIDE    (input) CHARACTER*1
 217:              // *          = 'L': apply H or H' from the Left
 218:              // *          = 'R': apply H or H' from the Right
 219:              // *
 220:              // *  TRANS   (input) CHARACTER*1
 221:              // *          = 'N': apply H (No transpose)
 222:              // *          = 'C': apply H' (Transpose)
 223:              // *
 224:              // *  DIRECT  (input) CHARACTER*1
 225:              // *          Indicates how H is formed from a product of elementary
 226:              // *          reflectors
 227:              // *          = 'F': H = H(1) H(2) . . . H(k) (Forward, not supported yet)
 228:              // *          = 'B': H = H(k) . . . H(2) H(1) (Backward)
 229:              // *
 230:              // *  STOREV  (input) CHARACTER*1
 231:              // *          Indicates how the vectors which define the elementary
 232:              // *          reflectors are stored:
 233:              // *          = 'C': Columnwise                        (not supported yet)
 234:              // *          = 'R': Rowwise
 235:              // *
 236:              // *  M       (input) INTEGER
 237:              // *          The number of rows of the matrix C.
 238:              // *
 239:              // *  N       (input) INTEGER
 240:              // *          The number of columns of the matrix C.
 241:              // *
 242:              // *  K       (input) INTEGER
 243:              // *          The order of the matrix T (= the number of elementary
 244:              // *          reflectors whose product defines the block reflector).
 245:              // *
 246:              // *  L       (input) INTEGER
 247:              // *          The number of columns of the matrix V containing the
 248:              // *          meaningful part of the Householder reflectors.
 249:              // *          If SIDE = 'L', M >= L >= 0, if SIDE = 'R', N >= L >= 0.
 250:              // *
 251:              // *  V       (input) DOUBLE PRECISION array, dimension (LDV,NV).
 252:              // *          If STOREV = 'C', NV = K; if STOREV = 'R', NV = L.
 253:              // *
 254:              // *  LDV     (input) INTEGER
 255:              // *          The leading dimension of the array V.
 256:              // *          If STOREV = 'C', LDV >= L; if STOREV = 'R', LDV >= K.
 257:              // *
 258:              // *  T       (input) DOUBLE PRECISION array, dimension (LDT,K)
 259:              // *          The triangular K-by-K matrix T in the representation of the
 260:              // *          block reflector.
 261:              // *
 262:              // *  LDT     (input) INTEGER
 263:              // *          The leading dimension of the array T. LDT >= K.
 264:              // *
 265:              // *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 266:              // *          On entry, the M-by-N matrix C.
 267:              // *          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
 268:              // *
 269:              // *  LDC     (input) INTEGER
 270:              // *          The leading dimension of the array C. LDC >= max(1,M).
 271:              // *
 272:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
 273:              // *
 274:              // *  LDWORK  (input) INTEGER
 275:              // *          The leading dimension of the array WORK.
 276:              // *          If SIDE = 'L', LDWORK >= max(1,N);
 277:              // *          if SIDE = 'R', LDWORK >= max(1,M).
 278:              // *
 279:              // *  Further Details
 280:              // *  ===============
 281:              // *
 282:              // *  Based on contributions by
 283:              // *    A. Petitet, Computer Science Dept., Univ. of Tenn., Knoxville, USA
 284:              // *
 285:              // *  =====================================================================
 286:              // *
 287:              // *     .. Parameters ..
 288:              // *     ..
 289:              // *     .. Local Scalars ..
 290:              // *     ..
 291:              // *     .. External Functions ..
 292:              // *     ..
 293:              // *     .. External Subroutines ..
 294:              // *     ..
 295:              // *     .. Executable Statements ..
 296:              // *
 297:              // *     Quick return if possible
 298:              // *
 299:   
 300:              #endregion
 301:   
 302:   
 303:              #region Body
 304:              
 305:              if (M <= 0 || N <= 0) return;
 306:              // *
 307:              // *     Check for currently supported options
 308:              // *
 309:              INFO = 0;
 310:              if (!this._lsame.Run(DIRECT, "B"))
 311:              {
 312:                  INFO =  - 3;
 313:              }
 314:              else
 315:              {
 316:                  if (!this._lsame.Run(STOREV, "R"))
 317:                  {
 318:                      INFO =  - 4;
 319:                  }
 320:              }
 321:              if (INFO != 0)
 322:              {
 323:                  this._xerbla.Run("DLARZB",  - INFO);
 324:                  return;
 325:              }
 326:              // *
 327:              if (this._lsame.Run(TRANS, "N"))
 328:              {
 329:                  FortranLib.Copy(ref TRANST , "T");
 330:              }
 331:              else
 332:              {
 333:                  FortranLib.Copy(ref TRANST , "N");
 334:              }
 335:              // *
 336:              if (this._lsame.Run(SIDE, "L"))
 337:              {
 338:                  // *
 339:                  // *        Form  H * C  or  H' * C
 340:                  // *
 341:                  // *        W( 1:n, 1:k ) = C( 1:k, 1:n )'
 342:                  // *
 343:                  for (J = 1; J <= K; J++)
 344:                  {
 345:                      this._dcopy.Run(N, C, J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
 346:                  }
 347:                  // *
 348:                  // *        W( 1:n, 1:k ) = W( 1:n, 1:k ) + ...
 349:                  // *                        C( m-l+1:m, 1:n )' * V( 1:k, 1:l )'
 350:                  // *
 351:                  if (L > 0)
 352:                  {
 353:                      this._dgemm.Run("Transpose", "Transpose", N, K, L, ONE
 354:                                      , C, M - L + 1+1 * LDC + o_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
 355:                                      , LDWORK);
 356:                  }
 357:                  // *
 358:                  // *        W( 1:n, 1:k ) = W( 1:n, 1:k ) * T'  or  W( 1:m, 1:k ) * T
 359:                  // *
 360:                  this._dtrmm.Run("Right", "Lower", TRANST, "Non-unit", N, K
 361:                                  , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 362:                  // *
 363:                  // *        C( 1:k, 1:n ) = C( 1:k, 1:n ) - W( 1:n, 1:k )'
 364:                  // *
 365:                  for (J = 1; J <= N; J++)
 366:                  {
 367:                      for (I = 1; I <= K; I++)
 368:                      {
 369:                          C[I+J * LDC + o_c] = C[I+J * LDC + o_c] - WORK[J+I * LDWORK + o_work];
 370:                      }
 371:                  }
 372:                  // *
 373:                  // *        C( m-l+1:m, 1:n ) = C( m-l+1:m, 1:n ) - ...
 374:                  // *                            V( 1:k, 1:l )' * W( 1:n, 1:k )'
 375:                  // *
 376:                  if (L > 0)
 377:                  {
 378:                      this._dgemm.Run("Transpose", "Transpose", L, N, K,  - ONE
 379:                                      , V, offset_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, M - L + 1+1 * LDC + o_c
 380:                                      , LDC);
 381:                  }
 382:                  // *
 383:              }
 384:              else
 385:              {
 386:                  if (this._lsame.Run(SIDE, "R"))
 387:                  {
 388:                      // *
 389:                      // *        Form  C * H  or  C * H'
 390:                      // *
 391:                      // *        W( 1:m, 1:k ) = C( 1:m, 1:k )
 392:                      // *
 393:                      for (J = 1; J <= K; J++)
 394:                      {
 395:                          this._dcopy.Run(M, C, 1+J * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
 396:                      }
 397:                      // *
 398:                      // *        W( 1:m, 1:k ) = W( 1:m, 1:k ) + ...
 399:                      // *                        C( 1:m, n-l+1:n ) * V( 1:k, 1:l )'
 400:                      // *
 401:                      if (L > 0)
 402:                      {
 403:                          this._dgemm.Run("No transpose", "Transpose", M, K, L, ONE
 404:                                          , C, 1+(N - L + 1) * LDC + o_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
 405:                                          , LDWORK);
 406:                      }
 407:                      // *
 408:                      // *        W( 1:m, 1:k ) = W( 1:m, 1:k ) * T  or  W( 1:m, 1:k ) * T'
 409:                      // *
 410:                      this._dtrmm.Run("Right", "Lower", TRANS, "Non-unit", M, K
 411:                                      , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 412:                      // *
 413:                      // *        C( 1:m, 1:k ) = C( 1:m, 1:k ) - W( 1:m, 1:k )
 414:                      // *
 415:                      for (J = 1; J <= K; J++)
 416:                      {
 417:                          for (I = 1; I <= M; I++)
 418:                          {
 419:                              C[I+J * LDC + o_c] = C[I+J * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 420:                          }
 421:                      }
 422:                      // *
 423:                      // *        C( 1:m, n-l+1:n ) = C( 1:m, n-l+1:n ) - ...
 424:                      // *                            W( 1:m, 1:k ) * V( 1:k, 1:l )
 425:                      // *
 426:                      if (L > 0)
 427:                      {
 428:                          this._dgemm.Run("No transpose", "No transpose", M, L, K,  - ONE
 429:                                          , WORK, offset_work, LDWORK, V, offset_v, LDV, ONE, ref C, 1+(N - L + 1) * LDC + o_c
 430:                                          , LDC);
 431:                      }
 432:                      // *
 433:                  }
 434:              }
 435:              // *
 436:              return;
 437:              // *
 438:              // *     End of DLARZB
 439:              // *
 440:   
 441:              #endregion
 442:   
 443:          }
 444:      }
 445:  }