  `   1:  #region Translated by Jose Antonio De Santiago-Castillo.`
`   2:   `
`   3:  //Translated by Jose Antonio De Santiago-Castillo. `
`   4:  //E-mail:JAntonioDeSantiago@gmail.com`
`   5:  //Web: www.DotNumerics.com`
`   6:  //`
`   7:  //Fortran to C# Translation.`
`   8:  //Translated by:`
`   9:  //F2CSharp Version 0.71 (November 10, 2009)`
`  10:  //Code Optimizations: None`
`  11:  //`
`  12:  #endregion`
`  13:   `
`  14:  using System;`
`  15:  using DotNumerics.FortranLibrary;`
`  16:   `
`  17:  namespace DotNumerics.CSLapack`
`  18:  {`
`  19:      /// <summary>`
`  20:      /// -- LAPACK routine (version 3.1) --`
`  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
`  22:      /// November 2006`
`  23:      /// Purpose`
`  24:      /// =======`
`  25:      /// `
`  26:      /// DGTSV  solves the equation`
`  27:      /// `
`  28:      /// A*X = B,`
`  29:      /// `
`  30:      /// where A is an n by n tridiagonal matrix, by Gaussian elimination with`
`  31:      /// partial pivoting.`
`  32:      /// `
`  33:      /// Note that the equation  A'*X = B  may be solved by interchanging the`
`  34:      /// order of the arguments DU and DL.`
`  35:      /// `
`  36:      ///</summary>`
`  37:      public class DGTSV`
`  38:      {`
`  39:      `
`  40:   `
`  41:          #region Dependencies`
`  42:          `
`  43:          XERBLA _xerbla; `
`  44:   `
`  45:          #endregion`
`  46:   `
`  47:   `
`  48:          #region Fields`
`  49:          `
`  50:          const double ZERO = 0.0E+0; int I = 0; int J = 0; double FACT = 0; double TEMP = 0; `
`  51:   `
`  52:          #endregion`
`  53:   `
`  54:          public DGTSV(XERBLA xerbla)`
`  55:          {`
`  56:      `
`  57:   `
`  58:              #region Set Dependencies`
`  59:              `
`  60:              this._xerbla = xerbla; `
`  61:   `
`  62:              #endregion`
`  63:   `
`  64:          }`
`  65:      `
`  66:          public DGTSV()`
`  67:          {`
`  68:      `
`  69:   `
`  70:              #region Dependencies (Initialization)`
`  71:              `
`  72:              XERBLA xerbla = new XERBLA();`
`  73:   `
`  74:              #endregion`
`  75:   `
`  76:   `
`  77:              #region Set Dependencies`
`  78:              `
`  79:              this._xerbla = xerbla; `
`  80:   `
`  81:              #endregion`
`  82:   `
`  83:          }`
`  84:          /// <summary>`
`  85:          /// Purpose`
`  86:          /// =======`
`  87:          /// `
`  88:          /// DGTSV  solves the equation`
`  89:          /// `
`  90:          /// A*X = B,`
`  91:          /// `
`  92:          /// where A is an n by n tridiagonal matrix, by Gaussian elimination with`
`  93:          /// partial pivoting.`
`  94:          /// `
`  95:          /// Note that the equation  A'*X = B  may be solved by interchanging the`
`  96:          /// order of the arguments DU and DL.`
`  97:          /// `
`  98:          ///</summary>`
`  99:          /// <param name="N">`
` 100:          /// (input) INTEGER`
` 101:          /// The order of the matrix A.  N .GE. 0.`
` 102:          ///</param>`
` 103:          /// <param name="NRHS">`
` 104:          /// (input) INTEGER`
` 105:          /// The number of right hand sides, i.e., the number of columns`
` 106:          /// of the matrix B.  NRHS .GE. 0.`
` 107:          ///</param>`
` 108:          /// <param name="DL">`
` 109:          /// (input/output) DOUBLE PRECISION array, dimension (N-1)`
` 110:          /// On entry, DL must contain the (n-1) sub-diagonal elements of`
` 111:          /// A.`
` 112:          /// `
` 113:          /// On exit, DL is overwritten by the (n-2) elements of the`
` 114:          /// second super-diagonal of the upper triangular matrix U from`
` 115:          /// the LU factorization of A, in DL(1), ..., DL(n-2).`
` 116:          ///</param>`
` 117:          /// <param name="D">`
` 118:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 119:          /// On entry, D must contain the diagonal elements of A.`
` 120:          /// `
` 121:          /// On exit, D is overwritten by the n diagonal elements of U.`
` 122:          ///</param>`
` 123:          /// <param name="DU">`
` 124:          /// (input/output) DOUBLE PRECISION array, dimension (N-1)`
` 125:          /// On entry, DU must contain the (n-1) super-diagonal elements`
` 126:          /// of A.`
` 127:          /// `
` 128:          /// On exit, DU is overwritten by the (n-1) elements of the first`
` 129:          /// super-diagonal of U.`
` 130:          ///</param>`
` 131:          /// <param name="B">`
` 132:          /// (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)`
` 133:          /// On entry, the N by NRHS matrix of right hand side matrix B.`
` 134:          /// On exit, if INFO = 0, the N by NRHS solution matrix X.`
` 135:          ///</param>`
` 136:          /// <param name="LDB">`
` 137:          /// (input) INTEGER`
` 138:          /// The leading dimension of the array B.  LDB .GE. max(1,N).`
` 139:          ///</param>`
` 140:          /// <param name="INFO">`
` 141:          /// (output) INTEGER`
` 142:          /// = 0: successful exit`
` 143:          /// .LT. 0: if INFO = -i, the i-th argument had an illegal value`
` 144:          /// .GT. 0: if INFO = i, U(i,i) is exactly zero, and the solution`
` 145:          /// has not been computed.  The factorization has not been`
` 146:          /// completed unless i = N.`
` 147:          ///</param>`
` 148:          public void Run(int N, int NRHS, ref double[] DL, int offset_dl, ref double[] D, int offset_d, ref double[] DU, int offset_du, ref double[] B, int offset_b`
` 149:                           , int LDB, ref int INFO)`
` 150:          {`
` 151:   `
` 152:              #region Array Index Correction`
` 153:              `
` 154:               int o_dl = -1 + offset_dl;  int o_d = -1 + offset_d;  int o_du = -1 + offset_du;  int o_b = -1 - LDB + offset_b; `
` 155:   `
` 156:              #endregion`
` 157:   `
` 158:   `
` 159:              #region Prolog`
` 160:              `
` 161:              // *`
` 162:              // *  -- LAPACK routine (version 3.1) --`
` 163:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 164:              // *     November 2006`
` 165:              // *`
` 166:              // *     .. Scalar Arguments ..`
` 167:              // *     ..`
` 168:              // *     .. Array Arguments ..`
` 169:              // *     ..`
` 170:              // *`
` 171:              // *  Purpose`
` 172:              // *  =======`
` 173:              // *`
` 174:              // *  DGTSV  solves the equation`
` 175:              // *`
` 176:              // *     A*X = B,`
` 177:              // *`
` 178:              // *  where A is an n by n tridiagonal matrix, by Gaussian elimination with`
` 179:              // *  partial pivoting.`
` 180:              // *`
` 181:              // *  Note that the equation  A'*X = B  may be solved by interchanging the`
` 182:              // *  order of the arguments DU and DL.`
` 183:              // *`
` 184:              // *  Arguments`
` 185:              // *  =========`
` 186:              // *`
` 187:              // *  N       (input) INTEGER`
` 188:              // *          The order of the matrix A.  N >= 0.`
` 189:              // *`
` 190:              // *  NRHS    (input) INTEGER`
` 191:              // *          The number of right hand sides, i.e., the number of columns`
` 192:              // *          of the matrix B.  NRHS >= 0.`
` 193:              // *`
` 194:              // *  DL      (input/output) DOUBLE PRECISION array, dimension (N-1)`
` 195:              // *          On entry, DL must contain the (n-1) sub-diagonal elements of`
` 196:              // *          A.`
` 197:              // *`
` 198:              // *          On exit, DL is overwritten by the (n-2) elements of the`
` 199:              // *          second super-diagonal of the upper triangular matrix U from`
` 200:              // *          the LU factorization of A, in DL(1), ..., DL(n-2).`
` 201:              // *`
` 202:              // *  D       (input/output) DOUBLE PRECISION array, dimension (N)`
` 203:              // *          On entry, D must contain the diagonal elements of A.`
` 204:              // *`
` 205:              // *          On exit, D is overwritten by the n diagonal elements of U.`
` 206:              // *`
` 207:              // *  DU      (input/output) DOUBLE PRECISION array, dimension (N-1)`
` 208:              // *          On entry, DU must contain the (n-1) super-diagonal elements`
` 209:              // *          of A.`
` 210:              // *`
` 211:              // *          On exit, DU is overwritten by the (n-1) elements of the first`
` 212:              // *          super-diagonal of U.`
` 213:              // *`
` 214:              // *  B       (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)`
` 215:              // *          On entry, the N by NRHS matrix of right hand side matrix B.`
` 216:              // *          On exit, if INFO = 0, the N by NRHS solution matrix X.`
` 217:              // *`
` 218:              // *  LDB     (input) INTEGER`
` 219:              // *          The leading dimension of the array B.  LDB >= max(1,N).`
` 220:              // *`
` 221:              // *  INFO    (output) INTEGER`
` 222:              // *          = 0: successful exit`
` 223:              // *          < 0: if INFO = -i, the i-th argument had an illegal value`
` 224:              // *          > 0: if INFO = i, U(i,i) is exactly zero, and the solution`
` 225:              // *               has not been computed.  The factorization has not been`
` 226:              // *               completed unless i = N.`
` 227:              // *`
` 228:              // *  =====================================================================`
` 229:              // *`
` 230:              // *     .. Parameters ..`
` 231:              // *     ..`
` 232:              // *     .. Local Scalars ..`
` 233:              // *     ..`
` 234:              // *     .. Intrinsic Functions ..`
` 235:              //      INTRINSIC          ABS, MAX;`
` 236:              // *     ..`
` 237:              // *     .. External Subroutines ..`
` 238:              // *     ..`
` 239:              // *     .. Executable Statements ..`
` 240:              // *`
` 241:   `
` 242:              #endregion`
` 243:   `
` 244:   `
` 245:              #region Body`
` 246:              `
` 247:              INFO = 0;`
` 248:              if (N < 0)`
` 249:              {`
` 250:                  INFO =  - 1;`
` 251:              }`
` 252:              else`
` 253:              {`
` 254:                  if (NRHS < 0)`
` 255:                  {`
` 256:                      INFO =  - 2;`
` 257:                  }`
` 258:                  else`
` 259:                  {`
` 260:                      if (LDB < Math.Max(1, N))`
` 261:                      {`
` 262:                          INFO =  - 7;`
` 263:                      }`
` 264:                  }`
` 265:              }`
` 266:              if (INFO != 0)`
` 267:              {`
` 268:                  this._xerbla.Run("DGTSV ",  - INFO);`
` 269:                  return;`
` 270:              }`
` 271:              // *`
` 272:              if (N == 0) return;`
` 273:              // *`
` 274:              if (NRHS == 1)`
` 275:              {`
` 276:                  for (I = 1; I <= N - 2; I++)`
` 277:                  {`
` 278:                      if (Math.Abs(D[I + o_d]) >= Math.Abs(DL[I + o_dl]))`
` 279:                      {`
` 280:                          // *`
` 281:                          // *              No row interchange required`
` 282:                          // *`
` 283:                          if (D[I + o_d] != ZERO)`
` 284:                          {`
` 285:                              FACT = DL[I + o_dl] / D[I + o_d];`
` 286:                              D[I + 1 + o_d] = D[I + 1 + o_d] - FACT * DU[I + o_du];`
` 287:                              B[I + 1+1 * LDB + o_b] = B[I + 1+1 * LDB + o_b] - FACT * B[I+1 * LDB + o_b];`
` 288:                          }`
` 289:                          else`
` 290:                          {`
` 291:                              INFO = I;`
` 292:                              return;`
` 293:                          }`
` 294:                          DL[I + o_dl] = ZERO;`
` 295:                      }`
` 296:                      else`
` 297:                      {`
` 298:                          // *`
` 299:                          // *              Interchange rows I and I+1`
` 300:                          // *`
` 301:                          FACT = D[I + o_d] / DL[I + o_dl];`
` 302:                          D[I + o_d] = DL[I + o_dl];`
` 303:                          TEMP = D[I + 1 + o_d];`
` 304:                          D[I + 1 + o_d] = DU[I + o_du] - FACT * TEMP;`
` 305:                          DL[I + o_dl] = DU[I + 1 + o_du];`
` 306:                          DU[I + 1 + o_du] =  - FACT * DL[I + o_dl];`
` 307:                          DU[I + o_du] = TEMP;`
` 308:                          TEMP = B[I+1 * LDB + o_b];`
` 309:                          B[I+1 * LDB + o_b] = B[I + 1+1 * LDB + o_b];`
` 310:                          B[I + 1+1 * LDB + o_b] = TEMP - FACT * B[I + 1+1 * LDB + o_b];`
` 311:                      }`
` 312:                  }`
` 313:                  if (N > 1)`
` 314:                  {`
` 315:                      I = N - 1;`
` 316:                      if (Math.Abs(D[I + o_d]) >= Math.Abs(DL[I + o_dl]))`
` 317:                      {`
` 318:                          if (D[I + o_d] != ZERO)`
` 319:                          {`
` 320:                              FACT = DL[I + o_dl] / D[I + o_d];`
` 321:                              D[I + 1 + o_d] = D[I + 1 + o_d] - FACT * DU[I + o_du];`
` 322:                              B[I + 1+1 * LDB + o_b] = B[I + 1+1 * LDB + o_b] - FACT * B[I+1 * LDB + o_b];`
` 323:                          }`
` 324:                          else`
` 325:                          {`
` 326:                              INFO = I;`
` 327:                              return;`
` 328:                          }`
` 329:                      }`
` 330:                      else`
` 331:                      {`
` 332:                          FACT = D[I + o_d] / DL[I + o_dl];`
` 333:                          D[I + o_d] = DL[I + o_dl];`
` 334:                          TEMP = D[I + 1 + o_d];`
` 335:                          D[I + 1 + o_d] = DU[I + o_du] - FACT * TEMP;`
` 336:                          DU[I + o_du] = TEMP;`
` 337:                          TEMP = B[I+1 * LDB + o_b];`
` 338:                          B[I+1 * LDB + o_b] = B[I + 1+1 * LDB + o_b];`
` 339:                          B[I + 1+1 * LDB + o_b] = TEMP - FACT * B[I + 1+1 * LDB + o_b];`
` 340:                      }`
` 341:                  }`
` 342:                  if (D[N + o_d] == ZERO)`
` 343:                  {`
` 344:                      INFO = N;`
` 345:                      return;`
` 346:                  }`
` 347:              }`
` 348:              else`
` 349:              {`
` 350:                  for (I = 1; I <= N - 2; I++)`
` 351:                  {`
` 352:                      if (Math.Abs(D[I + o_d]) >= Math.Abs(DL[I + o_dl]))`
` 353:                      {`
` 354:                          // *`
` 355:                          // *              No row interchange required`
` 356:                          // *`
` 357:                          if (D[I + o_d] != ZERO)`
` 358:                          {`
` 359:                              FACT = DL[I + o_dl] / D[I + o_d];`
` 360:                              D[I + 1 + o_d] = D[I + 1 + o_d] - FACT * DU[I + o_du];`
` 361:                              for (J = 1; J <= NRHS; J++)`
` 362:                              {`
` 363:                                  B[I + 1+J * LDB + o_b] = B[I + 1+J * LDB + o_b] - FACT * B[I+J * LDB + o_b];`
` 364:                              }`
` 365:                          }`
` 366:                          else`
` 367:                          {`
` 368:                              INFO = I;`
` 369:                              return;`
` 370:                          }`
` 371:                          DL[I + o_dl] = ZERO;`
` 372:                      }`
` 373:                      else`
` 374:                      {`
` 375:                          // *`
` 376:                          // *              Interchange rows I and I+1`
` 377:                          // *`
` 378:                          FACT = D[I + o_d] / DL[I + o_dl];`
` 379:                          D[I + o_d] = DL[I + o_dl];`
` 380:                          TEMP = D[I + 1 + o_d];`
` 381:                          D[I + 1 + o_d] = DU[I + o_du] - FACT * TEMP;`
` 382:                          DL[I + o_dl] = DU[I + 1 + o_du];`
` 383:                          DU[I + 1 + o_du] =  - FACT * DL[I + o_dl];`
` 384:                          DU[I + o_du] = TEMP;`
` 385:                          for (J = 1; J <= NRHS; J++)`
` 386:                          {`
` 387:                              TEMP = B[I+J * LDB + o_b];`
` 388:                              B[I+J * LDB + o_b] = B[I + 1+J * LDB + o_b];`
` 389:                              B[I + 1+J * LDB + o_b] = TEMP - FACT * B[I + 1+J * LDB + o_b];`
` 390:                          }`
` 391:                      }`
` 392:                  }`
` 393:                  if (N > 1)`
` 394:                  {`
` 395:                      I = N - 1;`
` 396:                      if (Math.Abs(D[I + o_d]) >= Math.Abs(DL[I + o_dl]))`
` 397:                      {`
` 398:                          if (D[I + o_d] != ZERO)`
` 399:                          {`
` 400:                              FACT = DL[I + o_dl] / D[I + o_d];`
` 401:                              D[I + 1 + o_d] = D[I + 1 + o_d] - FACT * DU[I + o_du];`
` 402:                              for (J = 1; J <= NRHS; J++)`
` 403:                              {`
` 404:                                  B[I + 1+J * LDB + o_b] = B[I + 1+J * LDB + o_b] - FACT * B[I+J * LDB + o_b];`
` 405:                              }`
` 406:                          }`
` 407:                          else`
` 408:                          {`
` 409:                              INFO = I;`
` 410:                              return;`
` 411:                          }`
` 412:                      }`
` 413:                      else`
` 414:                      {`
` 415:                          FACT = D[I + o_d] / DL[I + o_dl];`
` 416:                          D[I + o_d] = DL[I + o_dl];`
` 417:                          TEMP = D[I + 1 + o_d];`
` 418:                          D[I + 1 + o_d] = DU[I + o_du] - FACT * TEMP;`
` 419:                          DU[I + o_du] = TEMP;`
` 420:                          for (J = 1; J <= NRHS; J++)`
` 421:                          {`
` 422:                              TEMP = B[I+J * LDB + o_b];`
` 423:                              B[I+J * LDB + o_b] = B[I + 1+J * LDB + o_b];`
` 424:                              B[I + 1+J * LDB + o_b] = TEMP - FACT * B[I + 1+J * LDB + o_b];`
` 425:                          }`
` 426:                      }`
` 427:                  }`
` 428:                  if (D[N + o_d] == ZERO)`
` 429:                  {`
` 430:                      INFO = N;`
` 431:                      return;`
` 432:                  }`
` 433:              }`
` 434:              // *`
` 435:              // *     Back solve with the matrix U from the factorization.`
` 436:              // *`
` 437:              if (NRHS <= 2)`
` 438:              {`
` 439:                  J = 1;`
` 440:              LABEL70:;`
` 441:                  B[N+J * LDB + o_b] = B[N+J * LDB + o_b] / D[N + o_d];`
` 442:                  if (N > 1) B[N - 1+J * LDB + o_b] = (B[N - 1+J * LDB + o_b] - DU[N - 1 + o_du] * B[N+J * LDB + o_b]) / D[N - 1 + o_d];`
` 443:                  for (I = N - 2; I >= 1; I +=  - 1)`
` 444:                  {`
` 445:                      B[I+J * LDB + o_b] = (B[I+J * LDB + o_b] - DU[I + o_du] * B[I + 1+J * LDB + o_b] - DL[I + o_dl] * B[I + 2+J * LDB + o_b]) / D[I + o_d];`
` 446:                  }`
` 447:                  if (J < NRHS)`
` 448:                  {`
` 449:                      J = J + 1;`
` 450:                      goto LABEL70;`
` 451:                  }`
` 452:              }`
` 453:              else`
` 454:              {`
` 455:                  for (J = 1; J <= NRHS; J++)`
` 456:                  {`
` 457:                      B[N+J * LDB + o_b] = B[N+J * LDB + o_b] / D[N + o_d];`
` 458:                      if (N > 1) B[N - 1+J * LDB + o_b] = (B[N - 1+J * LDB + o_b] - DU[N - 1 + o_du] * B[N+J * LDB + o_b]) / D[N - 1 + o_d];`
` 459:                      for (I = N - 2; I >= 1; I +=  - 1)`
` 460:                      {`
` 461:                          B[I+J * LDB + o_b] = (B[I+J * LDB + o_b] - DU[I + o_du] * B[I + 1+J * LDB + o_b] - DL[I + o_dl] * B[I + 2+J * LDB + o_b]) / D[I + o_d];`
` 462:                      }`
` 463:                  }`
` 464:              }`
` 465:              // *`
` 466:              return;`
` 467:              // *`
` 468:              // *     End of DGTSV`
` 469:              // *`
` 470:   `
` 471:              #endregion`
` 472:   `
` 473:          }`
` 474:      }`
` 475:  }`