`   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:      /// DGEBAL balances a general real matrix A.  This involves, first,`
`  27:      /// permuting A by a similarity transformation to isolate eigenvalues`
`  28:      /// in the first 1 to ILO-1 and last IHI+1 to N elements on the`
`  29:      /// diagonal; and second, applying a diagonal similarity transformation`
`  30:      /// to rows and columns ILO to IHI to make the rows and columns as`
`  31:      /// close in norm as possible.  Both steps are optional.`
`  32:      /// `
`  33:      /// Balancing may reduce the 1-norm of the matrix, and improve the`
`  34:      /// accuracy of the computed eigenvalues and/or eigenvectors.`
`  35:      /// `
`  36:      ///</summary>`
`  37:      public class DGEBAL`
`  38:      {`
`  39:      `
`  40:   `
`  41:          #region Dependencies`
`  42:          `
`  43:          LSAME _lsame; IDAMAX _idamax; DLAMCH _dlamch; DSCAL _dscal; DSWAP _dswap; XERBLA _xerbla; `
`  44:   `
`  45:          #endregion`
`  46:   `
`  47:   `
`  48:          #region Fields`
`  49:          `
`  50:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double SCLFAC = 2.0E+0; const double FACTOR = 0.95E+0; `
`  51:          bool NOCONV = false;int I = 0; int ICA = 0; int IEXC = 0; int IRA = 0; int J = 0; int K = 0; int L = 0; int M = 0; `
`  52:          double C = 0;double CA = 0; double F = 0; double G = 0; double R = 0; double RA = 0; double S = 0; double SFMAX1 = 0; `
`  53:          double SFMAX2 = 0;double SFMIN1 = 0; double SFMIN2 = 0; `
`  54:   `
`  55:          #endregion`
`  56:   `
`  57:          public DGEBAL(LSAME lsame, IDAMAX idamax, DLAMCH dlamch, DSCAL dscal, DSWAP dswap, XERBLA xerbla)`
`  58:          {`
`  59:      `
`  60:   `
`  61:              #region Set Dependencies`
`  62:              `
`  63:              this._lsame = lsame; this._idamax = idamax; this._dlamch = dlamch; this._dscal = dscal; this._dswap = dswap; `
`  64:              this._xerbla = xerbla;`
`  65:   `
`  66:              #endregion`
`  67:   `
`  68:          }`
`  69:      `
`  70:          public DGEBAL()`
`  71:          {`
`  72:      `
`  73:   `
`  74:              #region Dependencies (Initialization)`
`  75:              `
`  76:              LSAME lsame = new LSAME();`
`  77:              IDAMAX idamax = new IDAMAX();`
`  78:              DLAMC3 dlamc3 = new DLAMC3();`
`  79:              DSCAL dscal = new DSCAL();`
`  80:              DSWAP dswap = new DSWAP();`
`  81:              XERBLA xerbla = new XERBLA();`
`  82:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  83:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  84:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  85:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
`  86:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
`  87:   `
`  88:              #endregion`
`  89:   `
`  90:   `
`  91:              #region Set Dependencies`
`  92:              `
`  93:              this._lsame = lsame; this._idamax = idamax; this._dlamch = dlamch; this._dscal = dscal; this._dswap = dswap; `
`  94:              this._xerbla = xerbla;`
`  95:   `
`  96:              #endregion`
`  97:   `
`  98:          }`
`  99:          /// <summary>`
` 100:          /// Purpose`
` 101:          /// =======`
` 102:          /// `
` 103:          /// DGEBAL balances a general real matrix A.  This involves, first,`
` 104:          /// permuting A by a similarity transformation to isolate eigenvalues`
` 105:          /// in the first 1 to ILO-1 and last IHI+1 to N elements on the`
` 106:          /// diagonal; and second, applying a diagonal similarity transformation`
` 107:          /// to rows and columns ILO to IHI to make the rows and columns as`
` 108:          /// close in norm as possible.  Both steps are optional.`
` 109:          /// `
` 110:          /// Balancing may reduce the 1-norm of the matrix, and improve the`
` 111:          /// accuracy of the computed eigenvalues and/or eigenvectors.`
` 112:          /// `
` 113:          ///</summary>`
` 114:          /// <param name="JOB">`
` 115:          /// (input) CHARACTER*1`
` 116:          /// Specifies the operations to be performed on A:`
` 117:          /// = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0`
` 118:          /// for i = 1,...,N;`
` 119:          /// = 'P':  permute only;`
` 120:          /// = 'S':  scale only;`
` 121:          /// = 'B':  both permute and scale.`
` 122:          ///</param>`
` 123:          /// <param name="N">`
` 124:          /// (input) INTEGER`
` 125:          /// The order of the matrix A.  N .GE. 0.`
` 126:          ///</param>`
` 127:          /// <param name="A">`
` 128:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 129:          /// On entry, the input matrix A.`
` 130:          /// On exit,  A is overwritten by the balanced matrix.`
` 131:          /// If JOB = 'N', A is not referenced.`
` 132:          /// See Further Details.`
` 133:          ///</param>`
` 134:          /// <param name="LDA">`
` 135:          /// (input) INTEGER`
` 136:          /// The leading dimension of the array A.  LDA .GE. max(1,N).`
` 137:          ///</param>`
` 138:          /// <param name="ILO">`
` 139:          /// (output) INTEGER`
` 140:          ///</param>`
` 141:          /// <param name="IHI">`
` 142:          /// (output) INTEGER`
` 143:          /// ILO and IHI are set to integers such that on exit`
` 144:          /// A(i,j) = 0 if i .GT. j and j = 1,...,ILO-1 or I = IHI+1,...,N.`
` 145:          /// If JOB = 'N' or 'S', ILO = 1 and IHI = N.`
` 146:          ///</param>`
` 147:          /// <param name="SCALE">`
` 148:          /// (output) DOUBLE PRECISION array, dimension (N)`
` 149:          /// Details of the permutations and scaling factors applied to`
` 150:          /// A.  If P(j) is the index of the row and column interchanged`
` 151:          /// with row and column j and D(j) is the scaling factor`
` 152:          /// applied to row and column j, then`
` 153:          /// SCALE(j) = P(j)    for j = 1,...,ILO-1`
` 154:          /// = D(j)    for j = ILO,...,IHI`
` 155:          /// = P(j)    for j = IHI+1,...,N.`
` 156:          /// The order in which the interchanges are made is N to IHI+1,`
` 157:          /// then 1 to ILO-1.`
` 158:          ///</param>`
` 159:          /// <param name="INFO">`
` 160:          /// (output) INTEGER`
` 161:          /// = 0:  successful exit.`
` 162:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.`
` 163:          ///</param>`
` 164:          public void Run(string JOB, int N, ref double[] A, int offset_a, int LDA, ref int ILO, ref int IHI`
` 165:                           , ref double[] SCALE, int offset_scale, ref int INFO)`
` 166:          {`
` 167:   `
` 168:              #region Array Index Correction`
` 169:              `
` 170:               int o_a = -1 - LDA + offset_a;  int o_scale = -1 + offset_scale; `
` 171:   `
` 172:              #endregion`
` 173:   `
` 174:   `
` 175:              #region Strings`
` 176:              `
` 177:              JOB = JOB.Substring(0, 1);  `
` 178:   `
` 179:              #endregion`
` 180:   `
` 181:   `
` 182:              #region Prolog`
` 183:              `
` 184:              // *`
` 185:              // *  -- LAPACK routine (version 3.1) --`
` 186:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 187:              // *     November 2006`
` 188:              // *`
` 189:              // *     .. Scalar Arguments ..`
` 190:              // *     ..`
` 191:              // *     .. Array Arguments ..`
` 192:              // *     ..`
` 193:              // *`
` 194:              // *  Purpose`
` 195:              // *  =======`
` 196:              // *`
` 197:              // *  DGEBAL balances a general real matrix A.  This involves, first,`
` 198:              // *  permuting A by a similarity transformation to isolate eigenvalues`
` 199:              // *  in the first 1 to ILO-1 and last IHI+1 to N elements on the`
` 200:              // *  diagonal; and second, applying a diagonal similarity transformation`
` 201:              // *  to rows and columns ILO to IHI to make the rows and columns as`
` 202:              // *  close in norm as possible.  Both steps are optional.`
` 203:              // *`
` 204:              // *  Balancing may reduce the 1-norm of the matrix, and improve the`
` 205:              // *  accuracy of the computed eigenvalues and/or eigenvectors.`
` 206:              // *`
` 207:              // *  Arguments`
` 208:              // *  =========`
` 209:              // *`
` 210:              // *  JOB     (input) CHARACTER*1`
` 211:              // *          Specifies the operations to be performed on A:`
` 212:              // *          = 'N':  none:  simply set ILO = 1, IHI = N, SCALE(I) = 1.0`
` 213:              // *                  for i = 1,...,N;`
` 214:              // *          = 'P':  permute only;`
` 215:              // *          = 'S':  scale only;`
` 216:              // *          = 'B':  both permute and scale.`
` 217:              // *`
` 218:              // *  N       (input) INTEGER`
` 219:              // *          The order of the matrix A.  N >= 0.`
` 220:              // *`
` 221:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)`
` 222:              // *          On entry, the input matrix A.`
` 223:              // *          On exit,  A is overwritten by the balanced matrix.`
` 224:              // *          If JOB = 'N', A is not referenced.`
` 225:              // *          See Further Details.`
` 226:              // *`
` 227:              // *  LDA     (input) INTEGER`
` 228:              // *          The leading dimension of the array A.  LDA >= max(1,N).`
` 229:              // *`
` 230:              // *  ILO     (output) INTEGER`
` 231:              // *  IHI     (output) INTEGER`
` 232:              // *          ILO and IHI are set to integers such that on exit`
` 233:              // *          A(i,j) = 0 if i > j and j = 1,...,ILO-1 or I = IHI+1,...,N.`
` 234:              // *          If JOB = 'N' or 'S', ILO = 1 and IHI = N.`
` 235:              // *`
` 236:              // *  SCALE   (output) DOUBLE PRECISION array, dimension (N)`
` 237:              // *          Details of the permutations and scaling factors applied to`
` 238:              // *          A.  If P(j) is the index of the row and column interchanged`
` 239:              // *          with row and column j and D(j) is the scaling factor`
` 240:              // *          applied to row and column j, then`
` 241:              // *          SCALE(j) = P(j)    for j = 1,...,ILO-1`
` 242:              // *                   = D(j)    for j = ILO,...,IHI`
` 243:              // *                   = P(j)    for j = IHI+1,...,N.`
` 244:              // *          The order in which the interchanges are made is N to IHI+1,`
` 245:              // *          then 1 to ILO-1.`
` 246:              // *`
` 247:              // *  INFO    (output) INTEGER`
` 248:              // *          = 0:  successful exit.`
` 249:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.`
` 250:              // *`
` 251:              // *  Further Details`
` 252:              // *  ===============`
` 253:              // *`
` 254:              // *  The permutations consist of row and column interchanges which put`
` 255:              // *  the matrix in the form`
` 256:              // *`
` 257:              // *             ( T1   X   Y  )`
` 258:              // *     P A P = (  0   B   Z  )`
` 259:              // *             (  0   0   T2 )`
` 260:              // *`
` 261:              // *  where T1 and T2 are upper triangular matrices whose eigenvalues lie`
` 262:              // *  along the diagonal.  The column indices ILO and IHI mark the starting`
` 263:              // *  and ending columns of the submatrix B. Balancing consists of applying`
` 264:              // *  a diagonal similarity transformation inv(D) * B * D to make the`
` 265:              // *  1-norms of each row of B and its corresponding column nearly equal.`
` 266:              // *  The output matrix is`
` 267:              // *`
` 268:              // *     ( T1     X*D          Y    )`
` 269:              // *     (  0  inv(D)*B*D  inv(D)*Z ).`
` 270:              // *     (  0      0           T2   )`
` 271:              // *`
` 272:              // *  Information about the permutations P and the diagonal matrix D is`
` 273:              // *  returned in the vector SCALE.`
` 274:              // *`
` 275:              // *  This subroutine is based on the EISPACK routine BALANC.`
` 276:              // *`
` 277:              // *  Modified by Tzu-Yi Chen, Computer Science Division, University of`
` 278:              // *    California at Berkeley, USA`
` 279:              // *`
` 280:              // *  =====================================================================`
` 281:              // *`
` 282:              // *     .. Parameters ..`
` 283:              // *     ..`
` 284:              // *     .. Local Scalars ..`
` 285:              // *     ..`
` 286:              // *     .. External Functions ..`
` 287:              // *     ..`
` 288:              // *     .. External Subroutines ..`
` 289:              // *     ..`
` 290:              // *     .. Intrinsic Functions ..`
` 291:              //      INTRINSIC          ABS, MAX, MIN;`
` 292:              // *     ..`
` 293:              // *     .. Executable Statements ..`
` 294:              // *`
` 295:              // *     Test the input parameters`
` 296:              // *`
` 297:   `
` 298:              #endregion`
` 299:   `
` 300:   `
` 301:              #region Body`
` 302:              `
` 303:              INFO = 0;`
` 304:              if (!this._lsame.Run(JOB, "N") && !this._lsame.Run(JOB, "P") && !this._lsame.Run(JOB, "S") && !this._lsame.Run(JOB, "B"))`
` 305:              {`
` 306:                  INFO =  - 1;`
` 307:              }`
` 308:              else`
` 309:              {`
` 310:                  if (N < 0)`
` 311:                  {`
` 312:                      INFO =  - 2;`
` 313:                  }`
` 314:                  else`
` 315:                  {`
` 316:                      if (LDA < Math.Max(1, N))`
` 317:                      {`
` 318:                          INFO =  - 4;`
` 319:                      }`
` 320:                  }`
` 321:              }`
` 322:              if (INFO != 0)`
` 323:              {`
` 324:                  this._xerbla.Run("DGEBAL",  - INFO);`
` 325:                  return;`
` 326:              }`
` 327:              // *`
` 328:              K = 1;`
` 329:              L = N;`
` 330:              // *`
` 331:              if (N == 0) goto LABEL210;`
` 332:              // *`
` 333:              if (this._lsame.Run(JOB, "N"))`
` 334:              {`
` 335:                  for (I = 1; I <= N; I++)`
` 336:                  {`
` 337:                      SCALE[I + o_scale] = ONE;`
` 338:                  }`
` 339:                  goto LABEL210;`
` 340:              }`
` 341:              // *`
` 342:              if (this._lsame.Run(JOB, "S")) goto LABEL120;`
` 343:              // *`
` 344:              // *     Permutation to isolate eigenvalues if possible`
` 345:              // *`
` 346:              goto LABEL50;`
` 347:              // *`
` 348:              // *     Row and column exchange.`
` 349:              // *`
` 350:          LABEL20:;`
` 351:              SCALE[M + o_scale] = J;`
` 352:              if (J == M) goto LABEL30;`
` 353:              // *`
` 354:              this._dswap.Run(L, ref A, 1+J * LDA + o_a, 1, ref A, 1+M * LDA + o_a, 1);`
` 355:              this._dswap.Run(N - K + 1, ref A, J+K * LDA + o_a, LDA, ref A, M+K * LDA + o_a, LDA);`
` 356:              // *`
` 357:          LABEL30:;`
` 358:              switch (IEXC)`
` 359:              {`
` 360:                  case 1: goto LABEL40;`
` 361:                  case 2: goto LABEL80;`
` 362:              }`
` 363:              // *`
` 364:              // *     Search for rows isolating an eigenvalue and push them down.`
` 365:              // *`
` 366:          LABEL40:;`
` 367:              if (L == 1) goto LABEL210;`
` 368:              L = L - 1;`
` 369:              // *`
` 370:          LABEL50:;`
` 371:              for (J = L; J >= 1; J +=  - 1)`
` 372:              {`
` 373:                  // *`
` 374:                  for (I = 1; I <= L; I++)`
` 375:                  {`
` 376:                      if (I == J) goto LABEL60;`
` 377:                      if (A[J+I * LDA + o_a] != ZERO) goto LABEL70;`
` 378:                  LABEL60:;`
` 379:                  }`
` 380:                  // *`
` 381:                  M = L;`
` 382:                  IEXC = 1;`
` 383:                  goto LABEL20;`
` 384:              LABEL70:;`
` 385:              }`
` 386:              // *`
` 387:              goto LABEL90;`
` 388:              // *`
` 389:              // *     Search for columns isolating an eigenvalue and push them left.`
` 390:              // *`
` 391:          LABEL80:;`
` 392:              K = K + 1;`
` 393:              // *`
` 394:          LABEL90:;`
` 395:              for (J = K; J <= L; J++)`
` 396:              {`
` 397:                  // *`
` 398:                  for (I = K; I <= L; I++)`
` 399:                  {`
` 400:                      if (I == J) goto LABEL100;`
` 401:                      if (A[I+J * LDA + o_a] != ZERO) goto LABEL110;`
` 402:                  LABEL100:;`
` 403:                  }`
` 404:                  // *`
` 405:                  M = K;`
` 406:                  IEXC = 2;`
` 407:                  goto LABEL20;`
` 408:              LABEL110:;`
` 409:              }`
` 410:              // *`
` 411:          LABEL120:;`
` 412:              for (I = K; I <= L; I++)`
` 413:              {`
` 414:                  SCALE[I + o_scale] = ONE;`
` 415:              }`
` 416:              // *`
` 417:              if (this._lsame.Run(JOB, "P")) goto LABEL210;`
` 418:              // *`
` 419:              // *     Balance the submatrix in rows K to L.`
` 420:              // *`
` 421:              // *     Iterative loop for norm reduction`
` 422:              // *`
` 423:              SFMIN1 = this._dlamch.Run("S") / this._dlamch.Run("P");`
` 424:              SFMAX1 = ONE / SFMIN1;`
` 425:              SFMIN2 = SFMIN1 * SCLFAC;`
` 426:              SFMAX2 = ONE / SFMIN2;`
` 427:          LABEL140:;`
` 428:              NOCONV = false;`
` 429:              // *`
` 430:              for (I = K; I <= L; I++)`
` 431:              {`
` 432:                  C = ZERO;`
` 433:                  R = ZERO;`
` 434:                  // *`
` 435:                  for (J = K; J <= L; J++)`
` 436:                  {`
` 437:                      if (J == I) goto LABEL150;`
` 438:                      C = C + Math.Abs(A[J+I * LDA + o_a]);`
` 439:                      R = R + Math.Abs(A[I+J * LDA + o_a]);`
` 440:                  LABEL150:;`
` 441:                  }`
` 442:                  ICA = this._idamax.Run(L, A, 1+I * LDA + o_a, 1);`
` 443:                  CA = Math.Abs(A[ICA+I * LDA + o_a]);`
` 444:                  IRA = this._idamax.Run(N - K + 1, A, I+K * LDA + o_a, LDA);`
` 445:                  RA = Math.Abs(A[I+(IRA + K - 1) * LDA + o_a]);`
` 446:                  // *`
` 447:                  // *        Guard against zero C or R due to underflow.`
` 448:                  // *`
` 449:                  if (C == ZERO || R == ZERO) goto LABEL200;`
` 450:                  G = R / SCLFAC;`
` 451:                  F = ONE;`
` 452:                  S = C + R;`
` 453:              LABEL160:;`
` 454:                  if (C >= G || Math.Max(F, Math.Max(C, CA)) >= SFMAX2 || Math.Min(R, Math.Min(G, RA)) <= SFMIN2) goto LABEL170;`
` 455:                  F = F * SCLFAC;`
` 456:                  C = C * SCLFAC;`
` 457:                  CA = CA * SCLFAC;`
` 458:                  R = R / SCLFAC;`
` 459:                  G = G / SCLFAC;`
` 460:                  RA = RA / SCLFAC;`
` 461:                  goto LABEL160;`
` 462:                  // *`
` 463:              LABEL170:;`
` 464:                  G = C / SCLFAC;`
` 465:              LABEL180:;`
` 466:                  if (G < R || Math.Max(R, RA) >= SFMAX2 || Math.Min(F, Math.Min(C, Math.Min(G, CA))) <= SFMIN2) goto LABEL190;`
` 467:                  F = F / SCLFAC;`
` 468:                  C = C / SCLFAC;`
` 469:                  G = G / SCLFAC;`
` 470:                  CA = CA / SCLFAC;`
` 471:                  R = R * SCLFAC;`
` 472:                  RA = RA * SCLFAC;`
` 473:                  goto LABEL180;`
` 474:                  // *`
` 475:                  // *        Now balance.`
` 476:                  // *`
` 477:              LABEL190:;`
` 478:                  if ((C + R) >= FACTOR * S) goto LABEL200;`
` 479:                  if (F < ONE && SCALE[I + o_scale] < ONE)`
` 480:                  {`
` 481:                      if (F * SCALE[I + o_scale] <= SFMIN1) goto LABEL200;`
` 482:                  }`
` 483:                  if (F > ONE && SCALE[I + o_scale] > ONE)`
` 484:                  {`
` 485:                      if (SCALE[I + o_scale] >= SFMAX1 / F) goto LABEL200;`
` 486:                  }`
` 487:                  G = ONE / F;`
` 488:                  SCALE[I + o_scale] = SCALE[I + o_scale] * F;`
` 489:                  NOCONV = true;`
` 490:                  // *`
` 491:                  this._dscal.Run(N - K + 1, G, ref A, I+K * LDA + o_a, LDA);`
` 492:                  this._dscal.Run(L, F, ref A, 1+I * LDA + o_a, 1);`
` 493:                  // *`
` 494:              LABEL200:;`
` 495:              }`
` 496:              // *`
` 497:              if (NOCONV) goto LABEL140;`
` 498:              // *`
` 499:          LABEL210:;`
` 500:              ILO = K;`
` 501:              IHI = L;`
` 502:              // *`
` 503:              return;`
` 504:              // *`
` 505:              // *     End of DGEBAL`
` 506:              // *`
` 507:   `
` 508:              #endregion`
` 509:   `
` 510:          }`
` 511:      }`
` 512:  }`