`   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:      /// DLASQ2 computes all the eigenvalues of the symmetric positive `
`  27:      /// definite tridiagonal matrix associated with the qd array Z to high`
`  28:      /// relative accuracy are computed to high relative accuracy, in the`
`  29:      /// absence of denormalization, underflow and overflow.`
`  30:      /// `
`  31:      /// To see the relation of Z to the tridiagonal matrix, let L be a`
`  32:      /// unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and`
`  33:      /// let U be an upper bidiagonal matrix with 1's above and diagonal`
`  34:      /// Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the`
`  35:      /// symmetric tridiagonal to which it is similar.`
`  36:      /// `
`  37:      /// Note : DLASQ2 defines a logical variable, IEEE, which is true`
`  38:      /// on machines which follow ieee-754 floating-point standard in their`
`  39:      /// handling of infinities and NaNs, and false otherwise. This variable`
`  40:      /// is passed to DLAZQ3.`
`  41:      /// `
`  42:      ///</summary>`
`  43:      public class DLASQ2`
`  44:      {`
`  45:      `
`  46:   `
`  47:          #region Dependencies`
`  48:          `
`  49:          DLAZQ3 _dlazq3; DLASRT _dlasrt; XERBLA _xerbla; DLAMCH _dlamch; ILAENV _ilaenv; `
`  50:   `
`  51:          #endregion`
`  52:   `
`  53:   `
`  54:          #region Fields`
`  55:          `
`  56:          const double CBIAS = 1.50E0; const double ZERO = 0.0E0; const double HALF = 0.5E0; const double ONE = 1.0E0; `
`  57:          const double TWO = 2.0E0;const double FOUR = 4.0E0; const double HUNDRD = 100.0E0; bool IEEE = false; int I0 = 0; `
`  58:          int I4 = 0;int IINFO = 0; int IPN4 = 0; int ITER = 0; int IWHILA = 0; int IWHILB = 0; int K = 0; int N0 = 0; int NBIG = 0; `
`  59:          int NDIV = 0;int NFAIL = 0; int PP = 0; int SPLT = 0; int TTYPE = 0; double D = 0; double DESIG = 0; double DMIN = 0; `
`  60:          double DMIN1 = 0;double DMIN2 = 0; double DN = 0; double DN1 = 0; double DN2 = 0; double E = 0; double EMAX = 0; `
`  61:          double EMIN = 0;double EPS = 0; double OLDEMN = 0; double QMAX = 0; double QMIN = 0; double S = 0; double SAFMIN = 0; `
`  62:          double SIGMA = 0;double T = 0; double TAU = 0; double TEMP = 0; double TOL = 0; double TOL2 = 0; double TRACE = 0; `
`  63:          double ZMAX = 0;`
`  64:   `
`  65:          #endregion`
`  66:   `
`  67:          public DLASQ2(DLAZQ3 dlazq3, DLASRT dlasrt, XERBLA xerbla, DLAMCH dlamch, ILAENV ilaenv)`
`  68:          {`
`  69:      `
`  70:   `
`  71:              #region Set Dependencies`
`  72:              `
`  73:              this._dlazq3 = dlazq3; this._dlasrt = dlasrt; this._xerbla = xerbla; this._dlamch = dlamch; this._ilaenv = ilaenv; `
`  74:   `
`  75:              #endregion`
`  76:   `
`  77:          }`
`  78:      `
`  79:          public DLASQ2()`
`  80:          {`
`  81:      `
`  82:   `
`  83:              #region Dependencies (Initialization)`
`  84:              `
`  85:              DLASQ5 dlasq5 = new DLASQ5();`
`  86:              LSAME lsame = new LSAME();`
`  87:              DLAMC3 dlamc3 = new DLAMC3();`
`  88:              DLAZQ4 dlazq4 = new DLAZQ4();`
`  89:              XERBLA xerbla = new XERBLA();`
`  90:              IEEECK ieeeck = new IEEECK();`
`  91:              IPARMQ iparmq = new IPARMQ();`
`  92:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  93:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  94:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  95:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
`  96:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
`  97:              DLASQ6 dlasq6 = new DLASQ6(dlamch);`
`  98:              DLAZQ3 dlazq3 = new DLAZQ3(dlasq5, dlasq6, dlazq4, dlamch);`
`  99:              DLASRT dlasrt = new DLASRT(lsame, xerbla);`
` 100:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);`
` 101:   `
` 102:              #endregion`
` 103:   `
` 104:   `
` 105:              #region Set Dependencies`
` 106:              `
` 107:              this._dlazq3 = dlazq3; this._dlasrt = dlasrt; this._xerbla = xerbla; this._dlamch = dlamch; this._ilaenv = ilaenv; `
` 108:   `
` 109:              #endregion`
` 110:   `
` 111:          }`
` 112:          /// <summary>`
` 113:          /// Purpose`
` 114:          /// =======`
` 115:          /// `
` 116:          /// DLASQ2 computes all the eigenvalues of the symmetric positive `
` 117:          /// definite tridiagonal matrix associated with the qd array Z to high`
` 118:          /// relative accuracy are computed to high relative accuracy, in the`
` 119:          /// absence of denormalization, underflow and overflow.`
` 120:          /// `
` 121:          /// To see the relation of Z to the tridiagonal matrix, let L be a`
` 122:          /// unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and`
` 123:          /// let U be an upper bidiagonal matrix with 1's above and diagonal`
` 124:          /// Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the`
` 125:          /// symmetric tridiagonal to which it is similar.`
` 126:          /// `
` 127:          /// Note : DLASQ2 defines a logical variable, IEEE, which is true`
` 128:          /// on machines which follow ieee-754 floating-point standard in their`
` 129:          /// handling of infinities and NaNs, and false otherwise. This variable`
` 130:          /// is passed to DLAZQ3.`
` 131:          /// `
` 132:          ///</summary>`
` 133:          /// <param name="N">`
` 134:          /// (input) INTEGER`
` 135:          /// The number of rows and columns in the matrix. N .GE. 0.`
` 136:          ///</param>`
` 137:          /// <param name="Z">`
` 138:          /// (workspace) DOUBLE PRECISION array, dimension ( 4*N )`
` 139:          /// On entry Z holds the qd array. On exit, entries 1 to N hold`
` 140:          /// the eigenvalues in decreasing order, Z( 2*N+1 ) holds the`
` 141:          /// trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If`
` 142:          /// N .GT. 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )`
` 143:          /// holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of`
` 144:          /// shifts that failed.`
` 145:          ///</param>`
` 146:          /// <param name="INFO">`
` 147:          /// (output) INTEGER`
` 148:          /// = 0: successful exit`
` 149:          /// .LT. 0: if the i-th argument is a scalar and had an illegal`
` 150:          /// value, then INFO = -i, if the i-th argument is an`
` 151:          /// array and the j-entry had an illegal value, then`
` 152:          /// INFO = -(i*100+j)`
` 153:          /// .GT. 0: the algorithm failed`
` 154:          /// = 1, a split was marked by a positive value in E`
` 155:          /// = 2, current block of Z not diagonalized after 30*N`
` 156:          /// iterations (in inner while loop)`
` 157:          /// = 3, termination criterion of outer while loop not met `
` 158:          /// (program created more than N unreduced blocks)`
` 159:          ///</param>`
` 160:          public void Run(int N, ref double[] Z, int offset_z, ref int INFO)`
` 161:          {`
` 162:   `
` 163:              #region Array Index Correction`
` 164:              `
` 165:               int o_z = -1 + offset_z; `
` 166:   `
` 167:              #endregion`
` 168:   `
` 169:   `
` 170:              #region Prolog`
` 171:              `
` 172:              // *`
` 173:              // *  -- LAPACK routine (version 3.1) --`
` 174:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 175:              // *     November 2006`
` 176:              // *`
` 177:              // *     Modified to call DLAZQ3 in place of DLASQ3, 13 Feb 03, SJH.`
` 178:              // *`
` 179:              // *     .. Scalar Arguments ..`
` 180:              // *     ..`
` 181:              // *     .. Array Arguments ..`
` 182:              // *     ..`
` 183:              // *`
` 184:              // *  Purpose`
` 185:              // *  =======`
` 186:              // *`
` 187:              // *  DLASQ2 computes all the eigenvalues of the symmetric positive `
` 188:              // *  definite tridiagonal matrix associated with the qd array Z to high`
` 189:              // *  relative accuracy are computed to high relative accuracy, in the`
` 190:              // *  absence of denormalization, underflow and overflow.`
` 191:              // *`
` 192:              // *  To see the relation of Z to the tridiagonal matrix, let L be a`
` 193:              // *  unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and`
` 194:              // *  let U be an upper bidiagonal matrix with 1's above and diagonal`
` 195:              // *  Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the`
` 196:              // *  symmetric tridiagonal to which it is similar.`
` 197:              // *`
` 198:              // *  Note : DLASQ2 defines a logical variable, IEEE, which is true`
` 199:              // *  on machines which follow ieee-754 floating-point standard in their`
` 200:              // *  handling of infinities and NaNs, and false otherwise. This variable`
` 201:              // *  is passed to DLAZQ3.`
` 202:              // *`
` 203:              // *  Arguments`
` 204:              // *  =========`
` 205:              // *`
` 206:              // *  N     (input) INTEGER`
` 207:              // *        The number of rows and columns in the matrix. N >= 0.`
` 208:              // *`
` 209:              // *  Z     (workspace) DOUBLE PRECISION array, dimension ( 4*N )`
` 210:              // *        On entry Z holds the qd array. On exit, entries 1 to N hold`
` 211:              // *        the eigenvalues in decreasing order, Z( 2*N+1 ) holds the`
` 212:              // *        trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If`
` 213:              // *        N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 )`
` 214:              // *        holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of`
` 215:              // *        shifts that failed.`
` 216:              // *`
` 217:              // *  INFO  (output) INTEGER`
` 218:              // *        = 0: successful exit`
` 219:              // *        < 0: if the i-th argument is a scalar and had an illegal`
` 220:              // *             value, then INFO = -i, if the i-th argument is an`
` 221:              // *             array and the j-entry had an illegal value, then`
` 222:              // *             INFO = -(i*100+j)`
` 223:              // *        > 0: the algorithm failed`
` 224:              // *              = 1, a split was marked by a positive value in E`
` 225:              // *              = 2, current block of Z not diagonalized after 30*N`
` 226:              // *                   iterations (in inner while loop)`
` 227:              // *              = 3, termination criterion of outer while loop not met `
` 228:              // *                   (program created more than N unreduced blocks)`
` 229:              // *`
` 230:              // *  Further Details`
` 231:              // *  ===============`
` 232:              // *  Local Variables: I0:N0 defines a current unreduced segment of Z.`
` 233:              // *  The shifts are accumulated in SIGMA. Iteration count is in ITER.`
` 234:              // *  Ping-pong is controlled by PP (alternates between 0 and 1).`
` 235:              // *`
` 236:              // *  =====================================================================`
` 237:              // *`
` 238:              // *     .. Parameters ..`
` 239:              // *     ..`
` 240:              // *     .. Local Scalars ..`
` 241:              // *     ..`
` 242:              // *     .. External Subroutines ..`
` 243:              // *     ..`
` 244:              // *     .. External Functions ..`
` 245:              // *     ..`
` 246:              // *     .. Intrinsic Functions ..`
` 247:              //      INTRINSIC          ABS, DBLE, MAX, MIN, SQRT;`
` 248:              // *     ..`
` 249:              // *     .. Executable Statements ..`
` 250:              // *      `
` 251:              // *     Test the input arguments.`
` 252:              // *     (in case DLASQ2 is not called by DLASQ1)`
` 253:              // *`
` 254:   `
` 255:              #endregion`
` 256:   `
` 257:   `
` 258:              #region Body`
` 259:              `
` 260:              INFO = 0;`
` 261:              EPS = this._dlamch.Run("Precision");`
` 262:              SAFMIN = this._dlamch.Run("Safe minimum");`
` 263:              TOL = EPS * HUNDRD;`
` 264:              TOL2 = Math.Pow(TOL,2);`
` 265:              // *`
` 266:              if (N < 0)`
` 267:              {`
` 268:                  INFO =  - 1;`
` 269:                  this._xerbla.Run("DLASQ2", 1);`
` 270:                  return;`
` 271:              }`
` 272:              else`
` 273:              {`
` 274:                  if (N == 0)`
` 275:                  {`
` 276:                      return;`
` 277:                  }`
` 278:                  else`
` 279:                  {`
` 280:                      if (N == 1)`
` 281:                      {`
` 282:                          // *`
` 283:                          // *        1-by-1 case.`
` 284:                          // *`
` 285:                          if (Z[1 + o_z] < ZERO)`
` 286:                          {`
` 287:                              INFO =  - 201;`
` 288:                              this._xerbla.Run("DLASQ2", 2);`
` 289:                          }`
` 290:                          return;`
` 291:                      }`
` 292:                      else`
` 293:                      {`
` 294:                          if (N == 2)`
` 295:                          {`
` 296:                              // *`
` 297:                              // *        2-by-2 case.`
` 298:                              // *`
` 299:                              if (Z[2 + o_z] < ZERO || Z[3 + o_z] < ZERO)`
` 300:                              {`
` 301:                                  INFO =  - 2;`
` 302:                                  this._xerbla.Run("DLASQ2", 2);`
` 303:                                  return;`
` 304:                              }`
` 305:                              else`
` 306:                              {`
` 307:                                  if (Z[3 + o_z] > Z[1 + o_z])`
` 308:                                  {`
` 309:                                      D = Z[3 + o_z];`
` 310:                                      Z[3 + o_z] = Z[1 + o_z];`
` 311:                                      Z[1 + o_z] = D;`
` 312:                                  }`
` 313:                              }`
` 314:                              Z[5 + o_z] = Z[1 + o_z] + Z[2 + o_z] + Z[3 + o_z];`
` 315:                              if (Z[2 + o_z] > Z[3 + o_z] * TOL2)`
` 316:                              {`
` 317:                                  T = HALF * ((Z[1 + o_z] - Z[3 + o_z]) + Z[2 + o_z]);`
` 318:                                  S = Z[3 + o_z] * (Z[2 + o_z] / T);`
` 319:                                  if (S <= T)`
` 320:                                  {`
` 321:                                      S = Z[3 + o_z] * (Z[2 + o_z] / (T * (ONE + Math.Sqrt(ONE + S / T))));`
` 322:                                  }`
` 323:                                  else`
` 324:                                  {`
` 325:                                      S = Z[3 + o_z] * (Z[2 + o_z] / (T + Math.Sqrt(T) * Math.Sqrt(T + S)));`
` 326:                                  }`
` 327:                                  T = Z[1 + o_z] + (S + Z[2 + o_z]);`
` 328:                                  Z[3 + o_z] = Z[3 + o_z] * (Z[1 + o_z] / T);`
` 329:                                  Z[1 + o_z] = T;`
` 330:                              }`
` 331:                              Z[2 + o_z] = Z[3 + o_z];`
` 332:                              Z[6 + o_z] = Z[2 + o_z] + Z[1 + o_z];`
` 333:                              return;`
` 334:                          }`
` 335:                      }`
` 336:                  }`
` 337:              }`
` 338:              // *`
` 339:              // *     Check for negative data and compute sums of q's and e's.`
` 340:              // *`
` 341:              Z[2 * N + o_z] = ZERO;`
` 342:              EMIN = Z[2 + o_z];`
` 343:              QMAX = ZERO;`
` 344:              ZMAX = ZERO;`
` 345:              D = ZERO;`
` 346:              E = ZERO;`
` 347:              // *`
` 348:              for (K = 1; K <= 2 * (N - 1); K += 2)`
` 349:              {`
` 350:                  if (Z[K + o_z] < ZERO)`
` 351:                  {`
` 352:                      INFO =  - (200 + K);`
` 353:                      this._xerbla.Run("DLASQ2", 2);`
` 354:                      return;`
` 355:                  }`
` 356:                  else`
` 357:                  {`
` 358:                      if (Z[K + 1 + o_z] < ZERO)`
` 359:                      {`
` 360:                          INFO =  - (200 + K + 1);`
` 361:                          this._xerbla.Run("DLASQ2", 2);`
` 362:                          return;`
` 363:                      }`
` 364:                  }`
` 365:                  D = D + Z[K + o_z];`
` 366:                  E = E + Z[K + 1 + o_z];`
` 367:                  QMAX = Math.Max(QMAX, Z[K + o_z]);`
` 368:                  EMIN = Math.Min(EMIN, Z[K + 1 + o_z]);`
` 369:                  ZMAX = Math.Max(QMAX, Math.Max(ZMAX, Z[K + 1 + o_z]));`
` 370:              }`
` 371:              if (Z[2 * N - 1 + o_z] < ZERO)`
` 372:              {`
` 373:                  INFO =  - (200 + 2 * N - 1);`
` 374:                  this._xerbla.Run("DLASQ2", 2);`
` 375:                  return;`
` 376:              }`
` 377:              D = D + Z[2 * N - 1 + o_z];`
` 378:              QMAX = Math.Max(QMAX, Z[2 * N - 1 + o_z]);`
` 379:              ZMAX = Math.Max(QMAX, ZMAX);`
` 380:              // *`
` 381:              // *     Check for diagonality.`
` 382:              // *`
` 383:              if (E == ZERO)`
` 384:              {`
` 385:                  for (K = 2; K <= N; K++)`
` 386:                  {`
` 387:                      Z[K + o_z] = Z[2 * K - 1 + o_z];`
` 388:                  }`
` 389:                  this._dlasrt.Run("D", N, ref Z, offset_z, ref IINFO);`
` 390:                  Z[2 * N - 1 + o_z] = D;`
` 391:                  return;`
` 392:              }`
` 393:              // *`
` 394:              TRACE = D + E;`
` 395:              // *`
` 396:              // *     Check for zero data.`
` 397:              // *`
` 398:              if (TRACE == ZERO)`
` 399:              {`
` 400:                  Z[2 * N - 1 + o_z] = ZERO;`
` 401:                  return;`
` 402:              }`
` 403:              // *         `
` 404:              // *     Check whether the machine is IEEE conformable.`
` 405:              // *         `
` 406:              IEEE = this._ilaenv.Run(10, "DLASQ2", "N", 1, 2, 3, 4) == 1 && this._ilaenv.Run(11, "DLASQ2", "N", 1, 2, 3, 4) == 1;`
` 407:              // *         `
` 408:              // *     Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).`
` 409:              // *`
` 410:              for (K = 2 * N; K >= 2; K +=  - 2)`
` 411:              {`
` 412:                  Z[2 * K + o_z] = ZERO;`
` 413:                  Z[2 * K - 1 + o_z] = Z[K + o_z];`
` 414:                  Z[2 * K - 2 + o_z] = ZERO;`
` 415:                  Z[2 * K - 3 + o_z] = Z[K - 1 + o_z];`
` 416:              }`
` 417:              // *`
` 418:              I0 = 1;`
` 419:              N0 = N;`
` 420:              // *`
` 421:              // *     Reverse the qd-array, if warranted.`
` 422:              // *`
` 423:              if (CBIAS * Z[4 * I0 - 3 + o_z] < Z[4 * N0 - 3 + o_z])`
` 424:              {`
` 425:                  IPN4 = 4 * (I0 + N0);`
` 426:                  for (I4 = 4 * I0; I4 <= 2 * (I0 + N0 - 1); I4 += 4)`
` 427:                  {`
` 428:                      TEMP = Z[I4 - 3 + o_z];`
` 429:                      Z[I4 - 3 + o_z] = Z[IPN4 - I4 - 3 + o_z];`
` 430:                      Z[IPN4 - I4 - 3 + o_z] = TEMP;`
` 431:                      TEMP = Z[I4 - 1 + o_z];`
` 432:                      Z[I4 - 1 + o_z] = Z[IPN4 - I4 - 5 + o_z];`
` 433:                      Z[IPN4 - I4 - 5 + o_z] = TEMP;`
` 434:                  }`
` 435:              }`
` 436:              // *`
` 437:              // *     Initial split checking via dqd and Li's test.`
` 438:              // *`
` 439:              PP = 0;`
` 440:              // *`
` 441:              for (K = 1; K <= 2; K++)`
` 442:              {`
` 443:                  // *`
` 444:                  D = Z[4 * N0 + PP - 3 + o_z];`
` 445:                  for (I4 = 4 * (N0 - 1) + PP; I4 >= 4 * I0 + PP; I4 +=  - 4)`
` 446:                  {`
` 447:                      if (Z[I4 - 1 + o_z] <= TOL2 * D)`
` 448:                      {`
` 449:                          Z[I4 - 1 + o_z] =  - ZERO;`
` 450:                          D = Z[I4 - 3 + o_z];`
` 451:                      }`
` 452:                      else`
` 453:                      {`
` 454:                          D = Z[I4 - 3 + o_z] * (D / (D + Z[I4 - 1 + o_z]));`
` 455:                      }`
` 456:                  }`
` 457:                  // *`
` 458:                  // *        dqd maps Z to ZZ plus Li's test.`
` 459:                  // *`
` 460:                  EMIN = Z[4 * I0 + PP + 1 + o_z];`
` 461:                  D = Z[4 * I0 + PP - 3 + o_z];`
` 462:                  for (I4 = 4 * I0 + PP; I4 <= 4 * (N0 - 1) + PP; I4 += 4)`
` 463:                  {`
` 464:                      Z[I4 - 2 * PP - 2 + o_z] = D + Z[I4 - 1 + o_z];`
` 465:                      if (Z[I4 - 1 + o_z] <= TOL2 * D)`
` 466:                      {`
` 467:                          Z[I4 - 1 + o_z] =  - ZERO;`
` 468:                          Z[I4 - 2 * PP - 2 + o_z] = D;`
` 469:                          Z[I4 - 2 * PP + o_z] = ZERO;`
` 470:                          D = Z[I4 + 1 + o_z];`
` 471:                      }`
` 472:                      else`
` 473:                      {`
` 474:                          if (SAFMIN * Z[I4 + 1 + o_z] < Z[I4 - 2 * PP - 2 + o_z] && SAFMIN * Z[I4 - 2 * PP - 2 + o_z] < Z[I4 + 1 + o_z])`
` 475:                          {`
` 476:                              TEMP = Z[I4 + 1 + o_z] / Z[I4 - 2 * PP - 2 + o_z];`
` 477:                              Z[I4 - 2 * PP + o_z] = Z[I4 - 1 + o_z] * TEMP;`
` 478:                              D = D * TEMP;`
` 479:                          }`
` 480:                          else`
` 481:                          {`
` 482:                              Z[I4 - 2 * PP + o_z] = Z[I4 + 1 + o_z] * (Z[I4 - 1 + o_z] / Z[I4 - 2 * PP - 2 + o_z]);`
` 483:                              D = Z[I4 + 1 + o_z] * (D / Z[I4 - 2 * PP - 2 + o_z]);`
` 484:                          }`
` 485:                      }`
` 486:                      EMIN = Math.Min(EMIN, Z[I4 - 2 * PP + o_z]);`
` 487:                  }`
` 488:                  Z[4 * N0 - PP - 2 + o_z] = D;`
` 489:                  // *`
` 490:                  // *        Now find qmax.`
` 491:                  // *`
` 492:                  QMAX = Z[4 * I0 - PP - 2 + o_z];`
` 493:                  for (I4 = 4 * I0 - PP + 2; I4 <= 4 * N0 - PP - 2; I4 += 4)`
` 494:                  {`
` 495:                      QMAX = Math.Max(QMAX, Z[I4 + o_z]);`
` 496:                  }`
` 497:                  // *`
` 498:                  // *        Prepare for the next iteration on K.`
` 499:                  // *`
` 500:                  PP = 1 - PP;`
` 501:              }`
` 502:              // *`
` 503:              // *     Initialise variables to pass to DLAZQ3`
` 504:              // *`
` 505:              TTYPE = 0;`
` 506:              DMIN1 = ZERO;`
` 507:              DMIN2 = ZERO;`
` 508:              DN = ZERO;`
` 509:              DN1 = ZERO;`
` 510:              DN2 = ZERO;`
` 511:              TAU = ZERO;`
` 512:              // *`
` 513:              ITER = 2;`
` 514:              NFAIL = 0;`
` 515:              NDIV = 2 * (N0 - I0);`
` 516:              // *`
` 517:              for (IWHILA = 1; IWHILA <= N + 1; IWHILA++)`
` 518:              {`
` 519:                  if (N0 < 1) goto LABEL150;`
` 520:                  // *`
` 521:                  // *        While array unfinished do `
` 522:                  // *`
` 523:                  // *        E(N0) holds the value of SIGMA when submatrix in I0:N0`
` 524:                  // *        splits from the rest of the array, but is negated.`
` 525:                  // *      `
` 526:                  DESIG = ZERO;`
` 527:                  if (N0 == N)`
` 528:                  {`
` 529:                      SIGMA = ZERO;`
` 530:                  }`
` 531:                  else`
` 532:                  {`
` 533:                      SIGMA =  - Z[4 * N0 - 1 + o_z];`
` 534:                  }`
` 535:                  if (SIGMA < ZERO)`
` 536:                  {`
` 537:                      INFO = 1;`
` 538:                      return;`
` 539:                  }`
` 540:                  // *`
` 541:                  // *        Find last unreduced submatrix's top index I0, find QMAX and`
` 542:                  // *        EMIN. Find Gershgorin-type bound if Q's much greater than E's.`
` 543:                  // *`
` 544:                  EMAX = ZERO;`
` 545:                  if (N0 > I0)`
` 546:                  {`
` 547:                      EMIN = Math.Abs(Z[4 * N0 - 5 + o_z]);`
` 548:                  }`
` 549:                  else`
` 550:                  {`
` 551:                      EMIN = ZERO;`
` 552:                  }`
` 553:                  QMIN = Z[4 * N0 - 3 + o_z];`
` 554:                  QMAX = QMIN;`
` 555:                  for (I4 = 4 * N0; I4 >= 8; I4 +=  - 4)`
` 556:                  {`
` 557:                      if (Z[I4 - 5 + o_z] <= ZERO) goto LABEL100;`
` 558:                      if (QMIN >= FOUR * EMAX)`
` 559:                      {`
` 560:                          QMIN = Math.Min(QMIN, Z[I4 - 3 + o_z]);`
` 561:                          EMAX = Math.Max(EMAX, Z[I4 - 5 + o_z]);`
` 562:                      }`
` 563:                      QMAX = Math.Max(QMAX, Z[I4 - 7 + o_z] + Z[I4 - 5 + o_z]);`
` 564:                      EMIN = Math.Min(EMIN, Z[I4 - 5 + o_z]);`
` 565:                  }`
` 566:                  I4 = 4;`
` 567:                  // *`
` 568:              LABEL100:;`
` 569:                  I0 = I4 / 4;`
` 570:                  // *`
` 571:                  // *        Store EMIN for passing to DLAZQ3.`
` 572:                  // *`
` 573:                  Z[4 * N0 - 1 + o_z] = EMIN;`
` 574:                  // *`
` 575:                  // *        Put -(initial shift) into DMIN.`
` 576:                  // *`
` 577:                  DMIN =  - Math.Max(ZERO, QMIN - TWO * Math.Sqrt(QMIN) * Math.Sqrt(EMAX));`
` 578:                  // *`
` 579:                  // *        Now I0:N0 is unreduced. PP = 0 for ping, PP = 1 for pong.`
` 580:                  // *`
` 581:                  PP = 0;`
` 582:                  // *`
` 583:                  NBIG = 30 * (N0 - I0 + 1);`
` 584:                  for (IWHILB = 1; IWHILB <= NBIG; IWHILB++)`
` 585:                  {`
` 586:                      if (I0 > N0) goto LABEL130;`
` 587:                      // *`
` 588:                      // *           While submatrix unfinished take a good dqds step.`
` 589:                      // *`
` 590:                      this._dlazq3.Run(I0, ref N0, ref Z, offset_z, PP, ref DMIN, ref SIGMA`
` 591:                                       , ref DESIG, ref QMAX, ref NFAIL, ref ITER, ref NDIV, IEEE`
` 592:                                       , ref TTYPE, ref DMIN1, ref DMIN2, ref DN, ref DN1, ref DN2`
` 593:                                       , ref TAU);`
` 594:                      // *`
` 595:                      PP = 1 - PP;`
` 596:                      // *`
` 597:                      // *           When EMIN is very small check for splits.`
` 598:                      // *`
` 599:                      if (PP == 0 && N0 - I0 >= 3)`
` 600:                      {`
` 601:                          if (Z[4 * N0 + o_z] <= TOL2 * QMAX || Z[4 * N0 - 1 + o_z] <= TOL2 * SIGMA)`
` 602:                          {`
` 603:                              SPLT = I0 - 1;`
` 604:                              QMAX = Z[4 * I0 - 3 + o_z];`
` 605:                              EMIN = Z[4 * I0 - 1 + o_z];`
` 606:                              OLDEMN = Z[4 * I0 + o_z];`
` 607:                              for (I4 = 4 * I0; I4 <= 4 * (N0 - 3); I4 += 4)`
` 608:                              {`
` 609:                                  if (Z[I4 + o_z] <= TOL2 * Z[I4 - 3 + o_z] || Z[I4 - 1 + o_z] <= TOL2 * SIGMA)`
` 610:                                  {`
` 611:                                      Z[I4 - 1 + o_z] =  - SIGMA;`
` 612:                                      SPLT = I4 / 4;`
` 613:                                      QMAX = ZERO;`
` 614:                                      EMIN = Z[I4 + 3 + o_z];`
` 615:                                      OLDEMN = Z[I4 + 4 + o_z];`
` 616:                                  }`
` 617:                                  else`
` 618:                                  {`
` 619:                                      QMAX = Math.Max(QMAX, Z[I4 + 1 + o_z]);`
` 620:                                      EMIN = Math.Min(EMIN, Z[I4 - 1 + o_z]);`
` 621:                                      OLDEMN = Math.Min(OLDEMN, Z[I4 + o_z]);`
` 622:                                  }`
` 623:                              }`
` 624:                              Z[4 * N0 - 1 + o_z] = EMIN;`
` 625:                              Z[4 * N0 + o_z] = OLDEMN;`
` 626:                              I0 = SPLT + 1;`
` 627:                          }`
` 628:                      }`
` 629:                      // *`
` 630:                  }`
` 631:                  // *`
` 632:                  INFO = 2;`
` 633:                  return;`
` 634:                  // *`
` 635:                  // *        end IWHILB`
` 636:                  // *`
` 637:              LABEL130:;`
` 638:                  // *`
` 639:              }`
` 640:              // *`
` 641:              INFO = 3;`
` 642:              return;`
` 643:              // *`
` 644:              // *     end IWHILA   `
` 645:              // *`
` 646:          LABEL150:;`
` 647:              // *      `
` 648:              // *     Move q's to the front.`
` 649:              // *      `
` 650:              for (K = 2; K <= N; K++)`
` 651:              {`
` 652:                  Z[K + o_z] = Z[4 * K - 3 + o_z];`
` 653:              }`
` 654:              // *      `
` 655:              // *     Sort and compute sum of eigenvalues.`
` 656:              // *`
` 657:              this._dlasrt.Run("D", N, ref Z, offset_z, ref IINFO);`
` 658:              // *`
` 659:              E = ZERO;`
` 660:              for (K = N; K >= 1; K +=  - 1)`
` 661:              {`
` 662:                  E = E + Z[K + o_z];`
` 663:              }`
` 664:              // *`
` 665:              // *     Store trace, sum(eigenvalues) and information on performance.`
` 666:              // *`
` 667:              Z[2 * N + 1 + o_z] = TRACE;`
` 668:              Z[2 * N + 2 + o_z] = E;`
` 669:              Z[2 * N + 3 + o_z] = Convert.ToDouble(ITER);`
` 670:              Z[2 * N + 4 + o_z] = Convert.ToDouble(NDIV) / Convert.ToDouble(Math.Pow(N,2));`
` 671:              Z[2 * N + 5 + o_z] = HUNDRD * NFAIL / Convert.ToDouble(ITER);`
` 672:              return;`
` 673:              // *`
` 674:              // *     End of DLASQ2`
` 675:              // *`
` 676:   `
` 677:              #endregion`
` 678:   `
` 679:          }`
` 680:      }`
` 681:  }`