`   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:      /// This subroutine computes the I-th updated eigenvalue of a symmetric`
`  27:      /// rank-one modification to a diagonal matrix whose elements are`
`  28:      /// given in the array d, and that`
`  29:      /// `
`  30:      /// D(i) .LT. D(j)  for  i .LT. j`
`  31:      /// `
`  32:      /// and that RHO .GT. 0.  This is arranged by the calling routine, and is`
`  33:      /// no loss in generality.  The rank-one modified system is thus`
`  34:      /// `
`  35:      /// diag( D )  +  RHO *  Z * Z_transpose.`
`  36:      /// `
`  37:      /// where we assume the Euclidean norm of Z is 1.`
`  38:      /// `
`  39:      /// The method consists of approximating the rational functions in the`
`  40:      /// secular equation by simpler interpolating rational functions.`
`  41:      /// `
`  42:      ///</summary>`
`  43:      public class DLAED4`
`  44:      {`
`  45:      `
`  46:   `
`  47:          #region Dependencies`
`  48:          `
`  49:          DLAMCH _dlamch; DLAED5 _dlaed5; DLAED6 _dlaed6; `
`  50:   `
`  51:          #endregion`
`  52:   `
`  53:   `
`  54:          #region Fields`
`  55:          `
`  56:          const int MAXIT = 30; const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; `
`  57:          const double THREE = 3.0E0;const double FOUR = 4.0E0; const double EIGHT = 8.0E0; const double TEN = 10.0E0; `
`  58:          bool ORGATI = false;bool SWTCH = false; bool SWTCH3 = false; int II = 0; int IIM1 = 0; int IIP1 = 0; int IP1 = 0; `
`  59:          int ITER = 0;int J = 0; int NITER = 0; double A = 0; double B = 0; double C = 0; double DEL = 0; double DLTLB = 0; `
`  60:          double DLTUB = 0;double DPHI = 0; double DPSI = 0; double DW = 0; double EPS = 0; double ERRETM = 0; double ETA = 0; `
`  61:          double MIDPT = 0;double PHI = 0; double PREW = 0; double PSI = 0; double RHOINV = 0; double TAU = 0; double TEMP = 0; `
`  62:          double TEMP1 = 0;double W = 0; double[] ZZ = new double[3]; int offset_zz = 0; int o_zz = -1; `
`  63:   `
`  64:          #endregion`
`  65:   `
`  66:          public DLAED4(DLAMCH dlamch, DLAED5 dlaed5, DLAED6 dlaed6)`
`  67:          {`
`  68:      `
`  69:   `
`  70:              #region Set Dependencies`
`  71:              `
`  72:              this._dlamch = dlamch; this._dlaed5 = dlaed5; this._dlaed6 = dlaed6; `
`  73:   `
`  74:              #endregion`
`  75:   `
`  76:          }`
`  77:      `
`  78:          public DLAED4()`
`  79:          {`
`  80:      `
`  81:   `
`  82:              #region Dependencies (Initialization)`
`  83:              `
`  84:              LSAME lsame = new LSAME();`
`  85:              DLAMC3 dlamc3 = new DLAMC3();`
`  86:              DLAED5 dlaed5 = new DLAED5();`
`  87:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
`  88:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
`  89:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
`  90:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
`  91:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
`  92:              DLAED6 dlaed6 = new DLAED6(dlamch);`
`  93:   `
`  94:              #endregion`
`  95:   `
`  96:   `
`  97:              #region Set Dependencies`
`  98:              `
`  99:              this._dlamch = dlamch; this._dlaed5 = dlaed5; this._dlaed6 = dlaed6; `
` 100:   `
` 101:              #endregion`
` 102:   `
` 103:          }`
` 104:          /// <summary>`
` 105:          /// Purpose`
` 106:          /// =======`
` 107:          /// `
` 108:          /// This subroutine computes the I-th updated eigenvalue of a symmetric`
` 109:          /// rank-one modification to a diagonal matrix whose elements are`
` 110:          /// given in the array d, and that`
` 111:          /// `
` 112:          /// D(i) .LT. D(j)  for  i .LT. j`
` 113:          /// `
` 114:          /// and that RHO .GT. 0.  This is arranged by the calling routine, and is`
` 115:          /// no loss in generality.  The rank-one modified system is thus`
` 116:          /// `
` 117:          /// diag( D )  +  RHO *  Z * Z_transpose.`
` 118:          /// `
` 119:          /// where we assume the Euclidean norm of Z is 1.`
` 120:          /// `
` 121:          /// The method consists of approximating the rational functions in the`
` 122:          /// secular equation by simpler interpolating rational functions.`
` 123:          /// `
` 124:          ///</summary>`
` 125:          /// <param name="N">`
` 126:          /// (input) INTEGER`
` 127:          /// The length of all arrays.`
` 128:          ///</param>`
` 129:          /// <param name="I">`
` 130:          /// (input) INTEGER`
` 131:          /// The index of the eigenvalue to be computed.  1 .LE. I .LE. N.`
` 132:          ///</param>`
` 133:          /// <param name="D">`
` 134:          /// (input) DOUBLE PRECISION array, dimension (N)`
` 135:          /// The original eigenvalues.  It is assumed that they are in`
` 136:          /// order, D(I) .LT. D(J)  for I .LT. J.`
` 137:          ///</param>`
` 138:          /// <param name="Z">`
` 139:          /// (input) DOUBLE PRECISION array, dimension (N)`
` 140:          /// The components of the updating vector.`
` 141:          ///</param>`
` 142:          /// <param name="DELTA">`
` 143:          /// (output) DOUBLE PRECISION array, dimension (N)`
` 144:          /// If N .GT. 2, DELTA contains (D(j) - lambda_I) in its  j-th`
` 145:          /// component.  If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5`
` 146:          /// for detail. The vector DELTA contains the information necessary`
` 147:          /// to construct the eigenvectors by DLAED3 and DLAED9.`
` 148:          ///</param>`
` 149:          /// <param name="RHO">`
` 150:          /// (input) DOUBLE PRECISION`
` 151:          /// The scalar in the symmetric updating formula.`
` 152:          ///</param>`
` 153:          /// <param name="DLAM">`
` 154:          /// (output) DOUBLE PRECISION`
` 155:          /// The computed lambda_I, the I-th updated eigenvalue.`
` 156:          ///</param>`
` 157:          /// <param name="INFO">`
` 158:          /// (output) INTEGER`
` 159:          /// = 0:  successful exit`
` 160:          /// .GT. 0:  if INFO = 1, the updating process failed.`
` 161:          ///</param>`
` 162:          public void Run(int N, int I, double[] D, int offset_d, double[] Z, int offset_z, ref double[] DELTA, int offset_delta, double RHO`
` 163:                           , ref double DLAM, ref int INFO)`
` 164:          {`
` 165:   `
` 166:              #region Array Index Correction`
` 167:              `
` 168:               int o_d = -1 + offset_d;  int o_z = -1 + offset_z;  int o_delta = -1 + offset_delta; `
` 169:   `
` 170:              #endregion`
` 171:   `
` 172:   `
` 173:              #region Prolog`
` 174:              `
` 175:              // *`
` 176:              // *  -- LAPACK routine (version 3.1) --`
` 177:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 178:              // *     November 2006`
` 179:              // *`
` 180:              // *     .. Scalar Arguments ..`
` 181:              // *     ..`
` 182:              // *     .. Array Arguments ..`
` 183:              // *     ..`
` 184:              // *`
` 185:              // *  Purpose`
` 186:              // *  =======`
` 187:              // *`
` 188:              // *  This subroutine computes the I-th updated eigenvalue of a symmetric`
` 189:              // *  rank-one modification to a diagonal matrix whose elements are`
` 190:              // *  given in the array d, and that`
` 191:              // *`
` 192:              // *             D(i) < D(j)  for  i < j`
` 193:              // *`
` 194:              // *  and that RHO > 0.  This is arranged by the calling routine, and is`
` 195:              // *  no loss in generality.  The rank-one modified system is thus`
` 196:              // *`
` 197:              // *             diag( D )  +  RHO *  Z * Z_transpose.`
` 198:              // *`
` 199:              // *  where we assume the Euclidean norm of Z is 1.`
` 200:              // *`
` 201:              // *  The method consists of approximating the rational functions in the`
` 202:              // *  secular equation by simpler interpolating rational functions.`
` 203:              // *`
` 204:              // *  Arguments`
` 205:              // *  =========`
` 206:              // *`
` 207:              // *  N      (input) INTEGER`
` 208:              // *         The length of all arrays.`
` 209:              // *`
` 210:              // *  I      (input) INTEGER`
` 211:              // *         The index of the eigenvalue to be computed.  1 <= I <= N.`
` 212:              // *`
` 213:              // *  D      (input) DOUBLE PRECISION array, dimension (N)`
` 214:              // *         The original eigenvalues.  It is assumed that they are in`
` 215:              // *         order, D(I) < D(J)  for I < J.`
` 216:              // *`
` 217:              // *  Z      (input) DOUBLE PRECISION array, dimension (N)`
` 218:              // *         The components of the updating vector.`
` 219:              // *`
` 220:              // *  DELTA  (output) DOUBLE PRECISION array, dimension (N)`
` 221:              // *         If N .GT. 2, DELTA contains (D(j) - lambda_I) in its  j-th`
` 222:              // *         component.  If N = 1, then DELTA(1) = 1. If N = 2, see DLAED5`
` 223:              // *         for detail. The vector DELTA contains the information necessary`
` 224:              // *         to construct the eigenvectors by DLAED3 and DLAED9.`
` 225:              // *`
` 226:              // *  RHO    (input) DOUBLE PRECISION`
` 227:              // *         The scalar in the symmetric updating formula.`
` 228:              // *`
` 229:              // *  DLAM   (output) DOUBLE PRECISION`
` 230:              // *         The computed lambda_I, the I-th updated eigenvalue.`
` 231:              // *`
` 232:              // *  INFO   (output) INTEGER`
` 233:              // *         = 0:  successful exit`
` 234:              // *         > 0:  if INFO = 1, the updating process failed.`
` 235:              // *`
` 236:              // *  Internal Parameters`
` 237:              // *  ===================`
` 238:              // *`
` 239:              // *  Logical variable ORGATI (origin-at-i?) is used for distinguishing`
` 240:              // *  whether D(i) or D(i+1) is treated as the origin.`
` 241:              // *`
` 242:              // *            ORGATI = .true.    origin at i`
` 243:              // *            ORGATI = .false.   origin at i+1`
` 244:              // *`
` 245:              // *   Logical variable SWTCH3 (switch-for-3-poles?) is for noting`
` 246:              // *   if we are working with THREE poles!`
` 247:              // *`
` 248:              // *   MAXIT is the maximum number of iterations allowed for each`
` 249:              // *   eigenvalue.`
` 250:              // *`
` 251:              // *  Further Details`
` 252:              // *  ===============`
` 253:              // *`
` 254:              // *  Based on contributions by`
` 255:              // *     Ren-Cang Li, Computer Science Division, University of California`
` 256:              // *     at Berkeley, USA`
` 257:              // *`
` 258:              // *  =====================================================================`
` 259:              // *`
` 260:              // *     .. Parameters ..`
` 261:              // *     ..`
` 262:              // *     .. Local Scalars ..`
` 263:              // *     ..`
` 264:              // *     .. Local Arrays ..`
` 265:              // *     ..`
` 266:              // *     .. External Functions ..`
` 267:              // *     ..`
` 268:              // *     .. External Subroutines ..`
` 269:              // *     ..`
` 270:              // *     .. Intrinsic Functions ..`
` 271:              //      INTRINSIC          ABS, MAX, MIN, SQRT;`
` 272:              // *     ..`
` 273:              // *     .. Executable Statements ..`
` 274:              // *`
` 275:              // *     Since this routine is called in an inner loop, we do no argument`
` 276:              // *     checking.`
` 277:              // *`
` 278:              // *     Quick return for N=1 and 2.`
` 279:              // *`
` 280:   `
` 281:              #endregion`
` 282:   `
` 283:   `
` 284:              #region Body`
` 285:              `
` 286:              INFO = 0;`
` 287:              if (N == 1)`
` 288:              {`
` 289:                  // *`
` 290:                  // *         Presumably, I=1 upon entry`
` 291:                  // *`
` 292:                  DLAM = D[1 + o_d] + RHO * Z[1 + o_z] * Z[1 + o_z];`
` 293:                  DELTA[1 + o_delta] = ONE;`
` 294:                  return;`
` 295:              }`
` 296:              if (N == 2)`
` 297:              {`
` 298:                  this._dlaed5.Run(I, D, offset_d, Z, offset_z, ref DELTA, offset_delta, RHO, ref DLAM);`
` 299:                  return;`
` 300:              }`
` 301:              // *`
` 302:              // *     Compute machine epsilon`
` 303:              // *`
` 304:              EPS = this._dlamch.Run("Epsilon");`
` 305:              RHOINV = ONE / RHO;`
` 306:              // *`
` 307:              // *     The case I = N`
` 308:              // *`
` 309:              if (I == N)`
` 310:              {`
` 311:                  // *`
` 312:                  // *        Initialize some basic variables`
` 313:                  // *`
` 314:                  II = N - 1;`
` 315:                  NITER = 1;`
` 316:                  // *`
` 317:                  // *        Calculate initial guess`
` 318:                  // *`
` 319:                  MIDPT = RHO / TWO;`
` 320:                  // *`
` 321:                  // *        If ||Z||_2 is not one, then TEMP should be set to`
` 322:                  // *        RHO * ||Z||_2^2 / TWO`
` 323:                  // *`
` 324:                  for (J = 1; J <= N; J++)`
` 325:                  {`
` 326:                      DELTA[J + o_delta] = (D[J + o_d] - D[I + o_d]) - MIDPT;`
` 327:                  }`
` 328:                  // *`
` 329:                  PSI = ZERO;`
` 330:                  for (J = 1; J <= N - 2; J++)`
` 331:                  {`
` 332:                      PSI = PSI + Z[J + o_z] * Z[J + o_z] / DELTA[J + o_delta];`
` 333:                  }`
` 334:                  // *`
` 335:                  C = RHOINV + PSI;`
` 336:                  W = C + Z[II + o_z] * Z[II + o_z] / DELTA[II + o_delta] + Z[N + o_z] * Z[N + o_z] / DELTA[N + o_delta];`
` 337:                  // *`
` 338:                  if (W <= ZERO)`
` 339:                  {`
` 340:                      TEMP = Z[N - 1 + o_z] * Z[N - 1 + o_z] / (D[N + o_d] - D[N - 1 + o_d] + RHO) + Z[N + o_z] * Z[N + o_z] / RHO;`
` 341:                      if (C <= TEMP)`
` 342:                      {`
` 343:                          TAU = RHO;`
` 344:                      }`
` 345:                      else`
` 346:                      {`
` 347:                          DEL = D[N + o_d] - D[N - 1 + o_d];`
` 348:                          A =  - C * DEL + Z[N - 1 + o_z] * Z[N - 1 + o_z] + Z[N + o_z] * Z[N + o_z];`
` 349:                          B = Z[N + o_z] * Z[N + o_z] * DEL;`
` 350:                          if (A < ZERO)`
` 351:                          {`
` 352:                              TAU = TWO * B / (Math.Sqrt(A * A + FOUR * B * C) - A);`
` 353:                          }`
` 354:                          else`
` 355:                          {`
` 356:                              TAU = (A + Math.Sqrt(A * A + FOUR * B * C)) / (TWO * C);`
` 357:                          }`
` 358:                      }`
` 359:                      // *`
` 360:                      // *           It can be proved that`
` 361:                      // *               D(N)+RHO/2 <= LAMBDA(N) < D(N)+TAU <= D(N)+RHO`
` 362:                      // *`
` 363:                      DLTLB = MIDPT;`
` 364:                      DLTUB = RHO;`
` 365:                  }`
` 366:                  else`
` 367:                  {`
` 368:                      DEL = D[N + o_d] - D[N - 1 + o_d];`
` 369:                      A =  - C * DEL + Z[N - 1 + o_z] * Z[N - 1 + o_z] + Z[N + o_z] * Z[N + o_z];`
` 370:                      B = Z[N + o_z] * Z[N + o_z] * DEL;`
` 371:                      if (A < ZERO)`
` 372:                      {`
` 373:                          TAU = TWO * B / (Math.Sqrt(A * A + FOUR * B * C) - A);`
` 374:                      }`
` 375:                      else`
` 376:                      {`
` 377:                          TAU = (A + Math.Sqrt(A * A + FOUR * B * C)) / (TWO * C);`
` 378:                      }`
` 379:                      // *`
` 380:                      // *           It can be proved that`
` 381:                      // *               D(N) < D(N)+TAU < LAMBDA(N) < D(N)+RHO/2`
` 382:                      // *`
` 383:                      DLTLB = ZERO;`
` 384:                      DLTUB = MIDPT;`
` 385:                  }`
` 386:                  // *`
` 387:                  for (J = 1; J <= N; J++)`
` 388:                  {`
` 389:                      DELTA[J + o_delta] = (D[J + o_d] - D[I + o_d]) - TAU;`
` 390:                  }`
` 391:                  // *`
` 392:                  // *        Evaluate PSI and the derivative DPSI`
` 393:                  // *`
` 394:                  DPSI = ZERO;`
` 395:                  PSI = ZERO;`
` 396:                  ERRETM = ZERO;`
` 397:                  for (J = 1; J <= II; J++)`
` 398:                  {`
` 399:                      TEMP = Z[J + o_z] / DELTA[J + o_delta];`
` 400:                      PSI = PSI + Z[J + o_z] * TEMP;`
` 401:                      DPSI = DPSI + TEMP * TEMP;`
` 402:                      ERRETM = ERRETM + PSI;`
` 403:                  }`
` 404:                  ERRETM = Math.Abs(ERRETM);`
` 405:                  // *`
` 406:                  // *        Evaluate PHI and the derivative DPHI`
` 407:                  // *`
` 408:                  TEMP = Z[N + o_z] / DELTA[N + o_delta];`
` 409:                  PHI = Z[N + o_z] * TEMP;`
` 410:                  DPHI = TEMP * TEMP;`
` 411:                  ERRETM = EIGHT * ( - PHI - PSI) + ERRETM - PHI + RHOINV + Math.Abs(TAU) * (DPSI + DPHI);`
` 412:                  // *`
` 413:                  W = RHOINV + PHI + PSI;`
` 414:                  // *`
` 415:                  // *        Test for convergence`
` 416:                  // *`
` 417:                  if (Math.Abs(W) <= EPS * ERRETM)`
` 418:                  {`
` 419:                      DLAM = D[I + o_d] + TAU;`
` 420:                      goto LABEL250;`
` 421:                  }`
` 422:                  // *`
` 423:                  if (W <= ZERO)`
` 424:                  {`
` 425:                      DLTLB = Math.Max(DLTLB, TAU);`
` 426:                  }`
` 427:                  else`
` 428:                  {`
` 429:                      DLTUB = Math.Min(DLTUB, TAU);`
` 430:                  }`
` 431:                  // *`
` 432:                  // *        Calculate the new step`
` 433:                  // *`
` 434:                  NITER = NITER + 1;`
` 435:                  C = W - DELTA[N - 1 + o_delta] * DPSI - DELTA[N + o_delta] * DPHI;`
` 436:                  A = (DELTA[N - 1 + o_delta] + DELTA[N + o_delta]) * W - DELTA[N - 1 + o_delta] * DELTA[N + o_delta] * (DPSI + DPHI);`
` 437:                  B = DELTA[N - 1 + o_delta] * DELTA[N + o_delta] * W;`
` 438:                  if (C < ZERO) C = Math.Abs(C);`
` 439:                  if (C == ZERO)`
` 440:                  {`
` 441:                      // *          ETA = B/A`
` 442:                      // *           ETA = RHO - TAU`
` 443:                      ETA = DLTUB - TAU;`
` 444:                  }`
` 445:                  else`
` 446:                  {`
` 447:                      if (A >= ZERO)`
` 448:                      {`
` 449:                          ETA = (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);`
` 450:                      }`
` 451:                      else`
` 452:                      {`
` 453:                          ETA = TWO * B / (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));`
` 454:                      }`
` 455:                  }`
` 456:                  // *`
` 457:                  // *        Note, eta should be positive if w is negative, and`
` 458:                  // *        eta should be negative otherwise. However,`
` 459:                  // *        if for some reason caused by roundoff, eta*w > 0,`
` 460:                  // *        we simply use one Newton step instead. This way`
` 461:                  // *        will guarantee eta*w < 0.`
` 462:                  // *`
` 463:                  if (W * ETA > ZERO) ETA =  - W / (DPSI + DPHI);`
` 464:                  TEMP = TAU + ETA;`
` 465:                  if (TEMP > DLTUB || TEMP < DLTLB)`
` 466:                  {`
` 467:                      if (W < ZERO)`
` 468:                      {`
` 469:                          ETA = (DLTUB - TAU) / TWO;`
` 470:                      }`
` 471:                      else`
` 472:                      {`
` 473:                          ETA = (DLTLB - TAU) / TWO;`
` 474:                      }`
` 475:                  }`
` 476:                  for (J = 1; J <= N; J++)`
` 477:                  {`
` 478:                      DELTA[J + o_delta] = DELTA[J + o_delta] - ETA;`
` 479:                  }`
` 480:                  // *`
` 481:                  TAU = TAU + ETA;`
` 482:                  // *`
` 483:                  // *        Evaluate PSI and the derivative DPSI`
` 484:                  // *`
` 485:                  DPSI = ZERO;`
` 486:                  PSI = ZERO;`
` 487:                  ERRETM = ZERO;`
` 488:                  for (J = 1; J <= II; J++)`
` 489:                  {`
` 490:                      TEMP = Z[J + o_z] / DELTA[J + o_delta];`
` 491:                      PSI = PSI + Z[J + o_z] * TEMP;`
` 492:                      DPSI = DPSI + TEMP * TEMP;`
` 493:                      ERRETM = ERRETM + PSI;`
` 494:                  }`
` 495:                  ERRETM = Math.Abs(ERRETM);`
` 496:                  // *`
` 497:                  // *        Evaluate PHI and the derivative DPHI`
` 498:                  // *`
` 499:                  TEMP = Z[N + o_z] / DELTA[N + o_delta];`
` 500:                  PHI = Z[N + o_z] * TEMP;`
` 501:                  DPHI = TEMP * TEMP;`
` 502:                  ERRETM = EIGHT * ( - PHI - PSI) + ERRETM - PHI + RHOINV + Math.Abs(TAU) * (DPSI + DPHI);`
` 503:                  // *`
` 504:                  W = RHOINV + PHI + PSI;`
` 505:                  // *`
` 506:                  // *        Main loop to update the values of the array   DELTA`
` 507:                  // *`
` 508:                  ITER = NITER + 1;`
` 509:                  // *`
` 510:                  for (NITER = ITER; NITER <= MAXIT; NITER++)`
` 511:                  {`
` 512:                      // *`
` 513:                      // *           Test for convergence`
` 514:                      // *`
` 515:                      if (Math.Abs(W) <= EPS * ERRETM)`
` 516:                      {`
` 517:                          DLAM = D[I + o_d] + TAU;`
` 518:                          goto LABEL250;`
` 519:                      }`
` 520:                      // *`
` 521:                      if (W <= ZERO)`
` 522:                      {`
` 523:                          DLTLB = Math.Max(DLTLB, TAU);`
` 524:                      }`
` 525:                      else`
` 526:                      {`
` 527:                          DLTUB = Math.Min(DLTUB, TAU);`
` 528:                      }`
` 529:                      // *`
` 530:                      // *           Calculate the new step`
` 531:                      // *`
` 532:                      C = W - DELTA[N - 1 + o_delta] * DPSI - DELTA[N + o_delta] * DPHI;`
` 533:                      A = (DELTA[N - 1 + o_delta] + DELTA[N + o_delta]) * W - DELTA[N - 1 + o_delta] * DELTA[N + o_delta] * (DPSI + DPHI);`
` 534:                      B = DELTA[N - 1 + o_delta] * DELTA[N + o_delta] * W;`
` 535:                      if (A >= ZERO)`
` 536:                      {`
` 537:                          ETA = (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);`
` 538:                      }`
` 539:                      else`
` 540:                      {`
` 541:                          ETA = TWO * B / (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));`
` 542:                      }`
` 543:                      // *`
` 544:                      // *           Note, eta should be positive if w is negative, and`
` 545:                      // *           eta should be negative otherwise. However,`
` 546:                      // *           if for some reason caused by roundoff, eta*w > 0,`
` 547:                      // *           we simply use one Newton step instead. This way`
` 548:                      // *           will guarantee eta*w < 0.`
` 549:                      // *`
` 550:                      if (W * ETA > ZERO) ETA =  - W / (DPSI + DPHI);`
` 551:                      TEMP = TAU + ETA;`
` 552:                      if (TEMP > DLTUB || TEMP < DLTLB)`
` 553:                      {`
` 554:                          if (W < ZERO)`
` 555:                          {`
` 556:                              ETA = (DLTUB - TAU) / TWO;`
` 557:                          }`
` 558:                          else`
` 559:                          {`
` 560:                              ETA = (DLTLB - TAU) / TWO;`
` 561:                          }`
` 562:                      }`
` 563:                      for (J = 1; J <= N; J++)`
` 564:                      {`
` 565:                          DELTA[J + o_delta] = DELTA[J + o_delta] - ETA;`
` 566:                      }`
` 567:                      // *`
` 568:                      TAU = TAU + ETA;`
` 569:                      // *`
` 570:                      // *           Evaluate PSI and the derivative DPSI`
` 571:                      // *`
` 572:                      DPSI = ZERO;`
` 573:                      PSI = ZERO;`
` 574:                      ERRETM = ZERO;`
` 575:                      for (J = 1; J <= II; J++)`
` 576:                      {`
` 577:                          TEMP = Z[J + o_z] / DELTA[J + o_delta];`
` 578:                          PSI = PSI + Z[J + o_z] * TEMP;`
` 579:                          DPSI = DPSI + TEMP * TEMP;`
` 580:                          ERRETM = ERRETM + PSI;`
` 581:                      }`
` 582:                      ERRETM = Math.Abs(ERRETM);`
` 583:                      // *`
` 584:                      // *           Evaluate PHI and the derivative DPHI`
` 585:                      // *`
` 586:                      TEMP = Z[N + o_z] / DELTA[N + o_delta];`
` 587:                      PHI = Z[N + o_z] * TEMP;`
` 588:                      DPHI = TEMP * TEMP;`
` 589:                      ERRETM = EIGHT * ( - PHI - PSI) + ERRETM - PHI + RHOINV + Math.Abs(TAU) * (DPSI + DPHI);`
` 590:                      // *`
` 591:                      W = RHOINV + PHI + PSI;`
` 592:                  }`
` 593:                  // *`
` 594:                  // *        Return with INFO = 1, NITER = MAXIT and not converged`
` 595:                  // *`
` 596:                  INFO = 1;`
` 597:                  DLAM = D[I + o_d] + TAU;`
` 598:                  goto LABEL250;`
` 599:                  // *`
` 600:                  // *        End for the case I = N`
` 601:                  // *`
` 602:              }`
` 603:              else`
` 604:              {`
` 605:                  // *`
` 606:                  // *        The case for I < N`
` 607:                  // *`
` 608:                  NITER = 1;`
` 609:                  IP1 = I + 1;`
` 610:                  // *`
` 611:                  // *        Calculate initial guess`
` 612:                  // *`
` 613:                  DEL = D[IP1 + o_d] - D[I + o_d];`
` 614:                  MIDPT = DEL / TWO;`
` 615:                  for (J = 1; J <= N; J++)`
` 616:                  {`
` 617:                      DELTA[J + o_delta] = (D[J + o_d] - D[I + o_d]) - MIDPT;`
` 618:                  }`
` 619:                  // *`
` 620:                  PSI = ZERO;`
` 621:                  for (J = 1; J <= I - 1; J++)`
` 622:                  {`
` 623:                      PSI = PSI + Z[J + o_z] * Z[J + o_z] / DELTA[J + o_delta];`
` 624:                  }`
` 625:                  // *`
` 626:                  PHI = ZERO;`
` 627:                  for (J = N; J >= I + 2; J +=  - 1)`
` 628:                  {`
` 629:                      PHI = PHI + Z[J + o_z] * Z[J + o_z] / DELTA[J + o_delta];`
` 630:                  }`
` 631:                  C = RHOINV + PSI + PHI;`
` 632:                  W = C + Z[I + o_z] * Z[I + o_z] / DELTA[I + o_delta] + Z[IP1 + o_z] * Z[IP1 + o_z] / DELTA[IP1 + o_delta];`
` 633:                  // *`
` 634:                  if (W > ZERO)`
` 635:                  {`
` 636:                      // *`
` 637:                      // *           d(i)< the ith eigenvalue < (d(i)+d(i+1))/2`
` 638:                      // *`
` 639:                      // *           We choose d(i) as origin.`
` 640:                      // *`
` 641:                      ORGATI = true;`
` 642:                      A = C * DEL + Z[I + o_z] * Z[I + o_z] + Z[IP1 + o_z] * Z[IP1 + o_z];`
` 643:                      B = Z[I + o_z] * Z[I + o_z] * DEL;`
` 644:                      if (A > ZERO)`
` 645:                      {`
` 646:                          TAU = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));`
` 647:                      }`
` 648:                      else`
` 649:                      {`
` 650:                          TAU = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);`
` 651:                      }`
` 652:                      DLTLB = ZERO;`
` 653:                      DLTUB = MIDPT;`
` 654:                  }`
` 655:                  else`
` 656:                  {`
` 657:                      // *`
` 658:                      // *           (d(i)+d(i+1))/2 <= the ith eigenvalue < d(i+1)`
` 659:                      // *`
` 660:                      // *           We choose d(i+1) as origin.`
` 661:                      // *`
` 662:                      ORGATI = false;`
` 663:                      A = C * DEL - Z[I + o_z] * Z[I + o_z] - Z[IP1 + o_z] * Z[IP1 + o_z];`
` 664:                      B = Z[IP1 + o_z] * Z[IP1 + o_z] * DEL;`
` 665:                      if (A < ZERO)`
` 666:                      {`
` 667:                          TAU = TWO * B / (A - Math.Sqrt(Math.Abs(A * A + FOUR * B * C)));`
` 668:                      }`
` 669:                      else`
` 670:                      {`
` 671:                          TAU =  - (A + Math.Sqrt(Math.Abs(A * A + FOUR * B * C))) / (TWO * C);`
` 672:                      }`
` 673:                      DLTLB =  - MIDPT;`
` 674:                      DLTUB = ZERO;`
` 675:                  }`
` 676:                  // *`
` 677:                  if (ORGATI)`
` 678:                  {`
` 679:                      for (J = 1; J <= N; J++)`
` 680:                      {`
` 681:                          DELTA[J + o_delta] = (D[J + o_d] - D[I + o_d]) - TAU;`
` 682:                      }`
` 683:                  }`
` 684:                  else`
` 685:                  {`
` 686:                      for (J = 1; J <= N; J++)`
` 687:                      {`
` 688:                          DELTA[J + o_delta] = (D[J + o_d] - D[IP1 + o_d]) - TAU;`
` 689:                      }`
` 690:                  }`
` 691:                  if (ORGATI)`
` 692:                  {`
` 693:                      II = I;`
` 694:                  }`
` 695:                  else`
` 696:                  {`
` 697:                      II = I + 1;`
` 698:                  }`
` 699:                  IIM1 = II - 1;`
` 700:                  IIP1 = II + 1;`
` 701:                  // *`
` 702:                  // *        Evaluate PSI and the derivative DPSI`
` 703:                  // *`
` 704:                  DPSI = ZERO;`
` 705:                  PSI = ZERO;`
` 706:                  ERRETM = ZERO;`
` 707:                  for (J = 1; J <= IIM1; J++)`
` 708:                  {`
` 709:                      TEMP = Z[J + o_z] / DELTA[J + o_delta];`
` 710:                      PSI = PSI + Z[J + o_z] * TEMP;`
` 711:                      DPSI = DPSI + TEMP * TEMP;`
` 712:                      ERRETM = ERRETM + PSI;`
` 713:                  }`
` 714:                  ERRETM = Math.Abs(ERRETM);`
` 715:                  // *`
` 716:                  // *        Evaluate PHI and the derivative DPHI`
` 717:                  // *`
` 718:                  DPHI = ZERO;`
` 719:                  PHI = ZERO;`
` 720:                  for (J = N; J >= IIP1; J +=  - 1)`
` 721:                  {`
` 722:                      TEMP = Z[J + o_z] / DELTA[J + o_delta];`
` 723:                      PHI = PHI + Z[J + o_z] * TEMP;`
` 724:                      DPHI = DPHI + TEMP * TEMP;`
` 725:                      ERRETM = ERRETM + PHI;`
` 726:                  }`
` 727:                  // *`
` 728:                  W = RHOINV + PHI + PSI;`
` 729:                  // *`
` 730:                  // *        W is the value of the secular function with`
` 731:                  // *        its ii-th element removed.`
` 732:                  // *`
` 733:                  SWTCH3 = false;`
` 734:                  if (ORGATI)`
` 735:                  {`
` 736:                      if (W < ZERO) SWTCH3 = true;`
` 737:                  }`
` 738:                  else`
` 739:                  {`
` 740:                      if (W > ZERO) SWTCH3 = true;`
` 741:                  }`
` 742:                  if (II == 1 || II == N) SWTCH3 = false;`
` 743:                  // *`
` 744:                  TEMP = Z[II + o_z] / DELTA[II + o_delta];`
` 745:                  DW = DPSI + DPHI + TEMP * TEMP;`
` 746:                  TEMP = Z[II + o_z] * TEMP;`
` 747:                  W = W + TEMP;`
` 748:                  ERRETM = EIGHT * (PHI - PSI) + ERRETM + TWO * RHOINV + THREE * Math.Abs(TEMP) + Math.Abs(TAU) * DW;`
` 749:                  // *`
` 750:                  // *        Test for convergence`
` 751:                  // *`
` 752:                  if (Math.Abs(W) <= EPS * ERRETM)`
` 753:                  {`
` 754:                      if (ORGATI)`
` 755:                      {`
` 756:                          DLAM = D[I + o_d] + TAU;`
` 757:                      }`
` 758:                      else`
` 759:                      {`
` 760:                          DLAM = D[IP1 + o_d] + TAU;`
` 761:                      }`
` 762:                      goto LABEL250;`
` 763:                  }`
` 764:                  // *`
` 765:                  if (W <= ZERO)`
` 766:                  {`
` 767:                      DLTLB = Math.Max(DLTLB, TAU);`
` 768:                  }`
` 769:                  else`
` 770:                  {`
` 771:                      DLTUB = Math.Min(DLTUB, TAU);`
` 772:                  }`
` 773:                  // *`
` 774:                  // *        Calculate the new step`
` 775:                  // *`
` 776:                  NITER = NITER + 1;`
` 777:                  if (!SWTCH3)`
` 778:                  {`
` 779:                      if (ORGATI)`
` 780:                      {`
` 781:                          C = W - DELTA[IP1 + o_delta] * DW - (D[I + o_d] - D[IP1 + o_d]) * Math.Pow(Z[I + o_z] / DELTA[I + o_delta],2);`
` 782:                      }`
` 783:                      else`
` 784:                      {`
` 785:                          C = W - DELTA[I + o_delta] * DW - (D[IP1 + o_d] - D[I + o_d]) * Math.Pow(Z[IP1 + o_z] / DELTA[IP1 + o_delta],2);`
` 786:                      }`
` 787:                      A = (DELTA[I + o_delta] + DELTA[IP1 + o_delta]) * W - DELTA[I + o_delta] * DELTA[IP1 + o_delta] * DW;`
` 788:                      B = DELTA[I + o_delta] * DELTA[IP1 + o_delta] * W;`
` 789:                      if (C == ZERO)`
` 790:                      {`
` 791:                          if (A == ZERO)`
` 792:                          {`
` 793:                              if (ORGATI)`
` 794:                              {`
` 795:                                  A = Z[I + o_z] * Z[I + o_z] + DELTA[IP1 + o_delta] * DELTA[IP1 + o_delta] * (DPSI + DPHI);`
` 796:                              }`
` 797:                              else`
` 798:                              {`
` 799:                                  A = Z[IP1 + o_z] * Z[IP1 + o_z] + DELTA[I + o_delta] * DELTA[I + o_delta] * (DPSI + DPHI);`
` 800:                              }`
` 801:                          }`
` 802:                          ETA = B / A;`
` 803:                      }`
` 804:                      else`
` 805:                      {`
` 806:                          if (A <= ZERO)`
` 807:                          {`
` 808:                              ETA = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);`
` 809:                          }`
` 810:                          else`
` 811:                          {`
` 812:                              ETA = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));`
` 813:                          }`
` 814:                      }`
` 815:                  }`
` 816:                  else`
` 817:                  {`
` 818:                      // *`
` 819:                      // *           Interpolation using THREE most relevant poles`
` 820:                      // *`
` 821:                      TEMP = RHOINV + PSI + PHI;`
` 822:                      if (ORGATI)`
` 823:                      {`
` 824:                          TEMP1 = Z[IIM1 + o_z] / DELTA[IIM1 + o_delta];`
` 825:                          TEMP1 = TEMP1 * TEMP1;`
` 826:                          C = TEMP - DELTA[IIP1 + o_delta] * (DPSI + DPHI) - (D[IIM1 + o_d] - D[IIP1 + o_d]) * TEMP1;`
` 827:                          ZZ[1 + o_zz] = Z[IIM1 + o_z] * Z[IIM1 + o_z];`
` 828:                          ZZ[3 + o_zz] = DELTA[IIP1 + o_delta] * DELTA[IIP1 + o_delta] * ((DPSI - TEMP1) + DPHI);`
` 829:                      }`
` 830:                      else`
` 831:                      {`
` 832:                          TEMP1 = Z[IIP1 + o_z] / DELTA[IIP1 + o_delta];`
` 833:                          TEMP1 = TEMP1 * TEMP1;`
` 834:                          C = TEMP - DELTA[IIM1 + o_delta] * (DPSI + DPHI) - (D[IIP1 + o_d] - D[IIM1 + o_d]) * TEMP1;`
` 835:                          ZZ[1 + o_zz] = DELTA[IIM1 + o_delta] * DELTA[IIM1 + o_delta] * (DPSI + (DPHI - TEMP1));`
` 836:                          ZZ[3 + o_zz] = Z[IIP1 + o_z] * Z[IIP1 + o_z];`
` 837:                      }`
` 838:                      ZZ[2 + o_zz] = Z[II + o_z] * Z[II + o_z];`
` 839:                      this._dlaed6.Run(NITER, ORGATI, C, DELTA, IIM1 + o_delta, ZZ, offset_zz, W`
` 840:                                       , ref ETA, ref INFO);`
` 841:                      if (INFO != 0) goto LABEL250;`
` 842:                  }`
` 843:                  // *`
` 844:                  // *        Note, eta should be positive if w is negative, and`
` 845:                  // *        eta should be negative otherwise. However,`
` 846:                  // *        if for some reason caused by roundoff, eta*w > 0,`
` 847:                  // *        we simply use one Newton step instead. This way`
` 848:                  // *        will guarantee eta*w < 0.`
` 849:                  // *`
` 850:                  if (W * ETA >= ZERO) ETA =  - W / DW;`
` 851:                  TEMP = TAU + ETA;`
` 852:                  if (TEMP > DLTUB || TEMP < DLTLB)`
` 853:                  {`
` 854:                      if (W < ZERO)`
` 855:                      {`
` 856:                          ETA = (DLTUB - TAU) / TWO;`
` 857:                      }`
` 858:                      else`
` 859:                      {`
` 860:                          ETA = (DLTLB - TAU) / TWO;`
` 861:                      }`
` 862:                  }`
` 863:                  // *`
` 864:                  PREW = W;`
` 865:                  // *`
` 866:                  for (J = 1; J <= N; J++)`
` 867:                  {`
` 868:                      DELTA[J + o_delta] = DELTA[J + o_delta] - ETA;`
` 869:                  }`
` 870:                  // *`
` 871:                  // *        Evaluate PSI and the derivative DPSI`
` 872:                  // *`
` 873:                  DPSI = ZERO;`
` 874:                  PSI = ZERO;`
` 875:                  ERRETM = ZERO;`
` 876:                  for (J = 1; J <= IIM1; J++)`
` 877:                  {`
` 878:                      TEMP = Z[J + o_z] / DELTA[J + o_delta];`
` 879:                      PSI = PSI + Z[J + o_z] * TEMP;`
` 880:                      DPSI = DPSI + TEMP * TEMP;`
` 881:                      ERRETM = ERRETM + PSI;`
` 882:                  }`
` 883:                  ERRETM = Math.Abs(ERRETM);`
` 884:                  // *`
` 885:                  // *        Evaluate PHI and the derivative DPHI`
` 886:                  // *`
` 887:                  DPHI = ZERO;`
` 888:                  PHI = ZERO;`
` 889:                  for (J = N; J >= IIP1; J +=  - 1)`
` 890:                  {`
` 891:                      TEMP = Z[J + o_z] / DELTA[J + o_delta];`
` 892:                      PHI = PHI + Z[J + o_z] * TEMP;`
` 893:                      DPHI = DPHI + TEMP * TEMP;`
` 894:                      ERRETM = ERRETM + PHI;`
` 895:                  }`
` 896:                  // *`
` 897:                  TEMP = Z[II + o_z] / DELTA[II + o_delta];`
` 898:                  DW = DPSI + DPHI + TEMP * TEMP;`
` 899:                  TEMP = Z[II + o_z] * TEMP;`
` 900:                  W = RHOINV + PHI + PSI + TEMP;`
` 901:                  ERRETM = EIGHT * (PHI - PSI) + ERRETM + TWO * RHOINV + THREE * Math.Abs(TEMP) + Math.Abs(TAU + ETA) * DW;`
` 902:                  // *`
` 903:                  SWTCH = false;`
` 904:                  if (ORGATI)`
` 905:                  {`
` 906:                      if ( - W > Math.Abs(PREW) / TEN) SWTCH = true;`
` 907:                  }`
` 908:                  else`
` 909:                  {`
` 910:                      if (W > Math.Abs(PREW) / TEN) SWTCH = true;`
` 911:                  }`
` 912:                  // *`
` 913:                  TAU = TAU + ETA;`
` 914:                  // *`
` 915:                  // *        Main loop to update the values of the array   DELTA`
` 916:                  // *`
` 917:                  ITER = NITER + 1;`
` 918:                  // *`
` 919:                  for (NITER = ITER; NITER <= MAXIT; NITER++)`
` 920:                  {`
` 921:                      // *`
` 922:                      // *           Test for convergence`
` 923:                      // *`
` 924:                      if (Math.Abs(W) <= EPS * ERRETM)`
` 925:                      {`
` 926:                          if (ORGATI)`
` 927:                          {`
` 928:                              DLAM = D[I + o_d] + TAU;`
` 929:                          }`
` 930:                          else`
` 931:                          {`
` 932:                              DLAM = D[IP1 + o_d] + TAU;`
` 933:                          }`
` 934:                          goto LABEL250;`
` 935:                      }`
` 936:                      // *`
` 937:                      if (W <= ZERO)`
` 938:                      {`
` 939:                          DLTLB = Math.Max(DLTLB, TAU);`
` 940:                      }`
` 941:                      else`
` 942:                      {`
` 943:                          DLTUB = Math.Min(DLTUB, TAU);`
` 944:                      }`
` 945:                      // *`
` 946:                      // *           Calculate the new step`
` 947:                      // *`
` 948:                      if (!SWTCH3)`
` 949:                      {`
` 950:                          if (!SWTCH)`
` 951:                          {`
` 952:                              if (ORGATI)`
` 953:                              {`
` 954:                                  C = W - DELTA[IP1 + o_delta] * DW - (D[I + o_d] - D[IP1 + o_d]) * Math.Pow(Z[I + o_z] / DELTA[I + o_delta],2);`
` 955:                              }`
` 956:                              else`
` 957:                              {`
` 958:                                  C = W - DELTA[I + o_delta] * DW - (D[IP1 + o_d] - D[I + o_d]) * Math.Pow(Z[IP1 + o_z] / DELTA[IP1 + o_delta],2);`
` 959:                              }`
` 960:                          }`
` 961:                          else`
` 962:                          {`
` 963:                              TEMP = Z[II + o_z] / DELTA[II + o_delta];`
` 964:                              if (ORGATI)`
` 965:                              {`
` 966:                                  DPSI = DPSI + TEMP * TEMP;`
` 967:                              }`
` 968:                              else`
` 969:                              {`
` 970:                                  DPHI = DPHI + TEMP * TEMP;`
` 971:                              }`
` 972:                              C = W - DELTA[I + o_delta] * DPSI - DELTA[IP1 + o_delta] * DPHI;`
` 973:                          }`
` 974:                          A = (DELTA[I + o_delta] + DELTA[IP1 + o_delta]) * W - DELTA[I + o_delta] * DELTA[IP1 + o_delta] * DW;`
` 975:                          B = DELTA[I + o_delta] * DELTA[IP1 + o_delta] * W;`
` 976:                          if (C == ZERO)`
` 977:                          {`
` 978:                              if (A == ZERO)`
` 979:                              {`
` 980:                                  if (!SWTCH)`
` 981:                                  {`
` 982:                                      if (ORGATI)`
` 983:                                      {`
` 984:                                          A = Z[I + o_z] * Z[I + o_z] + DELTA[IP1 + o_delta] * DELTA[IP1 + o_delta] * (DPSI + DPHI);`
` 985:                                      }`
` 986:                                      else`
` 987:                                      {`
` 988:                                          A = Z[IP1 + o_z] * Z[IP1 + o_z] + DELTA[I + o_delta] * DELTA[I + o_delta] * (DPSI + DPHI);`
` 989:                                      }`
` 990:                                  }`
` 991:                                  else`
` 992:                                  {`
` 993:                                      A = DELTA[I + o_delta] * DELTA[I + o_delta] * DPSI + DELTA[IP1 + o_delta] * DELTA[IP1 + o_delta] * DPHI;`
` 994:                                  }`
` 995:                              }`
` 996:                              ETA = B / A;`
` 997:                          }`
` 998:                          else`
` 999:                          {`
`1000:                              if (A <= ZERO)`
`1001:                              {`
`1002:                                  ETA = (A - Math.Sqrt(Math.Abs(A * A - FOUR * B * C))) / (TWO * C);`
`1003:                              }`
`1004:                              else`
`1005:                              {`
`1006:                                  ETA = TWO * B / (A + Math.Sqrt(Math.Abs(A * A - FOUR * B * C)));`
`1007:                              }`
`1008:                          }`
`1009:                      }`
`1010:                      else`
`1011:                      {`
`1012:                          // *`
`1013:                          // *              Interpolation using THREE most relevant poles`
`1014:                          // *`
`1015:                          TEMP = RHOINV + PSI + PHI;`
`1016:                          if (SWTCH)`
`1017:                          {`
`1018:                              C = TEMP - DELTA[IIM1 + o_delta] * DPSI - DELTA[IIP1 + o_delta] * DPHI;`
`1019:                              ZZ[1 + o_zz] = DELTA[IIM1 + o_delta] * DELTA[IIM1 + o_delta] * DPSI;`
`1020:                              ZZ[3 + o_zz] = DELTA[IIP1 + o_delta] * DELTA[IIP1 + o_delta] * DPHI;`
`1021:                          }`
`1022:                          else`
`1023:                          {`
`1024:                              if (ORGATI)`
`1025:                              {`
`1026:                                  TEMP1 = Z[IIM1 + o_z] / DELTA[IIM1 + o_delta];`
`1027:                                  TEMP1 = TEMP1 * TEMP1;`
`1028:                                  C = TEMP - DELTA[IIP1 + o_delta] * (DPSI + DPHI) - (D[IIM1 + o_d] - D[IIP1 + o_d]) * TEMP1;`
`1029:                                  ZZ[1 + o_zz] = Z[IIM1 + o_z] * Z[IIM1 + o_z];`
`1030:                                  ZZ[3 + o_zz] = DELTA[IIP1 + o_delta] * DELTA[IIP1 + o_delta] * ((DPSI - TEMP1) + DPHI);`
`1031:                              }`
`1032:                              else`
`1033:                              {`
`1034:                                  TEMP1 = Z[IIP1 + o_z] / DELTA[IIP1 + o_delta];`
`1035:                                  TEMP1 = TEMP1 * TEMP1;`
`1036:                                  C = TEMP - DELTA[IIM1 + o_delta] * (DPSI + DPHI) - (D[IIP1 + o_d] - D[IIM1 + o_d]) * TEMP1;`
`1037:                                  ZZ[1 + o_zz] = DELTA[IIM1 + o_delta] * DELTA[IIM1 + o_delta] * (DPSI + (DPHI - TEMP1));`
`1038:                                  ZZ[3 + o_zz] = Z[IIP1 + o_z] * Z[IIP1 + o_z];`
`1039:                              }`
`1040:                          }`
`1041:                          this._dlaed6.Run(NITER, ORGATI, C, DELTA, IIM1 + o_delta, ZZ, offset_zz, W`
`1042:                                           , ref ETA, ref INFO);`
`1043:                          if (INFO != 0) goto LABEL250;`
`1044:                      }`
`1045:                      // *`
`1046:                      // *           Note, eta should be positive if w is negative, and`
`1047:                      // *           eta should be negative otherwise. However,`
`1048:                      // *           if for some reason caused by roundoff, eta*w > 0,`
`1049:                      // *           we simply use one Newton step instead. This way`
`1050:                      // *           will guarantee eta*w < 0.`
`1051:                      // *`
`1052:                      if (W * ETA >= ZERO) ETA =  - W / DW;`
`1053:                      TEMP = TAU + ETA;`
`1054:                      if (TEMP > DLTUB || TEMP < DLTLB)`
`1055:                      {`
`1056:                          if (W < ZERO)`
`1057:                          {`
`1058:                              ETA = (DLTUB - TAU) / TWO;`
`1059:                          }`
`1060:                          else`
`1061:                          {`
`1062:                              ETA = (DLTLB - TAU) / TWO;`
`1063:                          }`
`1064:                      }`
`1065:                      // *`
`1066:                      for (J = 1; J <= N; J++)`
`1067:                      {`
`1068:                          DELTA[J + o_delta] = DELTA[J + o_delta] - ETA;`
`1069:                      }`
`1070:                      // *`
`1071:                      TAU = TAU + ETA;`
`1072:                      PREW = W;`
`1073:                      // *`
`1074:                      // *           Evaluate PSI and the derivative DPSI`
`1075:                      // *`
`1076:                      DPSI = ZERO;`
`1077:                      PSI = ZERO;`
`1078:                      ERRETM = ZERO;`
`1079:                      for (J = 1; J <= IIM1; J++)`
`1080:                      {`
`1081:                          TEMP = Z[J + o_z] / DELTA[J + o_delta];`
`1082:                          PSI = PSI + Z[J + o_z] * TEMP;`
`1083:                          DPSI = DPSI + TEMP * TEMP;`
`1084:                          ERRETM = ERRETM + PSI;`
`1085:                      }`
`1086:                      ERRETM = Math.Abs(ERRETM);`
`1087:                      // *`
`1088:                      // *           Evaluate PHI and the derivative DPHI`
`1089:                      // *`
`1090:                      DPHI = ZERO;`
`1091:                      PHI = ZERO;`
`1092:                      for (J = N; J >= IIP1; J +=  - 1)`
`1093:                      {`
`1094:                          TEMP = Z[J + o_z] / DELTA[J + o_delta];`
`1095:                          PHI = PHI + Z[J + o_z] * TEMP;`
`1096:                          DPHI = DPHI + TEMP * TEMP;`
`1097:                          ERRETM = ERRETM + PHI;`
`1098:                      }`
`1099:                      // *`
`1100:                      TEMP = Z[II + o_z] / DELTA[II + o_delta];`
`1101:                      DW = DPSI + DPHI + TEMP * TEMP;`
`1102:                      TEMP = Z[II + o_z] * TEMP;`
`1103:                      W = RHOINV + PHI + PSI + TEMP;`
`1104:                      ERRETM = EIGHT * (PHI - PSI) + ERRETM + TWO * RHOINV + THREE * Math.Abs(TEMP) + Math.Abs(TAU) * DW;`
`1105:                      if (W * PREW > ZERO && Math.Abs(W) > Math.Abs(PREW) / TEN) SWTCH = !SWTCH;`
`1106:                      // *`
`1107:                  }`
`1108:                  // *`
`1109:                  // *        Return with INFO = 1, NITER = MAXIT and not converged`
`1110:                  // *`
`1111:                  INFO = 1;`
`1112:                  if (ORGATI)`
`1113:                  {`
`1114:                      DLAM = D[I + o_d] + TAU;`
`1115:                  }`
`1116:                  else`
`1117:                  {`
`1118:                      DLAM = D[IP1 + o_d] + TAU;`
`1119:                  }`
`1120:                  // *`
`1121:              }`
`1122:              // *`
`1123:          LABEL250:;`
`1124:              // *`
`1125:              return;`
`1126:              // *`
`1127:              // *     End of DLAED4`
`1128:              // *`
`1129:   `
`1130:              #endregion`
`1131:   `
`1132:          }`
`1133:      }`
`1134:  }`