`   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.0) --`
`  21:      /// Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,`
`  22:      /// Courant Institute, Argonne National Lab, and Rice University`
`  23:      /// June 30, 1992`
`  24:      /// Purpose`
`  25:      /// =======`
`  26:      /// `
`  27:      /// DLATRS solves one of the triangular systems`
`  28:      /// `
`  29:      /// A *x = s*b  or  A'*x = s*b`
`  30:      /// `
`  31:      /// with scaling to prevent overflow.  Here A is an upper or lower`
`  32:      /// triangular matrix, A' denotes the transpose of A, x and b are`
`  33:      /// n-element vectors, and s is a scaling factor, usually less than`
`  34:      /// or equal to 1, chosen so that the components of x will be less than`
`  35:      /// the overflow threshold.  If the unscaled problem will not cause`
`  36:      /// overflow, the Level 2 BLAS routine DTRSV is called.  If the matrix A`
`  37:      /// is singular (A(j,j) = 0 for some j), then s is set to 0 and a`
`  38:      /// non-trivial solution to A*x = 0 is returned.`
`  39:      /// `
`  40:      ///</summary>`
`  41:      public class DLATRS`
`  42:      {`
`  43:      `
`  44:   `
`  45:          #region Dependencies`
`  46:          `
`  47:          LSAME _lsame; IDAMAX _idamax; DASUM _dasum; DDOT _ddot; DLAMCH _dlamch; DAXPY _daxpy; DSCAL _dscal; DTRSV _dtrsv; `
`  48:          XERBLA _xerbla;`
`  49:   `
`  50:          #endregion`
`  51:   `
`  52:   `
`  53:          #region Fields`
`  54:          `
`  55:          const double ZERO = 0.0E+0; const double HALF = 0.5E+0; const double ONE = 1.0E+0; bool NOTRAN = false; `
`  56:          bool NOUNIT = false;bool UPPER = false; int I = 0; int IMAX = 0; int J = 0; int JFIRST = 0; int JINC = 0; int JLAST = 0; `
`  57:          double BIGNUM = 0;double GROW = 0; double REC = 0; double SMLNUM = 0; double SUMJ = 0; double TJJ = 0; double TJJS = 0; `
`  58:          double TMAX = 0;double TSCAL = 0; double USCAL = 0; double XBND = 0; double XJ = 0; double XMAX = 0; `
`  59:   `
`  60:          #endregion`
`  61:   `
`  62:          public DLATRS(LSAME lsame, IDAMAX idamax, DASUM dasum, DDOT ddot, DLAMCH dlamch, DAXPY daxpy, DSCAL dscal, DTRSV dtrsv, XERBLA xerbla)`
`  63:          {`
`  64:      `
`  65:   `
`  66:              #region Set Dependencies`
`  67:              `
`  68:              this._lsame = lsame; this._idamax = idamax; this._dasum = dasum; this._ddot = ddot; this._dlamch = dlamch; `
`  69:              this._daxpy = daxpy;this._dscal = dscal; this._dtrsv = dtrsv; this._xerbla = xerbla; `
`  70:   `
`  71:              #endregion`
`  72:   `
`  73:          }`
`  74:      `
`  75:          public DLATRS()`
`  76:          {`
`  77:      `
`  78:   `
`  79:              #region Dependencies (Initialization)`
`  80:              `
`  81:              LSAME lsame = new LSAME();`
`  82:              IDAMAX idamax = new IDAMAX();`
`  83:              DASUM dasum = new DASUM();`
`  84:              DDOT ddot = new DDOT();`
`  85:              DLAMC3 dlamc3 = new DLAMC3();`
`  86:              DAXPY daxpy = new DAXPY();`
`  87:              DSCAL dscal = new DSCAL();`
`  88:              XERBLA xerbla = new XERBLA();`
`  89:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  90:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  91:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  92:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
`  93:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
`  94:              DTRSV dtrsv = new DTRSV(lsame, xerbla);`
`  95:   `
`  96:              #endregion`
`  97:   `
`  98:   `
`  99:              #region Set Dependencies`
` 100:              `
` 101:              this._lsame = lsame; this._idamax = idamax; this._dasum = dasum; this._ddot = ddot; this._dlamch = dlamch; `
` 102:              this._daxpy = daxpy;this._dscal = dscal; this._dtrsv = dtrsv; this._xerbla = xerbla; `
` 103:   `
` 104:              #endregion`
` 105:   `
` 106:          }`
` 107:          /// <summary>`
` 108:          /// Purpose`
` 109:          /// =======`
` 110:          /// `
` 111:          /// DLATRS solves one of the triangular systems`
` 112:          /// `
` 113:          /// A *x = s*b  or  A'*x = s*b`
` 114:          /// `
` 115:          /// with scaling to prevent overflow.  Here A is an upper or lower`
` 116:          /// triangular matrix, A' denotes the transpose of A, x and b are`
` 117:          /// n-element vectors, and s is a scaling factor, usually less than`
` 118:          /// or equal to 1, chosen so that the components of x will be less than`
` 119:          /// the overflow threshold.  If the unscaled problem will not cause`
` 120:          /// overflow, the Level 2 BLAS routine DTRSV is called.  If the matrix A`
` 121:          /// is singular (A(j,j) = 0 for some j), then s is set to 0 and a`
` 122:          /// non-trivial solution to A*x = 0 is returned.`
` 123:          /// `
` 124:          ///</summary>`
` 125:          /// <param name="UPLO">`
` 126:          /// (input) CHARACTER*1`
` 127:          /// Specifies whether the matrix A is upper or lower triangular.`
` 128:          /// = 'U':  Upper triangular`
` 129:          /// = 'L':  Lower triangular`
` 130:          ///</param>`
` 131:          /// <param name="TRANS">`
` 132:          /// (input) CHARACTER*1`
` 133:          /// Specifies the operation applied to A.`
` 134:          /// = 'N':  Solve A * x = s*b  (No transpose)`
` 135:          /// = 'T':  Solve A'* x = s*b  (Transpose)`
` 136:          /// = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)`
` 137:          ///</param>`
` 138:          /// <param name="DIAG">`
` 139:          /// (input) CHARACTER*1`
` 140:          /// Specifies whether or not the matrix A is unit triangular.`
` 141:          /// = 'N':  Non-unit triangular`
` 142:          /// = 'U':  Unit triangular`
` 143:          ///</param>`
` 144:          /// <param name="NORMIN">`
` 145:          /// (input) CHARACTER*1`
` 146:          /// Specifies whether CNORM has been set or not.`
` 147:          /// = 'Y':  CNORM contains the column norms on entry`
` 148:          /// = 'N':  CNORM is not set on entry.  On exit, the norms will`
` 149:          /// be computed and stored in CNORM.`
` 150:          ///</param>`
` 151:          /// <param name="N">`
` 152:          /// (input) INTEGER`
` 153:          /// The order of the matrix A.  N .GE. 0.`
` 154:          ///</param>`
` 155:          /// <param name="A">`
` 156:          /// *x = s*b  or  A'*x = s*b`
` 157:          ///</param>`
` 158:          /// <param name="LDA">`
` 159:          /// (input) INTEGER`
` 160:          /// The leading dimension of the array A.  LDA .GE. max (1,N).`
` 161:          ///</param>`
` 162:          /// <param name="X">`
` 163:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 164:          /// On entry, the right hand side b of the triangular system.`
` 165:          /// On exit, X is overwritten by the solution vector x.`
` 166:          ///</param>`
` 167:          /// <param name="SCALE">`
` 168:          /// (output) DOUBLE PRECISION`
` 169:          /// The scaling factor s for the triangular system`
` 170:          /// A * x = s*b  or  A'* x = s*b.`
` 171:          /// If SCALE = 0, the matrix A is singular or badly scaled, and`
` 172:          /// the vector x is an exact or approximate solution to A*x = 0.`
` 173:          ///</param>`
` 174:          /// <param name="CNORM">`
` 175:          /// (input or output) DOUBLE PRECISION array, dimension (N)`
` 176:          /// `
` 177:          /// If NORMIN = 'Y', CNORM is an input argument and CNORM(j)`
` 178:          /// contains the norm of the off-diagonal part of the j-th column`
` 179:          /// of A.  If TRANS = 'N', CNORM(j) must be greater than or equal`
` 180:          /// to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)`
` 181:          /// must be greater than or equal to the 1-norm.`
` 182:          /// `
` 183:          /// If NORMIN = 'N', CNORM is an output argument and CNORM(j)`
` 184:          /// returns the 1-norm of the offdiagonal part of the j-th column`
` 185:          /// of A.`
` 186:          ///</param>`
` 187:          /// <param name="INFO">`
` 188:          /// (output) INTEGER`
` 189:          /// = 0:  successful exit`
` 190:          /// .LT. 0:  if INFO = -k, the k-th argument had an illegal value`
` 191:          ///</param>`
` 192:          public void Run(string UPLO, string TRANS, string DIAG, string NORMIN, int N, double[] A, int offset_a`
` 193:                           , int LDA, ref double[] X, int offset_x, ref double SCALE, ref double[] CNORM, int offset_cnorm, ref int INFO)`
` 194:          {`
` 195:   `
` 196:              #region Array Index Correction`
` 197:              `
` 198:               int o_a = -1 - LDA + offset_a;  int o_x = -1 + offset_x;  int o_cnorm = -1 + offset_cnorm; `
` 199:   `
` 200:              #endregion`
` 201:   `
` 202:   `
` 203:              #region Strings`
` 204:              `
` 205:              UPLO = UPLO.Substring(0, 1);  TRANS = TRANS.Substring(0, 1);  DIAG = DIAG.Substring(0, 1);  `
` 206:              NORMIN = NORMIN.Substring(0, 1); `
` 207:   `
` 208:              #endregion`
` 209:   `
` 210:   `
` 211:              #region Prolog`
` 212:              `
` 213:              // *`
` 214:              // *  -- LAPACK auxiliary routine (version 3.0) --`
` 215:              // *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,`
` 216:              // *     Courant Institute, Argonne National Lab, and Rice University`
` 217:              // *     June 30, 1992`
` 218:              // *`
` 219:              // *     .. Scalar Arguments ..`
` 220:              // *     ..`
` 221:              // *     .. Array Arguments ..`
` 222:              // *     ..`
` 223:              // *`
` 224:              // *  Purpose`
` 225:              // *  =======`
` 226:              // *`
` 227:              // *  DLATRS solves one of the triangular systems`
` 228:              // *`
` 229:              // *     A *x = s*b  or  A'*x = s*b`
` 230:              // *`
` 231:              // *  with scaling to prevent overflow.  Here A is an upper or lower`
` 232:              // *  triangular matrix, A' denotes the transpose of A, x and b are`
` 233:              // *  n-element vectors, and s is a scaling factor, usually less than`
` 234:              // *  or equal to 1, chosen so that the components of x will be less than`
` 235:              // *  the overflow threshold.  If the unscaled problem will not cause`
` 236:              // *  overflow, the Level 2 BLAS routine DTRSV is called.  If the matrix A`
` 237:              // *  is singular (A(j,j) = 0 for some j), then s is set to 0 and a`
` 238:              // *  non-trivial solution to A*x = 0 is returned.`
` 239:              // *`
` 240:              // *  Arguments`
` 241:              // *  =========`
` 242:              // *`
` 243:              // *  UPLO    (input) CHARACTER*1`
` 244:              // *          Specifies whether the matrix A is upper or lower triangular.`
` 245:              // *          = 'U':  Upper triangular`
` 246:              // *          = 'L':  Lower triangular`
` 247:              // *`
` 248:              // *  TRANS   (input) CHARACTER*1`
` 249:              // *          Specifies the operation applied to A.`
` 250:              // *          = 'N':  Solve A * x = s*b  (No transpose)`
` 251:              // *          = 'T':  Solve A'* x = s*b  (Transpose)`
` 252:              // *          = 'C':  Solve A'* x = s*b  (Conjugate transpose = Transpose)`
` 253:              // *`
` 254:              // *  DIAG    (input) CHARACTER*1`
` 255:              // *          Specifies whether or not the matrix A is unit triangular.`
` 256:              // *          = 'N':  Non-unit triangular`
` 257:              // *          = 'U':  Unit triangular`
` 258:              // *`
` 259:              // *  NORMIN  (input) CHARACTER*1`
` 260:              // *          Specifies whether CNORM has been set or not.`
` 261:              // *          = 'Y':  CNORM contains the column norms on entry`
` 262:              // *          = 'N':  CNORM is not set on entry.  On exit, the norms will`
` 263:              // *                  be computed and stored in CNORM.`
` 264:              // *`
` 265:              // *  N       (input) INTEGER`
` 266:              // *          The order of the matrix A.  N >= 0.`
` 267:              // *`
` 268:              // *  A       (input) DOUBLE PRECISION array, dimension (LDA,N)`
` 269:              // *          The triangular matrix A.  If UPLO = 'U', the leading n by n`
` 270:              // *          upper triangular part of the array A contains the upper`
` 271:              // *          triangular matrix, and the strictly lower triangular part of`
` 272:              // *          A is not referenced.  If UPLO = 'L', the leading n by n lower`
` 273:              // *          triangular part of the array A contains the lower triangular`
` 274:              // *          matrix, and the strictly upper triangular part of A is not`
` 275:              // *          referenced.  If DIAG = 'U', the diagonal elements of A are`
` 276:              // *          also not referenced and are assumed to be 1.`
` 277:              // *`
` 278:              // *  LDA     (input) INTEGER`
` 279:              // *          The leading dimension of the array A.  LDA >= max (1,N).`
` 280:              // *`
` 281:              // *  X       (input/output) DOUBLE PRECISION array, dimension (N)`
` 282:              // *          On entry, the right hand side b of the triangular system.`
` 283:              // *          On exit, X is overwritten by the solution vector x.`
` 284:              // *`
` 285:              // *  SCALE   (output) DOUBLE PRECISION`
` 286:              // *          The scaling factor s for the triangular system`
` 287:              // *             A * x = s*b  or  A'* x = s*b.`
` 288:              // *          If SCALE = 0, the matrix A is singular or badly scaled, and`
` 289:              // *          the vector x is an exact or approximate solution to A*x = 0.`
` 290:              // *`
` 291:              // *  CNORM   (input or output) DOUBLE PRECISION array, dimension (N)`
` 292:              // *`
` 293:              // *          If NORMIN = 'Y', CNORM is an input argument and CNORM(j)`
` 294:              // *          contains the norm of the off-diagonal part of the j-th column`
` 295:              // *          of A.  If TRANS = 'N', CNORM(j) must be greater than or equal`
` 296:              // *          to the infinity-norm, and if TRANS = 'T' or 'C', CNORM(j)`
` 297:              // *          must be greater than or equal to the 1-norm.`
` 298:              // *`
` 299:              // *          If NORMIN = 'N', CNORM is an output argument and CNORM(j)`
` 300:              // *          returns the 1-norm of the offdiagonal part of the j-th column`
` 301:              // *          of A.`
` 302:              // *`
` 303:              // *  INFO    (output) INTEGER`
` 304:              // *          = 0:  successful exit`
` 305:              // *          < 0:  if INFO = -k, the k-th argument had an illegal value`
` 306:              // *`
` 307:              // *  Further Details`
` 308:              // *  ======= =======`
` 309:              // *`
` 310:              // *  A rough bound on x is computed; if that is less than overflow, DTRSV`
` 311:              // *  is called, otherwise, specific code is used which checks for possible`
` 312:              // *  overflow or divide-by-zero at every operation.`
` 313:              // *`
` 314:              // *  A columnwise scheme is used for solving A*x = b.  The basic algorithm`
` 315:              // *  if A is lower triangular is`
` 316:              // *`
` 317:              // *       x[1:n] := b[1:n]`
` 318:              // *       for j = 1, ..., n`
` 319:              // *            x(j) := x(j) / A(j,j)`
` 320:              // *            x[j+1:n] := x[j+1:n] - x(j) * A[j+1:n,j]`
` 321:              // *       end`
` 322:              // *`
` 323:              // *  Define bounds on the components of x after j iterations of the loop:`
` 324:              // *     M(j) = bound on x[1:j]`
` 325:              // *     G(j) = bound on x[j+1:n]`
` 326:              // *  Initially, let M(0) = 0 and G(0) = max{x(i), i=1,...,n}.`
` 327:              // *`
` 328:              // *  Then for iteration j+1 we have`
` 329:              // *     M(j+1) <= G(j) / | A(j+1,j+1) |`
` 330:              // *     G(j+1) <= G(j) + M(j+1) * | A[j+2:n,j+1] |`
` 331:              // *            <= G(j) ( 1 + CNORM(j+1) / | A(j+1,j+1) | )`
` 332:              // *`
` 333:              // *  where CNORM(j+1) is greater than or equal to the infinity-norm of`
` 334:              // *  column j+1 of A, not counting the diagonal.  Hence`
` 335:              // *`
` 336:              // *     G(j) <= G(0) product ( 1 + CNORM(i) / | A(i,i) | )`
` 337:              // *                  1<=i<=j`
` 338:              // *  and`
` 339:              // *`
` 340:              // *     |x(j)| <= ( G(0) / |A(j,j)| ) product ( 1 + CNORM(i) / |A(i,i)| )`
` 341:              // *                                   1<=i< j`
` 342:              // *`
` 343:              // *  Since |x(j)| <= M(j), we use the Level 2 BLAS routine DTRSV if the`
` 344:              // *  reciprocal of the largest M(j), j=1,..,n, is larger than`
` 345:              // *  max(underflow, 1/overflow).`
` 346:              // *`
` 347:              // *  The bound on x(j) is also used to determine when a step in the`
` 348:              // *  columnwise method can be performed without fear of overflow.  If`
` 349:              // *  the computed bound is greater than a large constant, x is scaled to`
` 350:              // *  prevent overflow, but if the bound overflows, x is set to 0, x(j) to`
` 351:              // *  1, and scale to 0, and a non-trivial solution to A*x = 0 is found.`
` 352:              // *`
` 353:              // *  Similarly, a row-wise scheme is used to solve A'*x = b.  The basic`
` 354:              // *  algorithm for A upper triangular is`
` 355:              // *`
` 356:              // *       for j = 1, ..., n`
` 357:              // *            x(j) := ( b(j) - A[1:j-1,j]' * x[1:j-1] ) / A(j,j)`
` 358:              // *       end`
` 359:              // *`
` 360:              // *  We simultaneously compute two bounds`
` 361:              // *       G(j) = bound on ( b(i) - A[1:i-1,i]' * x[1:i-1] ), 1<=i<=j`
` 362:              // *       M(j) = bound on x(i), 1<=i<=j`
` 363:              // *`
` 364:              // *  The initial values are G(0) = 0, M(0) = max{b(i), i=1,..,n}, and we`
` 365:              // *  add the constraint G(j) >= G(j-1) and M(j) >= M(j-1) for j >= 1.`
` 366:              // *  Then the bound on x(j) is`
` 367:              // *`
` 368:              // *       M(j) <= M(j-1) * ( 1 + CNORM(j) ) / | A(j,j) |`
` 369:              // *`
` 370:              // *            <= M(0) * product ( ( 1 + CNORM(i) ) / |A(i,i)| )`
` 371:              // *                      1<=i<=j`
` 372:              // *`
` 373:              // *  and we can safely call DTRSV if 1/M(n) and 1/G(n) are both greater`
` 374:              // *  than max(underflow, 1/overflow).`
` 375:              // *`
` 376:              // *  =====================================================================`
` 377:              // *`
` 378:              // *     .. Parameters ..`
` 379:              // *     ..`
` 380:              // *     .. Local Scalars ..`
` 381:              // *     ..`
` 382:              // *     .. External Functions ..`
` 383:              // *     ..`
` 384:              // *     .. External Subroutines ..`
` 385:              // *     ..`
` 386:              // *     .. Intrinsic Functions ..`
` 387:              //      INTRINSIC          ABS, MAX, MIN;`
` 388:              // *     ..`
` 389:              // *     .. Executable Statements ..`
` 390:              // *`
` 391:   `
` 392:              #endregion`
` 393:   `
` 394:   `
` 395:              #region Body`
` 396:              `
` 397:              INFO = 0;`
` 398:              UPPER = this._lsame.Run(UPLO, "U");`
` 399:              NOTRAN = this._lsame.Run(TRANS, "N");`
` 400:              NOUNIT = this._lsame.Run(DIAG, "N");`
` 401:              // *`
` 402:              // *     Test the input parameters.`
` 403:              // *`
` 404:              if (!UPPER && !this._lsame.Run(UPLO, "L"))`
` 405:              {`
` 406:                  INFO =  - 1;`
` 407:              }`
` 408:              else`
` 409:              {`
` 410:                  if (!NOTRAN && !this._lsame.Run(TRANS, "T") && !this._lsame.Run(TRANS, "C"))`
` 411:                  {`
` 412:                      INFO =  - 2;`
` 413:                  }`
` 414:                  else`
` 415:                  {`
` 416:                      if (!NOUNIT && !this._lsame.Run(DIAG, "U"))`
` 417:                      {`
` 418:                          INFO =  - 3;`
` 419:                      }`
` 420:                      else`
` 421:                      {`
` 422:                          if (!this._lsame.Run(NORMIN, "Y") && !this._lsame.Run(NORMIN, "N"))`
` 423:                          {`
` 424:                              INFO =  - 4;`
` 425:                          }`
` 426:                          else`
` 427:                          {`
` 428:                              if (N < 0)`
` 429:                              {`
` 430:                                  INFO =  - 5;`
` 431:                              }`
` 432:                              else`
` 433:                              {`
` 434:                                  if (LDA < Math.Max(1, N))`
` 435:                                  {`
` 436:                                      INFO =  - 7;`
` 437:                                  }`
` 438:                              }`
` 439:                          }`
` 440:                      }`
` 441:                  }`
` 442:              }`
` 443:              if (INFO != 0)`
` 444:              {`
` 445:                  this._xerbla.Run("DLATRS",  - INFO);`
` 446:                  return;`
` 447:              }`
` 448:              // *`
` 449:              // *     Quick return if possible`
` 450:              // *`
` 451:              if (N == 0) return;`
` 452:              // *`
` 453:              // *     Determine machine dependent parameters to control overflow.`
` 454:              // *`
` 455:              SMLNUM = this._dlamch.Run("Safe minimum") / this._dlamch.Run("Precision");`
` 456:              BIGNUM = ONE / SMLNUM;`
` 457:              SCALE = ONE;`
` 458:              // *`
` 459:              if (this._lsame.Run(NORMIN, "N"))`
` 460:              {`
` 461:                  // *`
` 462:                  // *        Compute the 1-norm of each column, not including the diagonal.`
` 463:                  // *`
` 464:                  if (UPPER)`
` 465:                  {`
` 466:                      // *`
` 467:                      // *           A is upper triangular.`
` 468:                      // *`
` 469:                      for (J = 1; J <= N; J++)`
` 470:                      {`
` 471:                          CNORM[J + o_cnorm] = this._dasum.Run(J - 1, A, 1+J * LDA + o_a, 1);`
` 472:                      }`
` 473:                  }`
` 474:                  else`
` 475:                  {`
` 476:                      // *`
` 477:                      // *           A is lower triangular.`
` 478:                      // *`
` 479:                      for (J = 1; J <= N - 1; J++)`
` 480:                      {`
` 481:                          CNORM[J + o_cnorm] = this._dasum.Run(N - J, A, J + 1+J * LDA + o_a, 1);`
` 482:                      }`
` 483:                      CNORM[N + o_cnorm] = ZERO;`
` 484:                  }`
` 485:              }`
` 486:              // *`
` 487:              // *     Scale the column norms by TSCAL if the maximum element in CNORM is`
` 488:              // *     greater than BIGNUM.`
` 489:              // *`
` 490:              IMAX = this._idamax.Run(N, CNORM, offset_cnorm, 1);`
` 491:              TMAX = CNORM[IMAX + o_cnorm];`
` 492:              if (TMAX <= BIGNUM)`
` 493:              {`
` 494:                  TSCAL = ONE;`
` 495:              }`
` 496:              else`
` 497:              {`
` 498:                  TSCAL = ONE / (SMLNUM * TMAX);`
` 499:                  this._dscal.Run(N, TSCAL, ref CNORM, offset_cnorm, 1);`
` 500:              }`
` 501:              // *`
` 502:              // *     Compute a bound on the computed solution vector to see if the`
` 503:              // *     Level 2 BLAS routine DTRSV can be used.`
` 504:              // *`
` 505:              J = this._idamax.Run(N, X, offset_x, 1);`
` 506:              XMAX = Math.Abs(X[J + o_x]);`
` 507:              XBND = XMAX;`
` 508:              if (NOTRAN)`
` 509:              {`
` 510:                  // *`
` 511:                  // *        Compute the growth in A * x = b.`
` 512:                  // *`
` 513:                  if (UPPER)`
` 514:                  {`
` 515:                      JFIRST = N;`
` 516:                      JLAST = 1;`
` 517:                      JINC =  - 1;`
` 518:                  }`
` 519:                  else`
` 520:                  {`
` 521:                      JFIRST = 1;`
` 522:                      JLAST = N;`
` 523:                      JINC = 1;`
` 524:                  }`
` 525:                  // *`
` 526:                  if (TSCAL != ONE)`
` 527:                  {`
` 528:                      GROW = ZERO;`
` 529:                      goto LABEL50;`
` 530:                  }`
` 531:                  // *`
` 532:                  if (NOUNIT)`
` 533:                  {`
` 534:                      // *`
` 535:                      // *           A is non-unit triangular.`
` 536:                      // *`
` 537:                      // *           Compute GROW = 1/G(j) and XBND = 1/M(j).`
` 538:                      // *           Initially, G(0) = max{x(i), i=1,...,n}.`
` 539:                      // *`
` 540:                      GROW = ONE / Math.Max(XBND, SMLNUM);`
` 541:                      XBND = GROW;`
` 542:                      for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)`
` 543:                      {`
` 544:                          // *`
` 545:                          // *              Exit the loop if the growth factor is too small.`
` 546:                          // *`
` 547:                          if (GROW <= SMLNUM) goto LABEL50;`
` 548:                          // *`
` 549:                          // *              M(j) = G(j-1) / abs(A(j,j))`
` 550:                          // *`
` 551:                          TJJ = Math.Abs(A[J+J * LDA + o_a]);`
` 552:                          XBND = Math.Min(XBND, Math.Min(ONE, TJJ) * GROW);`
` 553:                          if (TJJ + CNORM[J + o_cnorm] >= SMLNUM)`
` 554:                          {`
` 555:                              // *`
` 556:                              // *                 G(j) = G(j-1)*( 1 + CNORM(j) / abs(A(j,j)) )`
` 557:                              // *`
` 558:                              GROW = GROW * (TJJ / (TJJ + CNORM[J + o_cnorm]));`
` 559:                          }`
` 560:                          else`
` 561:                          {`
` 562:                              // *`
` 563:                              // *                 G(j) could overflow, set GROW to 0.`
` 564:                              // *`
` 565:                              GROW = ZERO;`
` 566:                          }`
` 567:                      }`
` 568:                      GROW = XBND;`
` 569:                  }`
` 570:                  else`
` 571:                  {`
` 572:                      // *`
` 573:                      // *           A is unit triangular.`
` 574:                      // *`
` 575:                      // *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.`
` 576:                      // *`
` 577:                      GROW = Math.Min(ONE, ONE / Math.Max(XBND, SMLNUM));`
` 578:                      for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)`
` 579:                      {`
` 580:                          // *`
` 581:                          // *              Exit the loop if the growth factor is too small.`
` 582:                          // *`
` 583:                          if (GROW <= SMLNUM) goto LABEL50;`
` 584:                          // *`
` 585:                          // *              G(j) = G(j-1)*( 1 + CNORM(j) )`
` 586:                          // *`
` 587:                          GROW = GROW * (ONE / (ONE + CNORM[J + o_cnorm]));`
` 588:                      }`
` 589:                  }`
` 590:              LABEL50:;`
` 591:                  // *`
` 592:              }`
` 593:              else`
` 594:              {`
` 595:                  // *`
` 596:                  // *        Compute the growth in A' * x = b.`
` 597:                  // *`
` 598:                  if (UPPER)`
` 599:                  {`
` 600:                      JFIRST = 1;`
` 601:                      JLAST = N;`
` 602:                      JINC = 1;`
` 603:                  }`
` 604:                  else`
` 605:                  {`
` 606:                      JFIRST = N;`
` 607:                      JLAST = 1;`
` 608:                      JINC =  - 1;`
` 609:                  }`
` 610:                  // *`
` 611:                  if (TSCAL != ONE)`
` 612:                  {`
` 613:                      GROW = ZERO;`
` 614:                      goto LABEL80;`
` 615:                  }`
` 616:                  // *`
` 617:                  if (NOUNIT)`
` 618:                  {`
` 619:                      // *`
` 620:                      // *           A is non-unit triangular.`
` 621:                      // *`
` 622:                      // *           Compute GROW = 1/G(j) and XBND = 1/M(j).`
` 623:                      // *           Initially, M(0) = max{x(i), i=1,...,n}.`
` 624:                      // *`
` 625:                      GROW = ONE / Math.Max(XBND, SMLNUM);`
` 626:                      XBND = GROW;`
` 627:                      for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)`
` 628:                      {`
` 629:                          // *`
` 630:                          // *              Exit the loop if the growth factor is too small.`
` 631:                          // *`
` 632:                          if (GROW <= SMLNUM) goto LABEL80;`
` 633:                          // *`
` 634:                          // *              G(j) = max( G(j-1), M(j-1)*( 1 + CNORM(j) ) )`
` 635:                          // *`
` 636:                          XJ = ONE + CNORM[J + o_cnorm];`
` 637:                          GROW = Math.Min(GROW, XBND / XJ);`
` 638:                          // *`
` 639:                          // *              M(j) = M(j-1)*( 1 + CNORM(j) ) / abs(A(j,j))`
` 640:                          // *`
` 641:                          TJJ = Math.Abs(A[J+J * LDA + o_a]);`
` 642:                          if (XJ > TJJ) XBND = XBND * (TJJ / XJ);`
` 643:                      }`
` 644:                      GROW = Math.Min(GROW, XBND);`
` 645:                  }`
` 646:                  else`
` 647:                  {`
` 648:                      // *`
` 649:                      // *           A is unit triangular.`
` 650:                      // *`
` 651:                      // *           Compute GROW = 1/G(j), where G(0) = max{x(i), i=1,...,n}.`
` 652:                      // *`
` 653:                      GROW = Math.Min(ONE, ONE / Math.Max(XBND, SMLNUM));`
` 654:                      for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)`
` 655:                      {`
` 656:                          // *`
` 657:                          // *              Exit the loop if the growth factor is too small.`
` 658:                          // *`
` 659:                          if (GROW <= SMLNUM) goto LABEL80;`
` 660:                          // *`
` 661:                          // *              G(j) = ( 1 + CNORM(j) )*G(j-1)`
` 662:                          // *`
` 663:                          XJ = ONE + CNORM[J + o_cnorm];`
` 664:                          GROW = GROW / XJ;`
` 665:                      }`
` 666:                  }`
` 667:              LABEL80:;`
` 668:              }`
` 669:              // *`
` 670:              if ((GROW * TSCAL) > SMLNUM)`
` 671:              {`
` 672:                  // *`
` 673:                  // *        Use the Level 2 BLAS solve if the reciprocal of the bound on`
` 674:                  // *        elements of X is not too small.`
` 675:                  // *`
` 676:                  this._dtrsv.Run(UPLO, TRANS, DIAG, N, A, offset_a, LDA`
` 677:                                  , ref X, offset_x, 1);`
` 678:              }`
` 679:              else`
` 680:              {`
` 681:                  // *`
` 682:                  // *        Use a Level 1 BLAS solve, scaling intermediate results.`
` 683:                  // *`
` 684:                  if (XMAX > BIGNUM)`
` 685:                  {`
` 686:                      // *`
` 687:                      // *           Scale X so that its components are less than or equal to`
` 688:                      // *           BIGNUM in absolute value.`
` 689:                      // *`
` 690:                      SCALE = BIGNUM / XMAX;`
` 691:                      this._dscal.Run(N, SCALE, ref X, offset_x, 1);`
` 692:                      XMAX = BIGNUM;`
` 693:                  }`
` 694:                  // *`
` 695:                  if (NOTRAN)`
` 696:                  {`
` 697:                      // *`
` 698:                      // *           Solve A * x = b`
` 699:                      // *`
` 700:                      for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)`
` 701:                      {`
` 702:                          // *`
` 703:                          // *              Compute x(j) = b(j) / A(j,j), scaling x if necessary.`
` 704:                          // *`
` 705:                          XJ = Math.Abs(X[J + o_x]);`
` 706:                          if (NOUNIT)`
` 707:                          {`
` 708:                              TJJS = A[J+J * LDA + o_a] * TSCAL;`
` 709:                          }`
` 710:                          else`
` 711:                          {`
` 712:                              TJJS = TSCAL;`
` 713:                              if (TSCAL == ONE) goto LABEL100;`
` 714:                          }`
` 715:                          TJJ = Math.Abs(TJJS);`
` 716:                          if (TJJ > SMLNUM)`
` 717:                          {`
` 718:                              // *`
` 719:                              // *                    abs(A(j,j)) > SMLNUM:`
` 720:                              // *`
` 721:                              if (TJJ < ONE)`
` 722:                              {`
` 723:                                  if (XJ > TJJ * BIGNUM)`
` 724:                                  {`
` 725:                                      // *`
` 726:                                      // *                          Scale x by 1/b(j).`
` 727:                                      // *`
` 728:                                      REC = ONE / XJ;`
` 729:                                      this._dscal.Run(N, REC, ref X, offset_x, 1);`
` 730:                                      SCALE = SCALE * REC;`
` 731:                                      XMAX = XMAX * REC;`
` 732:                                  }`
` 733:                              }`
` 734:                              X[J + o_x] = X[J + o_x] / TJJS;`
` 735:                              XJ = Math.Abs(X[J + o_x]);`
` 736:                          }`
` 737:                          else`
` 738:                          {`
` 739:                              if (TJJ > ZERO)`
` 740:                              {`
` 741:                                  // *`
` 742:                                  // *                    0 < abs(A(j,j)) <= SMLNUM:`
` 743:                                  // *`
` 744:                                  if (XJ > TJJ * BIGNUM)`
` 745:                                  {`
` 746:                                      // *`
` 747:                                      // *                       Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM`
` 748:                                      // *                       to avoid overflow when dividing by A(j,j).`
` 749:                                      // *`
` 750:                                      REC = (TJJ * BIGNUM) / XJ;`
` 751:                                      if (CNORM[J + o_cnorm] > ONE)`
` 752:                                      {`
` 753:                                          // *`
` 754:                                          // *                          Scale by 1/CNORM(j) to avoid overflow when`
` 755:                                          // *                          multiplying x(j) times column j.`
` 756:                                          // *`
` 757:                                          REC = REC / CNORM[J + o_cnorm];`
` 758:                                      }`
` 759:                                      this._dscal.Run(N, REC, ref X, offset_x, 1);`
` 760:                                      SCALE = SCALE * REC;`
` 761:                                      XMAX = XMAX * REC;`
` 762:                                  }`
` 763:                                  X[J + o_x] = X[J + o_x] / TJJS;`
` 764:                                  XJ = Math.Abs(X[J + o_x]);`
` 765:                              }`
` 766:                              else`
` 767:                              {`
` 768:                                  // *`
` 769:                                  // *                    A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and`
` 770:                                  // *                    scale = 0, and compute a solution to A*x = 0.`
` 771:                                  // *`
` 772:                                  for (I = 1; I <= N; I++)`
` 773:                                  {`
` 774:                                      X[I + o_x] = ZERO;`
` 775:                                  }`
` 776:                                  X[J + o_x] = ONE;`
` 777:                                  XJ = ONE;`
` 778:                                  SCALE = ZERO;`
` 779:                                  XMAX = ZERO;`
` 780:                              }`
` 781:                          }`
` 782:                      LABEL100:;`
` 783:                          // *`
` 784:                          // *              Scale x if necessary to avoid overflow when adding a`
` 785:                          // *              multiple of column j of A.`
` 786:                          // *`
` 787:                          if (XJ > ONE)`
` 788:                          {`
` 789:                              REC = ONE / XJ;`
` 790:                              if (CNORM[J + o_cnorm] > (BIGNUM - XMAX) * REC)`
` 791:                              {`
` 792:                                  // *`
` 793:                                  // *                    Scale x by 1/(2*abs(x(j))).`
` 794:                                  // *`
` 795:                                  REC = REC * HALF;`
` 796:                                  this._dscal.Run(N, REC, ref X, offset_x, 1);`
` 797:                                  SCALE = SCALE * REC;`
` 798:                              }`
` 799:                          }`
` 800:                          else`
` 801:                          {`
` 802:                              if (XJ * CNORM[J + o_cnorm] > (BIGNUM - XMAX))`
` 803:                              {`
` 804:                                  // *`
` 805:                                  // *                 Scale x by 1/2.`
` 806:                                  // *`
` 807:                                  this._dscal.Run(N, HALF, ref X, offset_x, 1);`
` 808:                                  SCALE = SCALE * HALF;`
` 809:                              }`
` 810:                          }`
` 811:                          // *`
` 812:                          if (UPPER)`
` 813:                          {`
` 814:                              if (J > 1)`
` 815:                              {`
` 816:                                  // *`
` 817:                                  // *                    Compute the update`
` 818:                                  // *                       x(1:j-1) := x(1:j-1) - x(j) * A(1:j-1,j)`
` 819:                                  // *`
` 820:                                  this._daxpy.Run(J - 1,  - X[J + o_x] * TSCAL, A, 1+J * LDA + o_a, 1, ref X, offset_x, 1);`
` 821:                                  I = this._idamax.Run(J - 1, X, offset_x, 1);`
` 822:                                  XMAX = Math.Abs(X[I + o_x]);`
` 823:                              }`
` 824:                          }`
` 825:                          else`
` 826:                          {`
` 827:                              if (J < N)`
` 828:                              {`
` 829:                                  // *`
` 830:                                  // *                    Compute the update`
` 831:                                  // *                       x(j+1:n) := x(j+1:n) - x(j) * A(j+1:n,j)`
` 832:                                  // *`
` 833:                                  this._daxpy.Run(N - J,  - X[J + o_x] * TSCAL, A, J + 1+J * LDA + o_a, 1, ref X, J + 1 + o_x, 1);`
` 834:                                  I = J + this._idamax.Run(N - J, X, J + 1 + o_x, 1);`
` 835:                                  XMAX = Math.Abs(X[I + o_x]);`
` 836:                              }`
` 837:                          }`
` 838:                      }`
` 839:                      // *`
` 840:                  }`
` 841:                  else`
` 842:                  {`
` 843:                      // *`
` 844:                      // *           Solve A' * x = b`
` 845:                      // *`
` 846:                      for (J = JFIRST; (JINC >= 0) ? (J <= JLAST) : (J >= JLAST); J += JINC)`
` 847:                      {`
` 848:                          // *`
` 849:                          // *              Compute x(j) = b(j) - sum A(k,j)*x(k).`
` 850:                          // *                                    k<>j`
` 851:                          // *`
` 852:                          XJ = Math.Abs(X[J + o_x]);`
` 853:                          USCAL = TSCAL;`
` 854:                          REC = ONE / Math.Max(XMAX, ONE);`
` 855:                          if (CNORM[J + o_cnorm] > (BIGNUM - XJ) * REC)`
` 856:                          {`
` 857:                              // *`
` 858:                              // *                 If x(j) could overflow, scale x by 1/(2*XMAX).`
` 859:                              // *`
` 860:                              REC = REC * HALF;`
` 861:                              if (NOUNIT)`
` 862:                              {`
` 863:                                  TJJS = A[J+J * LDA + o_a] * TSCAL;`
` 864:                              }`
` 865:                              else`
` 866:                              {`
` 867:                                  TJJS = TSCAL;`
` 868:                              }`
` 869:                              TJJ = Math.Abs(TJJS);`
` 870:                              if (TJJ > ONE)`
` 871:                              {`
` 872:                                  // *`
` 873:                                  // *                       Divide by A(j,j) when scaling x if A(j,j) > 1.`
` 874:                                  // *`
` 875:                                  REC = Math.Min(ONE, REC * TJJ);`
` 876:                                  USCAL = USCAL / TJJS;`
` 877:                              }`
` 878:                              if (REC < ONE)`
` 879:                              {`
` 880:                                  this._dscal.Run(N, REC, ref X, offset_x, 1);`
` 881:                                  SCALE = SCALE * REC;`
` 882:                                  XMAX = XMAX * REC;`
` 883:                              }`
` 884:                          }`
` 885:                          // *`
` 886:                          SUMJ = ZERO;`
` 887:                          if (USCAL == ONE)`
` 888:                          {`
` 889:                              // *`
` 890:                              // *                 If the scaling needed for A in the dot product is 1,`
` 891:                              // *                 call DDOT to perform the dot product.`
` 892:                              // *`
` 893:                              if (UPPER)`
` 894:                              {`
` 895:                                  SUMJ = this._ddot.Run(J - 1, A, 1+J * LDA + o_a, 1, X, offset_x, 1);`
` 896:                              }`
` 897:                              else`
` 898:                              {`
` 899:                                  if (J < N)`
` 900:                                  {`
` 901:                                      SUMJ = this._ddot.Run(N - J, A, J + 1+J * LDA + o_a, 1, X, J + 1 + o_x, 1);`
` 902:                                  }`
` 903:                              }`
` 904:                          }`
` 905:                          else`
` 906:                          {`
` 907:                              // *`
` 908:                              // *                 Otherwise, use in-line code for the dot product.`
` 909:                              // *`
` 910:                              if (UPPER)`
` 911:                              {`
` 912:                                  for (I = 1; I <= J - 1; I++)`
` 913:                                  {`
` 914:                                      SUMJ = SUMJ + (A[I+J * LDA + o_a] * USCAL) * X[I + o_x];`
` 915:                                  }`
` 916:                              }`
` 917:                              else`
` 918:                              {`
` 919:                                  if (J < N)`
` 920:                                  {`
` 921:                                      for (I = J + 1; I <= N; I++)`
` 922:                                      {`
` 923:                                          SUMJ = SUMJ + (A[I+J * LDA + o_a] * USCAL) * X[I + o_x];`
` 924:                                      }`
` 925:                                  }`
` 926:                              }`
` 927:                          }`
` 928:                          // *`
` 929:                          if (USCAL == TSCAL)`
` 930:                          {`
` 931:                              // *`
` 932:                              // *                 Compute x(j) := ( x(j) - sumj ) / A(j,j) if 1/A(j,j)`
` 933:                              // *                 was not used to scale the dotproduct.`
` 934:                              // *`
` 935:                              X[J + o_x] = X[J + o_x] - SUMJ;`
` 936:                              XJ = Math.Abs(X[J + o_x]);`
` 937:                              if (NOUNIT)`
` 938:                              {`
` 939:                                  TJJS = A[J+J * LDA + o_a] * TSCAL;`
` 940:                              }`
` 941:                              else`
` 942:                              {`
` 943:                                  TJJS = TSCAL;`
` 944:                                  if (TSCAL == ONE) goto LABEL150;`
` 945:                              }`
` 946:                              // *`
` 947:                              // *                    Compute x(j) = x(j) / A(j,j), scaling if necessary.`
` 948:                              // *`
` 949:                              TJJ = Math.Abs(TJJS);`
` 950:                              if (TJJ > SMLNUM)`
` 951:                              {`
` 952:                                  // *`
` 953:                                  // *                       abs(A(j,j)) > SMLNUM:`
` 954:                                  // *`
` 955:                                  if (TJJ < ONE)`
` 956:                                  {`
` 957:                                      if (XJ > TJJ * BIGNUM)`
` 958:                                      {`
` 959:                                          // *`
` 960:                                          // *                             Scale X by 1/abs(x(j)).`
` 961:                                          // *`
` 962:                                          REC = ONE / XJ;`
` 963:                                          this._dscal.Run(N, REC, ref X, offset_x, 1);`
` 964:                                          SCALE = SCALE * REC;`
` 965:                                          XMAX = XMAX * REC;`
` 966:                                      }`
` 967:                                  }`
` 968:                                  X[J + o_x] = X[J + o_x] / TJJS;`
` 969:                              }`
` 970:                              else`
` 971:                              {`
` 972:                                  if (TJJ > ZERO)`
` 973:                                  {`
` 974:                                      // *`
` 975:                                      // *                       0 < abs(A(j,j)) <= SMLNUM:`
` 976:                                      // *`
` 977:                                      if (XJ > TJJ * BIGNUM)`
` 978:                                      {`
` 979:                                          // *`
` 980:                                          // *                          Scale x by (1/abs(x(j)))*abs(A(j,j))*BIGNUM.`
` 981:                                          // *`
` 982:                                          REC = (TJJ * BIGNUM) / XJ;`
` 983:                                          this._dscal.Run(N, REC, ref X, offset_x, 1);`
` 984:                                          SCALE = SCALE * REC;`
` 985:                                          XMAX = XMAX * REC;`
` 986:                                      }`
` 987:                                      X[J + o_x] = X[J + o_x] / TJJS;`
` 988:                                  }`
` 989:                                  else`
` 990:                                  {`
` 991:                                      // *`
` 992:                                      // *                       A(j,j) = 0:  Set x(1:n) = 0, x(j) = 1, and`
` 993:                                      // *                       scale = 0, and compute a solution to A'*x = 0.`
` 994:                                      // *`
` 995:                                      for (I = 1; I <= N; I++)`
` 996:                                      {`
` 997:                                          X[I + o_x] = ZERO;`
` 998:                                      }`
` 999:                                      X[J + o_x] = ONE;`
`1000:                                      SCALE = ZERO;`
`1001:                                      XMAX = ZERO;`
`1002:                                  }`
`1003:                              }`
`1004:                          LABEL150:;`
`1005:                          }`
`1006:                          else`
`1007:                          {`
`1008:                              // *`
`1009:                              // *                 Compute x(j) := x(j) / A(j,j)  - sumj if the dot`
`1010:                              // *                 product has already been divided by 1/A(j,j).`
`1011:                              // *`
`1012:                              X[J + o_x] = X[J + o_x] / TJJS - SUMJ;`
`1013:                          }`
`1014:                          XMAX = Math.Max(XMAX, Math.Abs(X[J + o_x]));`
`1015:                      }`
`1016:                  }`
`1017:                  SCALE = SCALE / TSCAL;`
`1018:              }`
`1019:              // *`
`1020:              // *     Scale the column norms by 1/TSCAL for return.`
`1021:              // *`
`1022:              if (TSCAL != ONE)`
`1023:              {`
`1024:                  this._dscal.Run(N, ONE / TSCAL, ref CNORM, offset_cnorm, 1);`
`1025:              }`
`1026:              // *`
`1027:              return;`
`1028:              // *`
`1029:              // *     End of DLATRS`
`1030:              // *`
`1031:   `
`1032:              #endregion`
`1033:   `
`1034:          }`
`1035:      }`
`1036:  }`