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 auxiliary routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DLARFB applies a real block reflector H or its transpose H' to a
  27:      /// real m by n matrix C, from either the left or the right.
  28:      /// 
  29:      ///</summary>
  30:      public class DLARFB
  31:      {
  32:      
  33:   
  34:          #region Dependencies
  35:          
  36:          LSAME _lsame; DCOPY _dcopy; DGEMM _dgemm; DTRMM _dtrmm; 
  37:   
  38:          #endregion
  39:   
  40:   
  41:          #region Fields
  42:          
  43:          const double ONE = 1.0E+0; string TRANST = new string(' ', 1); int I = 0; int J = 0; 
  44:   
  45:          #endregion
  46:   
  47:          public DLARFB(LSAME lsame, DCOPY dcopy, DGEMM dgemm, DTRMM dtrmm)
  48:          {
  49:      
  50:   
  51:              #region Set Dependencies
  52:              
  53:              this._lsame = lsame; this._dcopy = dcopy; this._dgemm = dgemm; this._dtrmm = dtrmm; 
  54:   
  55:              #endregion
  56:   
  57:          }
  58:      
  59:          public DLARFB()
  60:          {
  61:      
  62:   
  63:              #region Dependencies (Initialization)
  64:              
  65:              LSAME lsame = new LSAME();
  66:              DCOPY dcopy = new DCOPY();
  67:              XERBLA xerbla = new XERBLA();
  68:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  69:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  70:   
  71:              #endregion
  72:   
  73:   
  74:              #region Set Dependencies
  75:              
  76:              this._lsame = lsame; this._dcopy = dcopy; this._dgemm = dgemm; this._dtrmm = dtrmm; 
  77:   
  78:              #endregion
  79:   
  80:          }
  81:          /// <summary>
  82:          /// Purpose
  83:          /// =======
  84:          /// 
  85:          /// DLARFB applies a real block reflector H or its transpose H' to a
  86:          /// real m by n matrix C, from either the left or the right.
  87:          /// 
  88:          ///</summary>
  89:          /// <param name="SIDE">
  90:          /// (input) CHARACTER*1
  91:          /// = 'L': apply H or H' from the Left
  92:          /// = 'R': apply H or H' from the Right
  93:          ///</param>
  94:          /// <param name="TRANS">
  95:          /// (input) CHARACTER*1
  96:          /// = 'N': apply H (No transpose)
  97:          /// = 'T': apply H' (Transpose)
  98:          ///</param>
  99:          /// <param name="DIRECT">
 100:          /// (input) CHARACTER*1
 101:          /// Indicates how H is formed from a product of elementary
 102:          /// reflectors
 103:          /// = 'F': H = H(1) H(2) . . . H(k) (Forward)
 104:          /// = 'B': H = H(k) . . . H(2) H(1) (Backward)
 105:          ///</param>
 106:          /// <param name="STOREV">
 107:          /// (input) CHARACTER*1
 108:          /// Indicates how the vectors which define the elementary
 109:          /// reflectors are stored:
 110:          /// = 'C': Columnwise
 111:          /// = 'R': Rowwise
 112:          ///</param>
 113:          /// <param name="M">
 114:          /// (input) INTEGER
 115:          /// The number of rows of the matrix C.
 116:          ///</param>
 117:          /// <param name="N">
 118:          /// (input) INTEGER
 119:          /// The number of columns of the matrix C.
 120:          ///</param>
 121:          /// <param name="K">
 122:          /// (input) INTEGER
 123:          /// The order of the matrix T (= the number of elementary
 124:          /// reflectors whose product defines the block reflector).
 125:          ///</param>
 126:          /// <param name="V">
 127:          /// (input) DOUBLE PRECISION array, dimension
 128:          /// (LDV,K) if STOREV = 'C'
 129:          /// (LDV,M) if STOREV = 'R' and SIDE = 'L'
 130:          /// (LDV,N) if STOREV = 'R' and SIDE = 'R'
 131:          /// The matrix V. See further details.
 132:          ///</param>
 133:          /// <param name="LDV">
 134:          /// (input) INTEGER
 135:          /// The leading dimension of the array V.
 136:          /// If STOREV = 'C' and SIDE = 'L', LDV .GE. max(1,M);
 137:          /// if STOREV = 'C' and SIDE = 'R', LDV .GE. max(1,N);
 138:          /// if STOREV = 'R', LDV .GE. K.
 139:          ///</param>
 140:          /// <param name="T">
 141:          /// (input) DOUBLE PRECISION array, dimension (LDT,K)
 142:          /// The triangular k by k matrix T in the representation of the
 143:          /// block reflector.
 144:          ///</param>
 145:          /// <param name="LDT">
 146:          /// (input) INTEGER
 147:          /// The leading dimension of the array T. LDT .GE. K.
 148:          ///</param>
 149:          /// <param name="C">
 150:          /// (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 151:          /// On entry, the m by n matrix C.
 152:          /// On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
 153:          ///</param>
 154:          /// <param name="LDC">
 155:          /// (input) INTEGER
 156:          /// The leading dimension of the array C. LDA .GE. max(1,M).
 157:          ///</param>
 158:          /// <param name="WORK">
 159:          /// (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
 160:          ///</param>
 161:          /// <param name="LDWORK">
 162:          /// (input) INTEGER
 163:          /// The leading dimension of the array WORK.
 164:          /// If SIDE = 'L', LDWORK .GE. max(1,N);
 165:          /// if SIDE = 'R', LDWORK .GE. max(1,M).
 166:          ///</param>
 167:          public void Run(string SIDE, string TRANS, string DIRECT, string STOREV, int M, int N
 168:                           , int K, double[] V, int offset_v, int LDV, double[] T, int offset_t, int LDT, ref double[] C, int offset_c
 169:                           , int LDC, ref double[] WORK, int offset_work, int LDWORK)
 170:          {
 171:   
 172:              #region Array Index Correction
 173:              
 174:               int o_v = -1 - LDV + offset_v;  int o_t = -1 - LDT + offset_t;  int o_c = -1 - LDC + offset_c; 
 175:               int o_work = -1 - LDWORK + offset_work;
 176:   
 177:              #endregion
 178:   
 179:   
 180:              #region Strings
 181:              
 182:              SIDE = SIDE.Substring(0, 1);  TRANS = TRANS.Substring(0, 1);  DIRECT = DIRECT.Substring(0, 1);  
 183:              STOREV = STOREV.Substring(0, 1); 
 184:   
 185:              #endregion
 186:   
 187:   
 188:              #region Prolog
 189:              
 190:              // *
 191:              // *  -- LAPACK auxiliary routine (version 3.1) --
 192:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 193:              // *     November 2006
 194:              // *
 195:              // *     .. Scalar Arguments ..
 196:              // *     ..
 197:              // *     .. Array Arguments ..
 198:              // *     ..
 199:              // *
 200:              // *  Purpose
 201:              // *  =======
 202:              // *
 203:              // *  DLARFB applies a real block reflector H or its transpose H' to a
 204:              // *  real m by n matrix C, from either the left or the right.
 205:              // *
 206:              // *  Arguments
 207:              // *  =========
 208:              // *
 209:              // *  SIDE    (input) CHARACTER*1
 210:              // *          = 'L': apply H or H' from the Left
 211:              // *          = 'R': apply H or H' from the Right
 212:              // *
 213:              // *  TRANS   (input) CHARACTER*1
 214:              // *          = 'N': apply H (No transpose)
 215:              // *          = 'T': apply H' (Transpose)
 216:              // *
 217:              // *  DIRECT  (input) CHARACTER*1
 218:              // *          Indicates how H is formed from a product of elementary
 219:              // *          reflectors
 220:              // *          = 'F': H = H(1) H(2) . . . H(k) (Forward)
 221:              // *          = 'B': H = H(k) . . . H(2) H(1) (Backward)
 222:              // *
 223:              // *  STOREV  (input) CHARACTER*1
 224:              // *          Indicates how the vectors which define the elementary
 225:              // *          reflectors are stored:
 226:              // *          = 'C': Columnwise
 227:              // *          = 'R': Rowwise
 228:              // *
 229:              // *  M       (input) INTEGER
 230:              // *          The number of rows of the matrix C.
 231:              // *
 232:              // *  N       (input) INTEGER
 233:              // *          The number of columns of the matrix C.
 234:              // *
 235:              // *  K       (input) INTEGER
 236:              // *          The order of the matrix T (= the number of elementary
 237:              // *          reflectors whose product defines the block reflector).
 238:              // *
 239:              // *  V       (input) DOUBLE PRECISION array, dimension
 240:              // *                                (LDV,K) if STOREV = 'C'
 241:              // *                                (LDV,M) if STOREV = 'R' and SIDE = 'L'
 242:              // *                                (LDV,N) if STOREV = 'R' and SIDE = 'R'
 243:              // *          The matrix V. See further details.
 244:              // *
 245:              // *  LDV     (input) INTEGER
 246:              // *          The leading dimension of the array V.
 247:              // *          If STOREV = 'C' and SIDE = 'L', LDV >= max(1,M);
 248:              // *          if STOREV = 'C' and SIDE = 'R', LDV >= max(1,N);
 249:              // *          if STOREV = 'R', LDV >= K.
 250:              // *
 251:              // *  T       (input) DOUBLE PRECISION array, dimension (LDT,K)
 252:              // *          The triangular k by k matrix T in the representation of the
 253:              // *          block reflector.
 254:              // *
 255:              // *  LDT     (input) INTEGER
 256:              // *          The leading dimension of the array T. LDT >= K.
 257:              // *
 258:              // *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
 259:              // *          On entry, the m by n matrix C.
 260:              // *          On exit, C is overwritten by H*C or H'*C or C*H or C*H'.
 261:              // *
 262:              // *  LDC     (input) INTEGER
 263:              // *          The leading dimension of the array C. LDA >= max(1,M).
 264:              // *
 265:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (LDWORK,K)
 266:              // *
 267:              // *  LDWORK  (input) INTEGER
 268:              // *          The leading dimension of the array WORK.
 269:              // *          If SIDE = 'L', LDWORK >= max(1,N);
 270:              // *          if SIDE = 'R', LDWORK >= max(1,M).
 271:              // *
 272:              // *  =====================================================================
 273:              // *
 274:              // *     .. Parameters ..
 275:              // *     ..
 276:              // *     .. Local Scalars ..
 277:              // *     ..
 278:              // *     .. External Functions ..
 279:              // *     ..
 280:              // *     .. External Subroutines ..
 281:              // *     ..
 282:              // *     .. Executable Statements ..
 283:              // *
 284:              // *     Quick return if possible
 285:              // *
 286:   
 287:              #endregion
 288:   
 289:   
 290:              #region Body
 291:              
 292:              if (M <= 0 || N <= 0) return;
 293:              // *
 294:              if (this._lsame.Run(TRANS, "N"))
 295:              {
 296:                  FortranLib.Copy(ref TRANST , "T");
 297:              }
 298:              else
 299:              {
 300:                  FortranLib.Copy(ref TRANST , "N");
 301:              }
 302:              // *
 303:              if (this._lsame.Run(STOREV, "C"))
 304:              {
 305:                  // *
 306:                  if (this._lsame.Run(DIRECT, "F"))
 307:                  {
 308:                      // *
 309:                      // *           Let  V =  ( V1 )    (first K rows)
 310:                      // *                     ( V2 )
 311:                      // *           where  V1  is unit lower triangular.
 312:                      // *
 313:                      if (this._lsame.Run(SIDE, "L"))
 314:                      {
 315:                          // *
 316:                          // *              Form  H * C  or  H' * C  where  C = ( C1 )
 317:                          // *                                                  ( C2 )
 318:                          // *
 319:                          // *              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)
 320:                          // *
 321:                          // *              W := C1'
 322:                          // *
 323:                          for (J = 1; J <= K; J++)
 324:                          {
 325:                              this._dcopy.Run(N, C, J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
 326:                          }
 327:                          // *
 328:                          // *              W := W * V1
 329:                          // *
 330:                          this._dtrmm.Run("Right", "Lower", "No transpose", "Unit", N, K
 331:                                          , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
 332:                          if (M > K)
 333:                          {
 334:                              // *
 335:                              // *                 W := W + C2'*V2
 336:                              // *
 337:                              this._dgemm.Run("Transpose", "No transpose", N, K, M - K, ONE
 338:                                              , C, K + 1+1 * LDC + o_c, LDC, V, K + 1+1 * LDV + o_v, LDV, ONE, ref WORK, offset_work
 339:                                              , LDWORK);
 340:                          }
 341:                          // *
 342:                          // *              W := W * T'  or  W * T
 343:                          // *
 344:                          this._dtrmm.Run("Right", "Upper", TRANST, "Non-unit", N, K
 345:                                          , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 346:                          // *
 347:                          // *              C := C - V * W'
 348:                          // *
 349:                          if (M > K)
 350:                          {
 351:                              // *
 352:                              // *                 C2 := C2 - V2 * W'
 353:                              // *
 354:                              this._dgemm.Run("No transpose", "Transpose", M - K, N, K,  - ONE
 355:                                              , V, K + 1+1 * LDV + o_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, K + 1+1 * LDC + o_c
 356:                                              , LDC);
 357:                          }
 358:                          // *
 359:                          // *              W := W * V1'
 360:                          // *
 361:                          this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", N, K
 362:                                          , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
 363:                          // *
 364:                          // *              C1 := C1 - W'
 365:                          // *
 366:                          for (J = 1; J <= K; J++)
 367:                          {
 368:                              for (I = 1; I <= N; I++)
 369:                              {
 370:                                  C[J+I * LDC + o_c] = C[J+I * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 371:                              }
 372:                          }
 373:                          // *
 374:                      }
 375:                      else
 376:                      {
 377:                          if (this._lsame.Run(SIDE, "R"))
 378:                          {
 379:                              // *
 380:                              // *              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 381:                              // *
 382:                              // *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
 383:                              // *
 384:                              // *              W := C1
 385:                              // *
 386:                              for (J = 1; J <= K; J++)
 387:                              {
 388:                                  this._dcopy.Run(M, C, 1+J * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
 389:                              }
 390:                              // *
 391:                              // *              W := W * V1
 392:                              // *
 393:                              this._dtrmm.Run("Right", "Lower", "No transpose", "Unit", M, K
 394:                                              , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
 395:                              if (N > K)
 396:                              {
 397:                                  // *
 398:                                  // *                 W := W + C2 * V2
 399:                                  // *
 400:                                  this._dgemm.Run("No transpose", "No transpose", M, K, N - K, ONE
 401:                                                  , C, 1+(K + 1) * LDC + o_c, LDC, V, K + 1+1 * LDV + o_v, LDV, ONE, ref WORK, offset_work
 402:                                                  , LDWORK);
 403:                              }
 404:                              // *
 405:                              // *              W := W * T  or  W * T'
 406:                              // *
 407:                              this._dtrmm.Run("Right", "Upper", TRANS, "Non-unit", M, K
 408:                                              , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 409:                              // *
 410:                              // *              C := C - W * V'
 411:                              // *
 412:                              if (N > K)
 413:                              {
 414:                                  // *
 415:                                  // *                 C2 := C2 - W * V2'
 416:                                  // *
 417:                                  this._dgemm.Run("No transpose", "Transpose", M, N - K, K,  - ONE
 418:                                                  , WORK, offset_work, LDWORK, V, K + 1+1 * LDV + o_v, LDV, ONE, ref C, 1+(K + 1) * LDC + o_c
 419:                                                  , LDC);
 420:                              }
 421:                              // *
 422:                              // *              W := W * V1'
 423:                              // *
 424:                              this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", M, K
 425:                                              , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
 426:                              // *
 427:                              // *              C1 := C1 - W
 428:                              // *
 429:                              for (J = 1; J <= K; J++)
 430:                              {
 431:                                  for (I = 1; I <= M; I++)
 432:                                  {
 433:                                      C[I+J * LDC + o_c] = C[I+J * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 434:                                  }
 435:                              }
 436:                          }
 437:                      }
 438:                      // *
 439:                  }
 440:                  else
 441:                  {
 442:                      // *
 443:                      // *           Let  V =  ( V1 )
 444:                      // *                     ( V2 )    (last K rows)
 445:                      // *           where  V2  is unit upper triangular.
 446:                      // *
 447:                      if (this._lsame.Run(SIDE, "L"))
 448:                      {
 449:                          // *
 450:                          // *              Form  H * C  or  H' * C  where  C = ( C1 )
 451:                          // *                                                  ( C2 )
 452:                          // *
 453:                          // *              W := C' * V  =  (C1'*V1 + C2'*V2)  (stored in WORK)
 454:                          // *
 455:                          // *              W := C2'
 456:                          // *
 457:                          for (J = 1; J <= K; J++)
 458:                          {
 459:                              this._dcopy.Run(N, C, M - K + J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
 460:                          }
 461:                          // *
 462:                          // *              W := W * V2
 463:                          // *
 464:                          this._dtrmm.Run("Right", "Upper", "No transpose", "Unit", N, K
 465:                                          , ONE, V, M - K + 1+1 * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
 466:                          if (M > K)
 467:                          {
 468:                              // *
 469:                              // *                 W := W + C1'*V1
 470:                              // *
 471:                              this._dgemm.Run("Transpose", "No transpose", N, K, M - K, ONE
 472:                                              , C, offset_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
 473:                                              , LDWORK);
 474:                          }
 475:                          // *
 476:                          // *              W := W * T'  or  W * T
 477:                          // *
 478:                          this._dtrmm.Run("Right", "Lower", TRANST, "Non-unit", N, K
 479:                                          , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 480:                          // *
 481:                          // *              C := C - V * W'
 482:                          // *
 483:                          if (M > K)
 484:                          {
 485:                              // *
 486:                              // *                 C1 := C1 - V1 * W'
 487:                              // *
 488:                              this._dgemm.Run("No transpose", "Transpose", M - K, N, K,  - ONE
 489:                                              , V, offset_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, offset_c
 490:                                              , LDC);
 491:                          }
 492:                          // *
 493:                          // *              W := W * V2'
 494:                          // *
 495:                          this._dtrmm.Run("Right", "Upper", "Transpose", "Unit", N, K
 496:                                          , ONE, V, M - K + 1+1 * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
 497:                          // *
 498:                          // *              C2 := C2 - W'
 499:                          // *
 500:                          for (J = 1; J <= K; J++)
 501:                          {
 502:                              for (I = 1; I <= N; I++)
 503:                              {
 504:                                  C[M - K + J+I * LDC + o_c] = C[M - K + J+I * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 505:                              }
 506:                          }
 507:                          // *
 508:                      }
 509:                      else
 510:                      {
 511:                          if (this._lsame.Run(SIDE, "R"))
 512:                          {
 513:                              // *
 514:                              // *              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 515:                              // *
 516:                              // *              W := C * V  =  (C1*V1 + C2*V2)  (stored in WORK)
 517:                              // *
 518:                              // *              W := C2
 519:                              // *
 520:                              for (J = 1; J <= K; J++)
 521:                              {
 522:                                  this._dcopy.Run(M, C, 1+(N - K + J) * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
 523:                              }
 524:                              // *
 525:                              // *              W := W * V2
 526:                              // *
 527:                              this._dtrmm.Run("Right", "Upper", "No transpose", "Unit", M, K
 528:                                              , ONE, V, N - K + 1+1 * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
 529:                              if (N > K)
 530:                              {
 531:                                  // *
 532:                                  // *                 W := W + C1 * V1
 533:                                  // *
 534:                                  this._dgemm.Run("No transpose", "No transpose", M, K, N - K, ONE
 535:                                                  , C, offset_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
 536:                                                  , LDWORK);
 537:                              }
 538:                              // *
 539:                              // *              W := W * T  or  W * T'
 540:                              // *
 541:                              this._dtrmm.Run("Right", "Lower", TRANS, "Non-unit", M, K
 542:                                              , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 543:                              // *
 544:                              // *              C := C - W * V'
 545:                              // *
 546:                              if (N > K)
 547:                              {
 548:                                  // *
 549:                                  // *                 C1 := C1 - W * V1'
 550:                                  // *
 551:                                  this._dgemm.Run("No transpose", "Transpose", M, N - K, K,  - ONE
 552:                                                  , WORK, offset_work, LDWORK, V, offset_v, LDV, ONE, ref C, offset_c
 553:                                                  , LDC);
 554:                              }
 555:                              // *
 556:                              // *              W := W * V2'
 557:                              // *
 558:                              this._dtrmm.Run("Right", "Upper", "Transpose", "Unit", M, K
 559:                                              , ONE, V, N - K + 1+1 * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
 560:                              // *
 561:                              // *              C2 := C2 - W
 562:                              // *
 563:                              for (J = 1; J <= K; J++)
 564:                              {
 565:                                  for (I = 1; I <= M; I++)
 566:                                  {
 567:                                      C[I+(N - K + J) * LDC + o_c] = C[I+(N - K + J) * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 568:                                  }
 569:                              }
 570:                          }
 571:                      }
 572:                  }
 573:                  // *
 574:              }
 575:              else
 576:              {
 577:                  if (this._lsame.Run(STOREV, "R"))
 578:                  {
 579:                      // *
 580:                      if (this._lsame.Run(DIRECT, "F"))
 581:                      {
 582:                          // *
 583:                          // *           Let  V =  ( V1  V2 )    (V1: first K columns)
 584:                          // *           where  V1  is unit upper triangular.
 585:                          // *
 586:                          if (this._lsame.Run(SIDE, "L"))
 587:                          {
 588:                              // *
 589:                              // *              Form  H * C  or  H' * C  where  C = ( C1 )
 590:                              // *                                                  ( C2 )
 591:                              // *
 592:                              // *              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)
 593:                              // *
 594:                              // *              W := C1'
 595:                              // *
 596:                              for (J = 1; J <= K; J++)
 597:                              {
 598:                                  this._dcopy.Run(N, C, J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
 599:                              }
 600:                              // *
 601:                              // *              W := W * V1'
 602:                              // *
 603:                              this._dtrmm.Run("Right", "Upper", "Transpose", "Unit", N, K
 604:                                              , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
 605:                              if (M > K)
 606:                              {
 607:                                  // *
 608:                                  // *                 W := W + C2'*V2'
 609:                                  // *
 610:                                  this._dgemm.Run("Transpose", "Transpose", N, K, M - K, ONE
 611:                                                  , C, K + 1+1 * LDC + o_c, LDC, V, 1+(K + 1) * LDV + o_v, LDV, ONE, ref WORK, offset_work
 612:                                                  , LDWORK);
 613:                              }
 614:                              // *
 615:                              // *              W := W * T'  or  W * T
 616:                              // *
 617:                              this._dtrmm.Run("Right", "Upper", TRANST, "Non-unit", N, K
 618:                                              , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 619:                              // *
 620:                              // *              C := C - V' * W'
 621:                              // *
 622:                              if (M > K)
 623:                              {
 624:                                  // *
 625:                                  // *                 C2 := C2 - V2' * W'
 626:                                  // *
 627:                                  this._dgemm.Run("Transpose", "Transpose", M - K, N, K,  - ONE
 628:                                                  , V, 1+(K + 1) * LDV + o_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, K + 1+1 * LDC + o_c
 629:                                                  , LDC);
 630:                              }
 631:                              // *
 632:                              // *              W := W * V1
 633:                              // *
 634:                              this._dtrmm.Run("Right", "Upper", "No transpose", "Unit", N, K
 635:                                              , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
 636:                              // *
 637:                              // *              C1 := C1 - W'
 638:                              // *
 639:                              for (J = 1; J <= K; J++)
 640:                              {
 641:                                  for (I = 1; I <= N; I++)
 642:                                  {
 643:                                      C[J+I * LDC + o_c] = C[J+I * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 644:                                  }
 645:                              }
 646:                              // *
 647:                          }
 648:                          else
 649:                          {
 650:                              if (this._lsame.Run(SIDE, "R"))
 651:                              {
 652:                                  // *
 653:                                  // *              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 654:                                  // *
 655:                                  // *              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)
 656:                                  // *
 657:                                  // *              W := C1
 658:                                  // *
 659:                                  for (J = 1; J <= K; J++)
 660:                                  {
 661:                                      this._dcopy.Run(M, C, 1+J * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
 662:                                  }
 663:                                  // *
 664:                                  // *              W := W * V1'
 665:                                  // *
 666:                                  this._dtrmm.Run("Right", "Upper", "Transpose", "Unit", M, K
 667:                                                  , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
 668:                                  if (N > K)
 669:                                  {
 670:                                      // *
 671:                                      // *                 W := W + C2 * V2'
 672:                                      // *
 673:                                      this._dgemm.Run("No transpose", "Transpose", M, K, N - K, ONE
 674:                                                      , C, 1+(K + 1) * LDC + o_c, LDC, V, 1+(K + 1) * LDV + o_v, LDV, ONE, ref WORK, offset_work
 675:                                                      , LDWORK);
 676:                                  }
 677:                                  // *
 678:                                  // *              W := W * T  or  W * T'
 679:                                  // *
 680:                                  this._dtrmm.Run("Right", "Upper", TRANS, "Non-unit", M, K
 681:                                                  , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 682:                                  // *
 683:                                  // *              C := C - W * V
 684:                                  // *
 685:                                  if (N > K)
 686:                                  {
 687:                                      // *
 688:                                      // *                 C2 := C2 - W * V2
 689:                                      // *
 690:                                      this._dgemm.Run("No transpose", "No transpose", M, N - K, K,  - ONE
 691:                                                      , WORK, offset_work, LDWORK, V, 1+(K + 1) * LDV + o_v, LDV, ONE, ref C, 1+(K + 1) * LDC + o_c
 692:                                                      , LDC);
 693:                                  }
 694:                                  // *
 695:                                  // *              W := W * V1
 696:                                  // *
 697:                                  this._dtrmm.Run("Right", "Upper", "No transpose", "Unit", M, K
 698:                                                  , ONE, V, offset_v, LDV, ref WORK, offset_work, LDWORK);
 699:                                  // *
 700:                                  // *              C1 := C1 - W
 701:                                  // *
 702:                                  for (J = 1; J <= K; J++)
 703:                                  {
 704:                                      for (I = 1; I <= M; I++)
 705:                                      {
 706:                                          C[I+J * LDC + o_c] = C[I+J * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 707:                                      }
 708:                                  }
 709:                                  // *
 710:                              }
 711:                          }
 712:                          // *
 713:                      }
 714:                      else
 715:                      {
 716:                          // *
 717:                          // *           Let  V =  ( V1  V2 )    (V2: last K columns)
 718:                          // *           where  V2  is unit lower triangular.
 719:                          // *
 720:                          if (this._lsame.Run(SIDE, "L"))
 721:                          {
 722:                              // *
 723:                              // *              Form  H * C  or  H' * C  where  C = ( C1 )
 724:                              // *                                                  ( C2 )
 725:                              // *
 726:                              // *              W := C' * V'  =  (C1'*V1' + C2'*V2') (stored in WORK)
 727:                              // *
 728:                              // *              W := C2'
 729:                              // *
 730:                              for (J = 1; J <= K; J++)
 731:                              {
 732:                                  this._dcopy.Run(N, C, M - K + J+1 * LDC + o_c, LDC, ref WORK, 1+J * LDWORK + o_work, 1);
 733:                              }
 734:                              // *
 735:                              // *              W := W * V2'
 736:                              // *
 737:                              this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", N, K
 738:                                              , ONE, V, 1+(M - K + 1) * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
 739:                              if (M > K)
 740:                              {
 741:                                  // *
 742:                                  // *                 W := W + C1'*V1'
 743:                                  // *
 744:                                  this._dgemm.Run("Transpose", "Transpose", N, K, M - K, ONE
 745:                                                  , C, offset_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
 746:                                                  , LDWORK);
 747:                              }
 748:                              // *
 749:                              // *              W := W * T'  or  W * T
 750:                              // *
 751:                              this._dtrmm.Run("Right", "Lower", TRANST, "Non-unit", N, K
 752:                                              , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 753:                              // *
 754:                              // *              C := C - V' * W'
 755:                              // *
 756:                              if (M > K)
 757:                              {
 758:                                  // *
 759:                                  // *                 C1 := C1 - V1' * W'
 760:                                  // *
 761:                                  this._dgemm.Run("Transpose", "Transpose", M - K, N, K,  - ONE
 762:                                                  , V, offset_v, LDV, WORK, offset_work, LDWORK, ONE, ref C, offset_c
 763:                                                  , LDC);
 764:                              }
 765:                              // *
 766:                              // *              W := W * V2
 767:                              // *
 768:                              this._dtrmm.Run("Right", "Lower", "No transpose", "Unit", N, K
 769:                                              , ONE, V, 1+(M - K + 1) * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
 770:                              // *
 771:                              // *              C2 := C2 - W'
 772:                              // *
 773:                              for (J = 1; J <= K; J++)
 774:                              {
 775:                                  for (I = 1; I <= N; I++)
 776:                                  {
 777:                                      C[M - K + J+I * LDC + o_c] = C[M - K + J+I * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 778:                                  }
 779:                              }
 780:                              // *
 781:                          }
 782:                          else
 783:                          {
 784:                              if (this._lsame.Run(SIDE, "R"))
 785:                              {
 786:                                  // *
 787:                                  // *              Form  C * H  or  C * H'  where  C = ( C1  C2 )
 788:                                  // *
 789:                                  // *              W := C * V'  =  (C1*V1' + C2*V2')  (stored in WORK)
 790:                                  // *
 791:                                  // *              W := C2
 792:                                  // *
 793:                                  for (J = 1; J <= K; J++)
 794:                                  {
 795:                                      this._dcopy.Run(M, C, 1+(N - K + J) * LDC + o_c, 1, ref WORK, 1+J * LDWORK + o_work, 1);
 796:                                  }
 797:                                  // *
 798:                                  // *              W := W * V2'
 799:                                  // *
 800:                                  this._dtrmm.Run("Right", "Lower", "Transpose", "Unit", M, K
 801:                                                  , ONE, V, 1+(N - K + 1) * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
 802:                                  if (N > K)
 803:                                  {
 804:                                      // *
 805:                                      // *                 W := W + C1 * V1'
 806:                                      // *
 807:                                      this._dgemm.Run("No transpose", "Transpose", M, K, N - K, ONE
 808:                                                      , C, offset_c, LDC, V, offset_v, LDV, ONE, ref WORK, offset_work
 809:                                                      , LDWORK);
 810:                                  }
 811:                                  // *
 812:                                  // *              W := W * T  or  W * T'
 813:                                  // *
 814:                                  this._dtrmm.Run("Right", "Lower", TRANS, "Non-unit", M, K
 815:                                                  , ONE, T, offset_t, LDT, ref WORK, offset_work, LDWORK);
 816:                                  // *
 817:                                  // *              C := C - W * V
 818:                                  // *
 819:                                  if (N > K)
 820:                                  {
 821:                                      // *
 822:                                      // *                 C1 := C1 - W * V1
 823:                                      // *
 824:                                      this._dgemm.Run("No transpose", "No transpose", M, N - K, K,  - ONE
 825:                                                      , WORK, offset_work, LDWORK, V, offset_v, LDV, ONE, ref C, offset_c
 826:                                                      , LDC);
 827:                                  }
 828:                                  // *
 829:                                  // *              W := W * V2
 830:                                  // *
 831:                                  this._dtrmm.Run("Right", "Lower", "No transpose", "Unit", M, K
 832:                                                  , ONE, V, 1+(N - K + 1) * LDV + o_v, LDV, ref WORK, offset_work, LDWORK);
 833:                                  // *
 834:                                  // *              C1 := C1 - W
 835:                                  // *
 836:                                  for (J = 1; J <= K; J++)
 837:                                  {
 838:                                      for (I = 1; I <= M; I++)
 839:                                      {
 840:                                          C[I+(N - K + J) * LDC + o_c] = C[I+(N - K + J) * LDC + o_c] - WORK[I+J * LDWORK + o_work];
 841:                                      }
 842:                                  }
 843:                                  // *
 844:                              }
 845:                          }
 846:                          // *
 847:                      }
 848:                  }
 849:              }
 850:              // *
 851:              return;
 852:              // *
 853:              // *     End of DLARFB
 854:              // *
 855:   
 856:              #endregion
 857:   
 858:          }
 859:      }
 860:  }