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