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:      /// DTRSV  solves one of the systems of equations
  24:      /// 
  25:      /// A*x = b,   or   A'*x = b,
  26:      /// 
  27:      /// where b and x are n element vectors and A is an n by n unit, or
  28:      /// non-unit, upper or lower triangular matrix.
  29:      /// 
  30:      /// No test for singularity or near-singularity is included in this
  31:      /// routine. Such tests must be performed before calling this routine.
  32:      /// 
  33:      ///</summary>
  34:      public class DTRSV
  35:      {
  36:      
  37:   
  38:          #region Dependencies
  39:          
  40:          LSAME _lsame; XERBLA _xerbla; 
  41:   
  42:          #endregion
  43:   
  44:   
  45:          #region Fields
  46:          
  47:          const double ZERO = 0.0E+0; double TEMP = 0; int I = 0; int INFO = 0; int IX = 0; int J = 0; int JX = 0; int KX = 0; 
  48:          bool NOUNIT = false;
  49:   
  50:          #endregion
  51:   
  52:          public DTRSV(LSAME lsame, XERBLA xerbla)
  53:          {
  54:      
  55:   
  56:              #region Set Dependencies
  57:              
  58:              this._lsame = lsame; this._xerbla = xerbla; 
  59:   
  60:              #endregion
  61:   
  62:          }
  63:      
  64:          public DTRSV()
  65:          {
  66:      
  67:   
  68:              #region Dependencies (Initialization)
  69:              
  70:              LSAME lsame = new LSAME();
  71:              XERBLA xerbla = new XERBLA();
  72:   
  73:              #endregion
  74:   
  75:   
  76:              #region Set Dependencies
  77:              
  78:              this._lsame = lsame; this._xerbla = xerbla; 
  79:   
  80:              #endregion
  81:   
  82:          }
  83:          /// <summary>
  84:          /// Purpose
  85:          /// =======
  86:          /// 
  87:          /// DTRSV  solves one of the systems of equations
  88:          /// 
  89:          /// A*x = b,   or   A'*x = b,
  90:          /// 
  91:          /// where b and x are n element vectors and A is an n by n unit, or
  92:          /// non-unit, upper or lower triangular matrix.
  93:          /// 
  94:          /// No test for singularity or near-singularity is included in this
  95:          /// routine. Such tests must be performed before calling this routine.
  96:          /// 
  97:          ///</summary>
  98:          /// <param name="UPLO">
  99:          /// - CHARACTER*1.
 100:          /// On entry, UPLO specifies whether the matrix is an upper or
 101:          /// lower triangular matrix as follows:
 102:          /// 
 103:          /// UPLO = 'U' or 'u'   A is an upper triangular matrix.
 104:          /// 
 105:          /// UPLO = 'L' or 'l'   A is a lower triangular matrix.
 106:          /// 
 107:          /// Unchanged on exit.
 108:          ///</param>
 109:          /// <param name="TRANS">
 110:          /// - CHARACTER*1.
 111:          /// On entry, TRANS specifies the equations to be solved as
 112:          /// follows:
 113:          /// 
 114:          /// TRANS = 'N' or 'n'   A*x = b.
 115:          /// 
 116:          /// TRANS = 'T' or 't'   A'*x = b.
 117:          /// 
 118:          /// TRANS = 'C' or 'c'   A'*x = b.
 119:          /// 
 120:          /// Unchanged on exit.
 121:          ///</param>
 122:          /// <param name="DIAG">
 123:          /// - CHARACTER*1.
 124:          /// On entry, DIAG specifies whether or not A is unit
 125:          /// triangular as follows:
 126:          /// 
 127:          /// DIAG = 'U' or 'u'   A is assumed to be unit triangular.
 128:          /// 
 129:          /// DIAG = 'N' or 'n'   A is not assumed to be unit
 130:          /// triangular.
 131:          /// 
 132:          /// Unchanged on exit.
 133:          ///</param>
 134:          /// <param name="N">
 135:          /// - INTEGER.
 136:          /// On entry, N specifies the order of the matrix A.
 137:          /// N must be at least zero.
 138:          /// Unchanged on exit.
 139:          ///</param>
 140:          /// <param name="A">
 141:          /// - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
 142:          /// Before entry with  UPLO = 'U' or 'u', the leading n by n
 143:          /// upper triangular part of the array A must contain the upper
 144:          /// triangular matrix and the strictly lower triangular part of
 145:          /// A is not referenced.
 146:          /// Before entry with UPLO = 'L' or 'l', the leading n by n
 147:          /// lower triangular part of the array A must contain the lower
 148:          /// triangular matrix and the strictly upper triangular part of
 149:          /// A is not referenced.
 150:          /// Note that when  DIAG = 'U' or 'u', the diagonal elements of
 151:          /// A are not referenced either, but are assumed to be unity.
 152:          /// Unchanged on exit.
 153:          ///</param>
 154:          /// <param name="LDA">
 155:          /// - INTEGER.
 156:          /// On entry, LDA specifies the first dimension of A as declared
 157:          /// in the calling (sub) program. LDA must be at least
 158:          /// max( 1, n ).
 159:          /// Unchanged on exit.
 160:          ///</param>
 161:          /// <param name="X">
 162:          /// - DOUBLE PRECISION array of dimension at least
 163:          /// ( 1 + ( n - 1 )*abs( INCX ) ).
 164:          /// Before entry, the incremented array X must contain the n
 165:          /// element right-hand side vector b. On exit, X is overwritten
 166:          /// with the solution vector x.
 167:          ///</param>
 168:          /// <param name="INCX">
 169:          /// - INTEGER.
 170:          /// On entry, INCX specifies the increment for the elements of
 171:          /// X. INCX must not be zero.
 172:          /// Unchanged on exit.
 173:          /// 
 174:          ///</param>
 175:          public void Run(string UPLO, string TRANS, string DIAG, int N, double[] A, int offset_a, int LDA
 176:                           , ref double[] X, int offset_x, int INCX)
 177:          {
 178:   
 179:              #region Array Index Correction
 180:              
 181:               int o_a = -1 - LDA + offset_a;  int o_x = -1 + offset_x; 
 182:   
 183:              #endregion
 184:   
 185:   
 186:              #region Strings
 187:              
 188:              UPLO = UPLO.Substring(0, 1);  TRANS = TRANS.Substring(0, 1);  DIAG = DIAG.Substring(0, 1);  
 189:   
 190:              #endregion
 191:   
 192:   
 193:              #region Prolog
 194:              
 195:              // *     .. Scalar Arguments ..
 196:              // *     ..
 197:              // *     .. Array Arguments ..
 198:              // *     ..
 199:              // *
 200:              // *  Purpose
 201:              // *  =======
 202:              // *
 203:              // *  DTRSV  solves one of the systems of equations
 204:              // *
 205:              // *     A*x = b,   or   A'*x = b,
 206:              // *
 207:              // *  where b and x are n element vectors and A is an n by n unit, or
 208:              // *  non-unit, upper or lower triangular matrix.
 209:              // *
 210:              // *  No test for singularity or near-singularity is included in this
 211:              // *  routine. Such tests must be performed before calling this routine.
 212:              // *
 213:              // *  Arguments
 214:              // *  ==========
 215:              // *
 216:              // *  UPLO   - CHARACTER*1.
 217:              // *           On entry, UPLO specifies whether the matrix is an upper or
 218:              // *           lower triangular matrix as follows:
 219:              // *
 220:              // *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
 221:              // *
 222:              // *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
 223:              // *
 224:              // *           Unchanged on exit.
 225:              // *
 226:              // *  TRANS  - CHARACTER*1.
 227:              // *           On entry, TRANS specifies the equations to be solved as
 228:              // *           follows:
 229:              // *
 230:              // *              TRANS = 'N' or 'n'   A*x = b.
 231:              // *
 232:              // *              TRANS = 'T' or 't'   A'*x = b.
 233:              // *
 234:              // *              TRANS = 'C' or 'c'   A'*x = b.
 235:              // *
 236:              // *           Unchanged on exit.
 237:              // *
 238:              // *  DIAG   - CHARACTER*1.
 239:              // *           On entry, DIAG specifies whether or not A is unit
 240:              // *           triangular as follows:
 241:              // *
 242:              // *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
 243:              // *
 244:              // *              DIAG = 'N' or 'n'   A is not assumed to be unit
 245:              // *                                  triangular.
 246:              // *
 247:              // *           Unchanged on exit.
 248:              // *
 249:              // *  N      - INTEGER.
 250:              // *           On entry, N specifies the order of the matrix A.
 251:              // *           N must be at least zero.
 252:              // *           Unchanged on exit.
 253:              // *
 254:              // *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
 255:              // *           Before entry with  UPLO = 'U' or 'u', the leading n by n
 256:              // *           upper triangular part of the array A must contain the upper
 257:              // *           triangular matrix and the strictly lower triangular part of
 258:              // *           A is not referenced.
 259:              // *           Before entry with UPLO = 'L' or 'l', the leading n by n
 260:              // *           lower triangular part of the array A must contain the lower
 261:              // *           triangular matrix and the strictly upper triangular part of
 262:              // *           A is not referenced.
 263:              // *           Note that when  DIAG = 'U' or 'u', the diagonal elements of
 264:              // *           A are not referenced either, but are assumed to be unity.
 265:              // *           Unchanged on exit.
 266:              // *
 267:              // *  LDA    - INTEGER.
 268:              // *           On entry, LDA specifies the first dimension of A as declared
 269:              // *           in the calling (sub) program. LDA must be at least
 270:              // *           max( 1, n ).
 271:              // *           Unchanged on exit.
 272:              // *
 273:              // *  X      - DOUBLE PRECISION array of dimension at least
 274:              // *           ( 1 + ( n - 1 )*abs( INCX ) ).
 275:              // *           Before entry, the incremented array X must contain the n
 276:              // *           element right-hand side vector b. On exit, X is overwritten
 277:              // *           with the solution vector x.
 278:              // *
 279:              // *  INCX   - INTEGER.
 280:              // *           On entry, INCX specifies the increment for the elements of
 281:              // *           X. INCX must not be zero.
 282:              // *           Unchanged on exit.
 283:              // *
 284:              // *
 285:              // *  Level 2 Blas routine.
 286:              // *
 287:              // *  -- Written on 22-October-1986.
 288:              // *     Jack Dongarra, Argonne National Lab.
 289:              // *     Jeremy Du Croz, Nag Central Office.
 290:              // *     Sven Hammarling, Nag Central Office.
 291:              // *     Richard Hanson, Sandia National Labs.
 292:              // *
 293:              // *
 294:              // *     .. Parameters ..
 295:              // *     ..
 296:              // *     .. Local Scalars ..
 297:              // *     ..
 298:              // *     .. External Functions ..
 299:              // *     ..
 300:              // *     .. External Subroutines ..
 301:              // *     ..
 302:              // *     .. Intrinsic Functions ..
 303:              //      INTRINSIC MAX;
 304:              // *     ..
 305:              // *
 306:              // *     Test the input parameters.
 307:              // *
 308:   
 309:              #endregion
 310:   
 311:   
 312:              #region Body
 313:              
 314:              INFO = 0;
 315:              if (!this._lsame.Run(UPLO, "U") && !this._lsame.Run(UPLO, "L"))
 316:              {
 317:                  INFO = 1;
 318:              }
 319:              else
 320:              {
 321:                  if (!this._lsame.Run(TRANS, "N") && !this._lsame.Run(TRANS, "T") && !this._lsame.Run(TRANS, "C"))
 322:                  {
 323:                      INFO = 2;
 324:                  }
 325:                  else
 326:                  {
 327:                      if (!this._lsame.Run(DIAG, "U") && !this._lsame.Run(DIAG, "N"))
 328:                      {
 329:                          INFO = 3;
 330:                      }
 331:                      else
 332:                      {
 333:                          if (N < 0)
 334:                          {
 335:                              INFO = 4;
 336:                          }
 337:                          else
 338:                          {
 339:                              if (LDA < Math.Max(1, N))
 340:                              {
 341:                                  INFO = 6;
 342:                              }
 343:                              else
 344:                              {
 345:                                  if (INCX == 0)
 346:                                  {
 347:                                      INFO = 8;
 348:                                  }
 349:                              }
 350:                          }
 351:                      }
 352:                  }
 353:              }
 354:              if (INFO != 0)
 355:              {
 356:                  this._xerbla.Run("DTRSV ", INFO);
 357:                  return;
 358:              }
 359:              // *
 360:              // *     Quick return if possible.
 361:              // *
 362:              if (N == 0) return;
 363:              // *
 364:              NOUNIT = this._lsame.Run(DIAG, "N");
 365:              // *
 366:              // *     Set up the start point in X if the increment is not unity. This
 367:              // *     will be  ( N - 1 )*INCX  too small for descending loops.
 368:              // *
 369:              if (INCX <= 0)
 370:              {
 371:                  KX = 1 - (N - 1) * INCX;
 372:              }
 373:              else
 374:              {
 375:                  if (INCX != 1)
 376:                  {
 377:                      KX = 1;
 378:                  }
 379:              }
 380:              // *
 381:              // *     Start the operations. In this version the elements of A are
 382:              // *     accessed sequentially with one pass through A.
 383:              // *
 384:              if (this._lsame.Run(TRANS, "N"))
 385:              {
 386:                  // *
 387:                  // *        Form  x := inv( A )*x.
 388:                  // *
 389:                  if (this._lsame.Run(UPLO, "U"))
 390:                  {
 391:                      if (INCX == 1)
 392:                      {
 393:                          for (J = N; J >= 1; J +=  - 1)
 394:                          {
 395:                              if (X[J + o_x] != ZERO)
 396:                              {
 397:                                  if (NOUNIT) X[J + o_x] = X[J + o_x] / A[J+J * LDA + o_a];
 398:                                  TEMP = X[J + o_x];
 399:                                  for (I = J - 1; I >= 1; I +=  - 1)
 400:                                  {
 401:                                      X[I + o_x] = X[I + o_x] - TEMP * A[I+J * LDA + o_a];
 402:                                  }
 403:                              }
 404:                          }
 405:                      }
 406:                      else
 407:                      {
 408:                          JX = KX + (N - 1) * INCX;
 409:                          for (J = N; J >= 1; J +=  - 1)
 410:                          {
 411:                              if (X[JX + o_x] != ZERO)
 412:                              {
 413:                                  if (NOUNIT) X[JX + o_x] = X[JX + o_x] / A[J+J * LDA + o_a];
 414:                                  TEMP = X[JX + o_x];
 415:                                  IX = JX;
 416:                                  for (I = J - 1; I >= 1; I +=  - 1)
 417:                                  {
 418:                                      IX = IX - INCX;
 419:                                      X[IX + o_x] = X[IX + o_x] - TEMP * A[I+J * LDA + o_a];
 420:                                  }
 421:                              }
 422:                              JX = JX - INCX;
 423:                          }
 424:                      }
 425:                  }
 426:                  else
 427:                  {
 428:                      if (INCX == 1)
 429:                      {
 430:                          for (J = 1; J <= N; J++)
 431:                          {
 432:                              if (X[J + o_x] != ZERO)
 433:                              {
 434:                                  if (NOUNIT) X[J + o_x] = X[J + o_x] / A[J+J * LDA + o_a];
 435:                                  TEMP = X[J + o_x];
 436:                                  for (I = J + 1; I <= N; I++)
 437:                                  {
 438:                                      X[I + o_x] = X[I + o_x] - TEMP * A[I+J * LDA + o_a];
 439:                                  }
 440:                              }
 441:                          }
 442:                      }
 443:                      else
 444:                      {
 445:                          JX = KX;
 446:                          for (J = 1; J <= N; J++)
 447:                          {
 448:                              if (X[JX + o_x] != ZERO)
 449:                              {
 450:                                  if (NOUNIT) X[JX + o_x] = X[JX + o_x] / A[J+J * LDA + o_a];
 451:                                  TEMP = X[JX + o_x];
 452:                                  IX = JX;
 453:                                  for (I = J + 1; I <= N; I++)
 454:                                  {
 455:                                      IX = IX + INCX;
 456:                                      X[IX + o_x] = X[IX + o_x] - TEMP * A[I+J * LDA + o_a];
 457:                                  }
 458:                              }
 459:                              JX = JX + INCX;
 460:                          }
 461:                      }
 462:                  }
 463:              }
 464:              else
 465:              {
 466:                  // *
 467:                  // *        Form  x := inv( A' )*x.
 468:                  // *
 469:                  if (this._lsame.Run(UPLO, "U"))
 470:                  {
 471:                      if (INCX == 1)
 472:                      {
 473:                          for (J = 1; J <= N; J++)
 474:                          {
 475:                              TEMP = X[J + o_x];
 476:                              for (I = 1; I <= J - 1; I++)
 477:                              {
 478:                                  TEMP = TEMP - A[I+J * LDA + o_a] * X[I + o_x];
 479:                              }
 480:                              if (NOUNIT) TEMP = TEMP / A[J+J * LDA + o_a];
 481:                              X[J + o_x] = TEMP;
 482:                          }
 483:                      }
 484:                      else
 485:                      {
 486:                          JX = KX;
 487:                          for (J = 1; J <= N; J++)
 488:                          {
 489:                              TEMP = X[JX + o_x];
 490:                              IX = KX;
 491:                              for (I = 1; I <= J - 1; I++)
 492:                              {
 493:                                  TEMP = TEMP - A[I+J * LDA + o_a] * X[IX + o_x];
 494:                                  IX = IX + INCX;
 495:                              }
 496:                              if (NOUNIT) TEMP = TEMP / A[J+J * LDA + o_a];
 497:                              X[JX + o_x] = TEMP;
 498:                              JX = JX + INCX;
 499:                          }
 500:                      }
 501:                  }
 502:                  else
 503:                  {
 504:                      if (INCX == 1)
 505:                      {
 506:                          for (J = N; J >= 1; J +=  - 1)
 507:                          {
 508:                              TEMP = X[J + o_x];
 509:                              for (I = N; I >= J + 1; I +=  - 1)
 510:                              {
 511:                                  TEMP = TEMP - A[I+J * LDA + o_a] * X[I + o_x];
 512:                              }
 513:                              if (NOUNIT) TEMP = TEMP / A[J+J * LDA + o_a];
 514:                              X[J + o_x] = TEMP;
 515:                          }
 516:                      }
 517:                      else
 518:                      {
 519:                          KX = KX + (N - 1) * INCX;
 520:                          JX = KX;
 521:                          for (J = N; J >= 1; J +=  - 1)
 522:                          {
 523:                              TEMP = X[JX + o_x];
 524:                              IX = KX;
 525:                              for (I = N; I >= J + 1; I +=  - 1)
 526:                              {
 527:                                  TEMP = TEMP - A[I+J * LDA + o_a] * X[IX + o_x];
 528:                                  IX = IX - INCX;
 529:                              }
 530:                              if (NOUNIT) TEMP = TEMP / A[J+J * LDA + o_a];
 531:                              X[JX + o_x] = TEMP;
 532:                              JX = JX - INCX;
 533:                          }
 534:                      }
 535:                  }
 536:              }
 537:              // *
 538:              return;
 539:              // *
 540:              // *     End of DTRSV .
 541:              // *
 542:   
 543:              #endregion
 544:   
 545:          }
 546:      }
 547:  }