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:      /// Purpose
  21:      /// =======
  22:      /// 
  23:      /// DTRSM  solves one of the matrix equations
  24:      /// 
  25:      /// op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
  26:      /// 
  27:      /// where alpha is a scalar, X and B are m by n matrices, A is a unit, or
  28:      /// non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
  29:      /// 
  30:      /// op( A ) = A   or   op( A ) = A'.
  31:      /// 
  32:      /// The matrix X is overwritten on B.
  33:      /// 
  34:      ///</summary>
  35:      public class DTRSM
  36:      {
  37:      
  38:   
  39:          #region Dependencies
  40:          
  41:          LSAME _lsame; XERBLA _xerbla; 
  42:   
  43:          #endregion
  44:   
  45:   
  46:          #region Fields
  47:          
  48:          double TEMP = 0; int I = 0; int INFO = 0; int J = 0; int K = 0; int NROWA = 0; bool LSIDE = false; bool NOUNIT = false; 
  49:          bool UPPER = false;const double ONE = 1.0E+0; const double ZERO = 0.0E+0; 
  50:   
  51:          #endregion
  52:   
  53:          public DTRSM(LSAME lsame, XERBLA xerbla)
  54:          {
  55:      
  56:   
  57:              #region Set Dependencies
  58:              
  59:              this._lsame = lsame; this._xerbla = xerbla; 
  60:   
  61:              #endregion
  62:   
  63:          }
  64:      
  65:          public DTRSM()
  66:          {
  67:      
  68:   
  69:              #region Dependencies (Initialization)
  70:              
  71:              LSAME lsame = new LSAME();
  72:              XERBLA xerbla = new XERBLA();
  73:   
  74:              #endregion
  75:   
  76:   
  77:              #region Set Dependencies
  78:              
  79:              this._lsame = lsame; this._xerbla = xerbla; 
  80:   
  81:              #endregion
  82:   
  83:          }
  84:          /// <summary>
  85:          /// Purpose
  86:          /// =======
  87:          /// 
  88:          /// DTRSM  solves one of the matrix equations
  89:          /// 
  90:          /// op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
  91:          /// 
  92:          /// where alpha is a scalar, X and B are m by n matrices, A is a unit, or
  93:          /// non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
  94:          /// 
  95:          /// op( A ) = A   or   op( A ) = A'.
  96:          /// 
  97:          /// The matrix X is overwritten on B.
  98:          /// 
  99:          ///</summary>
 100:          /// <param name="SIDE">
 101:          /// - CHARACTER*1.
 102:          /// On entry, SIDE specifies whether op( A ) appears on the left
 103:          /// or right of X as follows:
 104:          /// 
 105:          /// SIDE = 'L' or 'l'   op( A )*X = alpha*B.
 106:          /// 
 107:          /// SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
 108:          /// 
 109:          /// Unchanged on exit.
 110:          ///</param>
 111:          /// <param name="UPLO">
 112:          /// - CHARACTER*1.
 113:          /// On entry, UPLO specifies whether the matrix A is an upper or
 114:          /// lower triangular matrix as follows:
 115:          /// 
 116:          /// UPLO = 'U' or 'u'   A is an upper triangular matrix.
 117:          /// 
 118:          /// UPLO = 'L' or 'l'   A is a lower triangular matrix.
 119:          /// 
 120:          /// Unchanged on exit.
 121:          ///</param>
 122:          /// <param name="TRANSA">
 123:          /// - CHARACTER*1.
 124:          /// On entry, TRANSA specifies the form of op( A ) to be used in
 125:          /// the matrix multiplication as follows:
 126:          /// 
 127:          /// TRANSA = 'N' or 'n'   op( A ) = A.
 128:          /// 
 129:          /// TRANSA = 'T' or 't'   op( A ) = A'.
 130:          /// 
 131:          /// TRANSA = 'C' or 'c'   op( A ) = A'.
 132:          /// 
 133:          /// Unchanged on exit.
 134:          ///</param>
 135:          /// <param name="DIAG">
 136:          /// - CHARACTER*1.
 137:          /// On entry, DIAG specifies whether or not A is unit triangular
 138:          /// as follows:
 139:          /// 
 140:          /// DIAG = 'U' or 'u'   A is assumed to be unit triangular.
 141:          /// 
 142:          /// DIAG = 'N' or 'n'   A is not assumed to be unit
 143:          /// triangular.
 144:          /// 
 145:          /// Unchanged on exit.
 146:          ///</param>
 147:          /// <param name="M">
 148:          /// - INTEGER.
 149:          /// On entry, M specifies the number of rows of B. M must be at
 150:          /// least zero.
 151:          /// Unchanged on exit.
 152:          ///</param>
 153:          /// <param name="N">
 154:          /// - INTEGER.
 155:          /// On entry, N specifies the number of columns of B.  N must be
 156:          /// at least zero.
 157:          /// Unchanged on exit.
 158:          ///</param>
 159:          /// <param name="ALPHA">
 160:          /// - DOUBLE PRECISION.
 161:          /// On entry,  ALPHA specifies the scalar  alpha. When  alpha is
 162:          /// zero then  A is not referenced and  B need not be set before
 163:          /// entry.
 164:          /// Unchanged on exit.
 165:          ///</param>
 166:          /// <param name="A">
 167:          /// - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
 168:          /// when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
 169:          /// Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
 170:          /// upper triangular part of the array  A must contain the upper
 171:          /// triangular matrix  and the strictly lower triangular part of
 172:          /// A is not referenced.
 173:          /// Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
 174:          /// lower triangular part of the array  A must contain the lower
 175:          /// triangular matrix  and the strictly upper triangular part of
 176:          /// A is not referenced.
 177:          /// Note that when  DIAG = 'U' or 'u',  the diagonal elements of
 178:          /// A  are not referenced either,  but are assumed to be  unity.
 179:          /// Unchanged on exit.
 180:          ///</param>
 181:          /// <param name="LDA">
 182:          /// - INTEGER.
 183:          /// On entry, LDA specifies the first dimension of A as declared
 184:          /// in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
 185:          /// LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
 186:          /// then LDA must be at least max( 1, n ).
 187:          /// Unchanged on exit.
 188:          ///</param>
 189:          /// <param name="B">
 190:          /// - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
 191:          /// Before entry,  the leading  m by n part of the array  B must
 192:          /// contain  the  right-hand  side  matrix  B,  and  on exit  is
 193:          /// overwritten by the solution matrix  X.
 194:          ///</param>
 195:          /// <param name="LDB">
 196:          /// - INTEGER.
 197:          /// On entry, LDB specifies the first dimension of B as declared
 198:          /// in  the  calling  (sub)  program.   LDB  must  be  at  least
 199:          /// max( 1, m ).
 200:          /// Unchanged on exit.
 201:          /// 
 202:          ///</param>
 203:          public void Run(string SIDE, string UPLO, string TRANSA, string DIAG, int M, int N
 204:                           , double ALPHA, double[] A, int offset_a, int LDA, ref double[] B, int offset_b, int LDB)
 205:          {
 206:   
 207:              #region Array Index Correction
 208:              
 209:               int o_a = -1 - LDA + offset_a;  int o_b = -1 - LDB + offset_b; 
 210:   
 211:              #endregion
 212:   
 213:   
 214:              #region Strings
 215:              
 216:              SIDE = SIDE.Substring(0, 1);  UPLO = UPLO.Substring(0, 1);  TRANSA = TRANSA.Substring(0, 1);  
 217:              DIAG = DIAG.Substring(0, 1); 
 218:   
 219:              #endregion
 220:   
 221:   
 222:              #region Prolog
 223:              
 224:              // *     .. Scalar Arguments ..
 225:              // *     ..
 226:              // *     .. Array Arguments ..
 227:              // *     ..
 228:              // *
 229:              // *  Purpose
 230:              // *  =======
 231:              // *
 232:              // *  DTRSM  solves one of the matrix equations
 233:              // *
 234:              // *     op( A )*X = alpha*B,   or   X*op( A ) = alpha*B,
 235:              // *
 236:              // *  where alpha is a scalar, X and B are m by n matrices, A is a unit, or
 237:              // *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
 238:              // *
 239:              // *     op( A ) = A   or   op( A ) = A'.
 240:              // *
 241:              // *  The matrix X is overwritten on B.
 242:              // *
 243:              // *  Arguments
 244:              // *  ==========
 245:              // *
 246:              // *  SIDE   - CHARACTER*1.
 247:              // *           On entry, SIDE specifies whether op( A ) appears on the left
 248:              // *           or right of X as follows:
 249:              // *
 250:              // *              SIDE = 'L' or 'l'   op( A )*X = alpha*B.
 251:              // *
 252:              // *              SIDE = 'R' or 'r'   X*op( A ) = alpha*B.
 253:              // *
 254:              // *           Unchanged on exit.
 255:              // *
 256:              // *  UPLO   - CHARACTER*1.
 257:              // *           On entry, UPLO specifies whether the matrix A is an upper or
 258:              // *           lower triangular matrix as follows:
 259:              // *
 260:              // *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
 261:              // *
 262:              // *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
 263:              // *
 264:              // *           Unchanged on exit.
 265:              // *
 266:              // *  TRANSA - CHARACTER*1.
 267:              // *           On entry, TRANSA specifies the form of op( A ) to be used in
 268:              // *           the matrix multiplication as follows:
 269:              // *
 270:              // *              TRANSA = 'N' or 'n'   op( A ) = A.
 271:              // *
 272:              // *              TRANSA = 'T' or 't'   op( A ) = A'.
 273:              // *
 274:              // *              TRANSA = 'C' or 'c'   op( A ) = A'.
 275:              // *
 276:              // *           Unchanged on exit.
 277:              // *
 278:              // *  DIAG   - CHARACTER*1.
 279:              // *           On entry, DIAG specifies whether or not A is unit triangular
 280:              // *           as follows:
 281:              // *
 282:              // *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
 283:              // *
 284:              // *              DIAG = 'N' or 'n'   A is not assumed to be unit
 285:              // *                                  triangular.
 286:              // *
 287:              // *           Unchanged on exit.
 288:              // *
 289:              // *  M      - INTEGER.
 290:              // *           On entry, M specifies the number of rows of B. M must be at
 291:              // *           least zero.
 292:              // *           Unchanged on exit.
 293:              // *
 294:              // *  N      - INTEGER.
 295:              // *           On entry, N specifies the number of columns of B.  N must be
 296:              // *           at least zero.
 297:              // *           Unchanged on exit.
 298:              // *
 299:              // *  ALPHA  - DOUBLE PRECISION.
 300:              // *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
 301:              // *           zero then  A is not referenced and  B need not be set before
 302:              // *           entry.
 303:              // *           Unchanged on exit.
 304:              // *
 305:              // *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
 306:              // *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
 307:              // *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
 308:              // *           upper triangular part of the array  A must contain the upper
 309:              // *           triangular matrix  and the strictly lower triangular part of
 310:              // *           A is not referenced.
 311:              // *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
 312:              // *           lower triangular part of the array  A must contain the lower
 313:              // *           triangular matrix  and the strictly upper triangular part of
 314:              // *           A is not referenced.
 315:              // *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
 316:              // *           A  are not referenced either,  but are assumed to be  unity.
 317:              // *           Unchanged on exit.
 318:              // *
 319:              // *  LDA    - INTEGER.
 320:              // *           On entry, LDA specifies the first dimension of A as declared
 321:              // *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
 322:              // *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
 323:              // *           then LDA must be at least max( 1, n ).
 324:              // *           Unchanged on exit.
 325:              // *
 326:              // *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
 327:              // *           Before entry,  the leading  m by n part of the array  B must
 328:              // *           contain  the  right-hand  side  matrix  B,  and  on exit  is
 329:              // *           overwritten by the solution matrix  X.
 330:              // *
 331:              // *  LDB    - INTEGER.
 332:              // *           On entry, LDB specifies the first dimension of B as declared
 333:              // *           in  the  calling  (sub)  program.   LDB  must  be  at  least
 334:              // *           max( 1, m ).
 335:              // *           Unchanged on exit.
 336:              // *
 337:              // *
 338:              // *  Level 3 Blas routine.
 339:              // *
 340:              // *
 341:              // *  -- Written on 8-February-1989.
 342:              // *     Jack Dongarra, Argonne National Laboratory.
 343:              // *     Iain Duff, AERE Harwell.
 344:              // *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
 345:              // *     Sven Hammarling, Numerical Algorithms Group Ltd.
 346:              // *
 347:              // *
 348:              // *     .. External Functions ..
 349:              // *     ..
 350:              // *     .. External Subroutines ..
 351:              // *     ..
 352:              // *     .. Intrinsic Functions ..
 353:              //      INTRINSIC MAX;
 354:              // *     ..
 355:              // *     .. Local Scalars ..
 356:              // *     ..
 357:              // *     .. Parameters ..
 358:              // *     ..
 359:              // *
 360:              // *     Test the input parameters.
 361:              // *
 362:   
 363:              #endregion
 364:   
 365:   
 366:              #region Body
 367:              
 368:              LSIDE = this._lsame.Run(SIDE, "L");
 369:              if (LSIDE)
 370:              {
 371:                  NROWA = M;
 372:              }
 373:              else
 374:              {
 375:                  NROWA = N;
 376:              }
 377:              NOUNIT = this._lsame.Run(DIAG, "N");
 378:              UPPER = this._lsame.Run(UPLO, "U");
 379:              // *
 380:              INFO = 0;
 381:              if ((!LSIDE) && (!this._lsame.Run(SIDE, "R")))
 382:              {
 383:                  INFO = 1;
 384:              }
 385:              else
 386:              {
 387:                  if ((!UPPER) && (!this._lsame.Run(UPLO, "L")))
 388:                  {
 389:                      INFO = 2;
 390:                  }
 391:                  else
 392:                  {
 393:                      if ((!this._lsame.Run(TRANSA, "N")) && (!this._lsame.Run(TRANSA, "T")) && (!this._lsame.Run(TRANSA, "C")))
 394:                      {
 395:                          INFO = 3;
 396:                      }
 397:                      else
 398:                      {
 399:                          if ((!this._lsame.Run(DIAG, "U")) && (!this._lsame.Run(DIAG, "N")))
 400:                          {
 401:                              INFO = 4;
 402:                          }
 403:                          else
 404:                          {
 405:                              if (M < 0)
 406:                              {
 407:                                  INFO = 5;
 408:                              }
 409:                              else
 410:                              {
 411:                                  if (N < 0)
 412:                                  {
 413:                                      INFO = 6;
 414:                                  }
 415:                                  else
 416:                                  {
 417:                                      if (LDA < Math.Max(1, NROWA))
 418:                                      {
 419:                                          INFO = 9;
 420:                                      }
 421:                                      else
 422:                                      {
 423:                                          if (LDB < Math.Max(1, M))
 424:                                          {
 425:                                              INFO = 11;
 426:                                          }
 427:                                      }
 428:                                  }
 429:                              }
 430:                          }
 431:                      }
 432:                  }
 433:              }
 434:              if (INFO != 0)
 435:              {
 436:                  this._xerbla.Run("DTRSM ", INFO);
 437:                  return;
 438:              }
 439:              // *
 440:              // *     Quick return if possible.
 441:              // *
 442:              if (N == 0) return;
 443:              // *
 444:              // *     And when  alpha.eq.zero.
 445:              // *
 446:              if (ALPHA == ZERO)
 447:              {
 448:                  for (J = 1; J <= N; J++)
 449:                  {
 450:                      for (I = 1; I <= M; I++)
 451:                      {
 452:                          B[I+J * LDB + o_b] = ZERO;
 453:                      }
 454:                  }
 455:                  return;
 456:              }
 457:              // *
 458:              // *     Start the operations.
 459:              // *
 460:              if (LSIDE)
 461:              {
 462:                  if (this._lsame.Run(TRANSA, "N"))
 463:                  {
 464:                      // *
 465:                      // *           Form  B := alpha*inv( A )*B.
 466:                      // *
 467:                      if (UPPER)
 468:                      {
 469:                          for (J = 1; J <= N; J++)
 470:                          {
 471:                              if (ALPHA != ONE)
 472:                              {
 473:                                  for (I = 1; I <= M; I++)
 474:                                  {
 475:                                      B[I+J * LDB + o_b] = ALPHA * B[I+J * LDB + o_b];
 476:                                  }
 477:                              }
 478:                              for (K = M; K >= 1; K +=  - 1)
 479:                              {
 480:                                  if (B[K+J * LDB + o_b] != ZERO)
 481:                                  {
 482:                                      if (NOUNIT) B[K+J * LDB + o_b] = B[K+J * LDB + o_b] / A[K+K * LDA + o_a];
 483:                                      for (I = 1; I <= K - 1; I++)
 484:                                      {
 485:                                          B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - B[K+J * LDB + o_b] * A[I+K * LDA + o_a];
 486:                                      }
 487:                                  }
 488:                              }
 489:                          }
 490:                      }
 491:                      else
 492:                      {
 493:                          for (J = 1; J <= N; J++)
 494:                          {
 495:                              if (ALPHA != ONE)
 496:                              {
 497:                                  for (I = 1; I <= M; I++)
 498:                                  {
 499:                                      B[I+J * LDB + o_b] = ALPHA * B[I+J * LDB + o_b];
 500:                                  }
 501:                              }
 502:                              for (K = 1; K <= M; K++)
 503:                              {
 504:                                  if (B[K+J * LDB + o_b] != ZERO)
 505:                                  {
 506:                                      if (NOUNIT) B[K+J * LDB + o_b] = B[K+J * LDB + o_b] / A[K+K * LDA + o_a];
 507:                                      for (I = K + 1; I <= M; I++)
 508:                                      {
 509:                                          B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - B[K+J * LDB + o_b] * A[I+K * LDA + o_a];
 510:                                      }
 511:                                  }
 512:                              }
 513:                          }
 514:                      }
 515:                  }
 516:                  else
 517:                  {
 518:                      // *
 519:                      // *           Form  B := alpha*inv( A' )*B.
 520:                      // *
 521:                      if (UPPER)
 522:                      {
 523:                          for (J = 1; J <= N; J++)
 524:                          {
 525:                              for (I = 1; I <= M; I++)
 526:                              {
 527:                                  TEMP = ALPHA * B[I+J * LDB + o_b];
 528:                                  for (K = 1; K <= I - 1; K++)
 529:                                  {
 530:                                      TEMP = TEMP - A[K+I * LDA + o_a] * B[K+J * LDB + o_b];
 531:                                  }
 532:                                  if (NOUNIT) TEMP = TEMP / A[I+I * LDA + o_a];
 533:                                  B[I+J * LDB + o_b] = TEMP;
 534:                              }
 535:                          }
 536:                      }
 537:                      else
 538:                      {
 539:                          for (J = 1; J <= N; J++)
 540:                          {
 541:                              for (I = M; I >= 1; I +=  - 1)
 542:                              {
 543:                                  TEMP = ALPHA * B[I+J * LDB + o_b];
 544:                                  for (K = I + 1; K <= M; K++)
 545:                                  {
 546:                                      TEMP = TEMP - A[K+I * LDA + o_a] * B[K+J * LDB + o_b];
 547:                                  }
 548:                                  if (NOUNIT) TEMP = TEMP / A[I+I * LDA + o_a];
 549:                                  B[I+J * LDB + o_b] = TEMP;
 550:                              }
 551:                          }
 552:                      }
 553:                  }
 554:              }
 555:              else
 556:              {
 557:                  if (this._lsame.Run(TRANSA, "N"))
 558:                  {
 559:                      // *
 560:                      // *           Form  B := alpha*B*inv( A ).
 561:                      // *
 562:                      if (UPPER)
 563:                      {
 564:                          for (J = 1; J <= N; J++)
 565:                          {
 566:                              if (ALPHA != ONE)
 567:                              {
 568:                                  for (I = 1; I <= M; I++)
 569:                                  {
 570:                                      B[I+J * LDB + o_b] = ALPHA * B[I+J * LDB + o_b];
 571:                                  }
 572:                              }
 573:                              for (K = 1; K <= J - 1; K++)
 574:                              {
 575:                                  if (A[K+J * LDA + o_a] != ZERO)
 576:                                  {
 577:                                      for (I = 1; I <= M; I++)
 578:                                      {
 579:                                          B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - A[K+J * LDA + o_a] * B[I+K * LDB + o_b];
 580:                                      }
 581:                                  }
 582:                              }
 583:                              if (NOUNIT)
 584:                              {
 585:                                  TEMP = ONE / A[J+J * LDA + o_a];
 586:                                  for (I = 1; I <= M; I++)
 587:                                  {
 588:                                      B[I+J * LDB + o_b] = TEMP * B[I+J * LDB + o_b];
 589:                                  }
 590:                              }
 591:                          }
 592:                      }
 593:                      else
 594:                      {
 595:                          for (J = N; J >= 1; J +=  - 1)
 596:                          {
 597:                              if (ALPHA != ONE)
 598:                              {
 599:                                  for (I = 1; I <= M; I++)
 600:                                  {
 601:                                      B[I+J * LDB + o_b] = ALPHA * B[I+J * LDB + o_b];
 602:                                  }
 603:                              }
 604:                              for (K = J + 1; K <= N; K++)
 605:                              {
 606:                                  if (A[K+J * LDA + o_a] != ZERO)
 607:                                  {
 608:                                      for (I = 1; I <= M; I++)
 609:                                      {
 610:                                          B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - A[K+J * LDA + o_a] * B[I+K * LDB + o_b];
 611:                                      }
 612:                                  }
 613:                              }
 614:                              if (NOUNIT)
 615:                              {
 616:                                  TEMP = ONE / A[J+J * LDA + o_a];
 617:                                  for (I = 1; I <= M; I++)
 618:                                  {
 619:                                      B[I+J * LDB + o_b] = TEMP * B[I+J * LDB + o_b];
 620:                                  }
 621:                              }
 622:                          }
 623:                      }
 624:                  }
 625:                  else
 626:                  {
 627:                      // *
 628:                      // *           Form  B := alpha*B*inv( A' ).
 629:                      // *
 630:                      if (UPPER)
 631:                      {
 632:                          for (K = N; K >= 1; K +=  - 1)
 633:                          {
 634:                              if (NOUNIT)
 635:                              {
 636:                                  TEMP = ONE / A[K+K * LDA + o_a];
 637:                                  for (I = 1; I <= M; I++)
 638:                                  {
 639:                                      B[I+K * LDB + o_b] = TEMP * B[I+K * LDB + o_b];
 640:                                  }
 641:                              }
 642:                              for (J = 1; J <= K - 1; J++)
 643:                              {
 644:                                  if (A[J+K * LDA + o_a] != ZERO)
 645:                                  {
 646:                                      TEMP = A[J+K * LDA + o_a];
 647:                                      for (I = 1; I <= M; I++)
 648:                                      {
 649:                                          B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - TEMP * B[I+K * LDB + o_b];
 650:                                      }
 651:                                  }
 652:                              }
 653:                              if (ALPHA != ONE)
 654:                              {
 655:                                  for (I = 1; I <= M; I++)
 656:                                  {
 657:                                      B[I+K * LDB + o_b] = ALPHA * B[I+K * LDB + o_b];
 658:                                  }
 659:                              }
 660:                          }
 661:                      }
 662:                      else
 663:                      {
 664:                          for (K = 1; K <= N; K++)
 665:                          {
 666:                              if (NOUNIT)
 667:                              {
 668:                                  TEMP = ONE / A[K+K * LDA + o_a];
 669:                                  for (I = 1; I <= M; I++)
 670:                                  {
 671:                                      B[I+K * LDB + o_b] = TEMP * B[I+K * LDB + o_b];
 672:                                  }
 673:                              }
 674:                              for (J = K + 1; J <= N; J++)
 675:                              {
 676:                                  if (A[J+K * LDA + o_a] != ZERO)
 677:                                  {
 678:                                      TEMP = A[J+K * LDA + o_a];
 679:                                      for (I = 1; I <= M; I++)
 680:                                      {
 681:                                          B[I+J * LDB + o_b] = B[I+J * LDB + o_b] - TEMP * B[I+K * LDB + o_b];
 682:                                      }
 683:                                  }
 684:                              }
 685:                              if (ALPHA != ONE)
 686:                              {
 687:                                  for (I = 1; I <= M; I++)
 688:                                  {
 689:                                      B[I+K * LDB + o_b] = ALPHA * B[I+K * LDB + o_b];
 690:                                  }
 691:                              }
 692:                          }
 693:                      }
 694:                  }
 695:              }
 696:              // *
 697:              return;
 698:              // *
 699:              // *     End of DTRSM .
 700:              // *
 701:   
 702:              #endregion
 703:   
 704:          }
 705:      }
 706:  }