`   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 eigenvalue`
`  27:      /// of a positive symmetric rank-one modification of a 2-by-2 diagonal`
`  28:      /// matrix`
`  29:      /// `
`  30:      /// diag( D ) * diag( D ) +  RHO *  Z * transpose(Z) .`
`  31:      /// `
`  32:      /// The diagonal entries in the array D are assumed to satisfy`
`  33:      /// `
`  34:      /// 0 .LE. D(i) .LT. D(j)  for  i .LT. j .`
`  35:      /// `
`  36:      /// We also assume RHO .GT. 0 and that the Euclidean norm of the vector`
`  37:      /// Z is one.`
`  38:      /// `
`  39:      ///</summary>`
`  40:      public class DLASD5`
`  41:      {`
`  42:      `
`  43:   `
`  44:          #region Fields`
`  45:          `
`  46:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; const double TWO = 2.0E+0; const double THREE = 3.0E+0; `
`  47:          const double FOUR = 4.0E+0;double B = 0; double C = 0; double DEL = 0; double DELSQ = 0; double TAU = 0; double W = 0; `
`  48:   `
`  49:          #endregion`
`  50:   `
`  51:          public DLASD5()`
`  52:          {`
`  53:      `
`  54:          }`
`  55:      `
`  56:          /// <summary>`
`  57:          /// Purpose`
`  58:          /// =======`
`  59:          /// `
`  60:          /// This subroutine computes the square root of the I-th eigenvalue`
`  61:          /// of a positive symmetric rank-one modification of a 2-by-2 diagonal`
`  62:          /// matrix`
`  63:          /// `
`  64:          /// diag( D ) * diag( D ) +  RHO *  Z * transpose(Z) .`
`  65:          /// `
`  66:          /// The diagonal entries in the array D are assumed to satisfy`
`  67:          /// `
`  68:          /// 0 .LE. D(i) .LT. D(j)  for  i .LT. j .`
`  69:          /// `
`  70:          /// We also assume RHO .GT. 0 and that the Euclidean norm of the vector`
`  71:          /// Z is one.`
`  72:          /// `
`  73:          ///</summary>`
`  74:          /// <param name="I">`
`  75:          /// (input) INTEGER`
`  76:          /// The index of the eigenvalue to be computed.  I = 1 or I = 2.`
`  77:          ///</param>`
`  78:          /// <param name="D">`
`  79:          /// (input) DOUBLE PRECISION array, dimension ( 2 )`
`  80:          /// The original eigenvalues.  We assume 0 .LE. D(1) .LT. D(2).`
`  81:          ///</param>`
`  82:          /// <param name="Z">`
`  83:          /// (input) DOUBLE PRECISION array, dimension ( 2 )`
`  84:          /// The components of the updating vector.`
`  85:          ///</param>`
`  86:          /// <param name="DELTA">`
`  87:          /// (output) DOUBLE PRECISION array, dimension ( 2 )`
`  88:          /// Contains (D(j) - sigma_I) in its  j-th component.`
`  89:          /// The vector DELTA contains the information necessary`
`  90:          /// to construct the eigenvectors.`
`  91:          ///</param>`
`  92:          /// <param name="RHO">`
`  93:          /// (input) DOUBLE PRECISION`
`  94:          /// The scalar in the symmetric updating formula.`
`  95:          ///</param>`
`  96:          /// <param name="DSIGMA">`
`  97:          /// (output) DOUBLE PRECISION`
`  98:          /// The computed sigma_I, the I-th updated eigenvalue.`
`  99:          ///</param>`
` 100:          /// <param name="WORK">`
` 101:          /// (workspace) DOUBLE PRECISION array, dimension ( 2 )`
` 102:          /// WORK contains (D(j) + sigma_I) in its  j-th component.`
` 103:          ///</param>`
` 104:          public void Run(int I, double[] D, int offset_d, double[] Z, int offset_z, ref double[] DELTA, int offset_delta, double RHO, ref double DSIGMA`
` 105:                           , ref double[] WORK, int offset_work)`
` 106:          {`
` 107:   `
` 108:              #region Array Index Correction`
` 109:              `
` 110:               int o_d = -1 + offset_d;  int o_z = -1 + offset_z;  int o_delta = -1 + offset_delta;  int o_work = -1 + offset_work; `
` 111:   `
` 112:              #endregion`
` 113:   `
` 114:   `
` 115:              #region Prolog`
` 116:              `
` 117:              // *`
` 118:              // *  -- LAPACK auxiliary routine (version 3.1) --`
` 119:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 120:              // *     November 2006`
` 121:              // *`
` 122:              // *     .. Scalar Arguments ..`
` 123:              // *     ..`
` 124:              // *     .. Array Arguments ..`
` 125:              // *     ..`
` 126:              // *`
` 127:              // *  Purpose`
` 128:              // *  =======`
` 129:              // *`
` 130:              // *  This subroutine computes the square root of the I-th eigenvalue`
` 131:              // *  of a positive symmetric rank-one modification of a 2-by-2 diagonal`
` 132:              // *  matrix`
` 133:              // *`
` 134:              // *             diag( D ) * diag( D ) +  RHO *  Z * transpose(Z) .`
` 135:              // *`
` 136:              // *  The diagonal entries in the array D are assumed to satisfy`
` 137:              // *`
` 138:              // *             0 <= D(i) < D(j)  for  i < j .`
` 139:              // *`
` 140:              // *  We also assume RHO > 0 and that the Euclidean norm of the vector`
` 141:              // *  Z is one.`
` 142:              // *`
` 143:              // *  Arguments`
` 144:              // *  =========`
` 145:              // *`
` 146:              // *  I      (input) INTEGER`
` 147:              // *         The index of the eigenvalue to be computed.  I = 1 or I = 2.`
` 148:              // *`
` 149:              // *  D      (input) DOUBLE PRECISION array, dimension ( 2 )`
` 150:              // *         The original eigenvalues.  We assume 0 <= D(1) < D(2).`
` 151:              // *`
` 152:              // *  Z      (input) DOUBLE PRECISION array, dimension ( 2 )`
` 153:              // *         The components of the updating vector.`
` 154:              // *`
` 155:              // *  DELTA  (output) DOUBLE PRECISION array, dimension ( 2 )`
` 156:              // *         Contains (D(j) - sigma_I) in its  j-th component.`
` 157:              // *         The vector DELTA contains the information necessary`
` 158:              // *         to construct the eigenvectors.`
` 159:              // *`
` 160:              // *  RHO    (input) DOUBLE PRECISION`
` 161:              // *         The scalar in the symmetric updating formula.`
` 162:              // *`
` 163:              // *  DSIGMA (output) DOUBLE PRECISION`
` 164:              // *         The computed sigma_I, the I-th updated eigenvalue.`
` 165:              // *`
` 166:              // *  WORK   (workspace) DOUBLE PRECISION array, dimension ( 2 )`
` 167:              // *         WORK contains (D(j) + sigma_I) in its  j-th component.`
` 168:              // *`
` 169:              // *  Further Details`
` 170:              // *  ===============`
` 171:              // *`
` 172:              // *  Based on contributions by`
` 173:              // *     Ren-Cang Li, Computer Science Division, University of California`
` 174:              // *     at Berkeley, USA`
` 175:              // *`
` 176:              // *  =====================================================================`
` 177:              // *`
` 178:              // *     .. Parameters ..`
` 179:              // *     ..`
` 180:              // *     .. Local Scalars ..`
` 181:              // *     ..`
` 182:              // *     .. Intrinsic Functions ..`
` 183:              //      INTRINSIC          ABS, SQRT;`
` 184:              // *     ..`
` 185:              // *     .. Executable Statements ..`
` 186:              // *`
` 187:   `
` 188:              #endregion`
` 189:   `
` 190:   `
` 191:              #region Body`
` 192:              `
` 193:              DEL = D[2 + o_d] - D[1 + o_d];`
` 194:              DELSQ = DEL * (D[2 + o_d] + D[1 + o_d]);`
` 195:              if (I == 1)`
` 196:              {`
` 197:                  W = ONE + FOUR * RHO * (Z[2 + o_z] * Z[2 + o_z] / (D[1 + o_d] + THREE * D[2 + o_d]) - Z[1 + o_z] * Z[1 + o_z] / (THREE * D[1 + o_d] + D[2 + o_d])) / DEL;`
` 198:                  if (W > ZERO)`
` 199:                  {`
` 200:                      B = DELSQ + RHO * (Z[1 + o_z] * Z[1 + o_z] + Z[2 + o_z] * Z[2 + o_z]);`
` 201:                      C = RHO * Z[1 + o_z] * Z[1 + o_z] * DELSQ;`
` 202:                      // *`
` 203:                      // *           B > ZERO, always`
` 204:                      // *`
` 205:                      // *           The following TAU is DSIGMA * DSIGMA - D( 1 ) * D( 1 )`
` 206:                      // *`
` 207:                      TAU = TWO * C / (B + Math.Sqrt(Math.Abs(B * B - FOUR * C)));`
` 208:                      // *`
` 209:                      // *           The following TAU is DSIGMA - D( 1 )`
` 210:                      // *`
` 211:                      TAU = TAU / (D[1 + o_d] + Math.Sqrt(D[1 + o_d] * D[1 + o_d] + TAU));`
` 212:                      DSIGMA = D[1 + o_d] + TAU;`
` 213:                      DELTA[1 + o_delta] =  - TAU;`
` 214:                      DELTA[2 + o_delta] = DEL - TAU;`
` 215:                      WORK[1 + o_work] = TWO * D[1 + o_d] + TAU;`
` 216:                      WORK[2 + o_work] = (D[1 + o_d] + TAU) + D[2 + o_d];`
` 217:                      // *           DELTA( 1 ) = -Z( 1 ) / TAU`
` 218:                      // *           DELTA( 2 ) = Z( 2 ) / ( DEL-TAU )`
` 219:                  }`
` 220:                  else`
` 221:                  {`
` 222:                      B =  - DELSQ + RHO * (Z[1 + o_z] * Z[1 + o_z] + Z[2 + o_z] * Z[2 + o_z]);`
` 223:                      C = RHO * Z[2 + o_z] * Z[2 + o_z] * DELSQ;`
` 224:                      // *`
` 225:                      // *           The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )`
` 226:                      // *`
` 227:                      if (B > ZERO)`
` 228:                      {`
` 229:                          TAU =  - TWO * C / (B + Math.Sqrt(B * B + FOUR * C));`
` 230:                      }`
` 231:                      else`
` 232:                      {`
` 233:                          TAU = (B - Math.Sqrt(B * B + FOUR * C)) / TWO;`
` 234:                      }`
` 235:                      // *`
` 236:                      // *           The following TAU is DSIGMA - D( 2 )`
` 237:                      // *`
` 238:                      TAU = TAU / (D[2 + o_d] + Math.Sqrt(Math.Abs(D[2 + o_d] * D[2 + o_d] + TAU)));`
` 239:                      DSIGMA = D[2 + o_d] + TAU;`
` 240:                      DELTA[1 + o_delta] =  - (DEL + TAU);`
` 241:                      DELTA[2 + o_delta] =  - TAU;`
` 242:                      WORK[1 + o_work] = D[1 + o_d] + TAU + D[2 + o_d];`
` 243:                      WORK[2 + o_work] = TWO * D[2 + o_d] + TAU;`
` 244:                      // *           DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )`
` 245:                      // *           DELTA( 2 ) = -Z( 2 ) / TAU`
` 246:                  }`
` 247:                  // *        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )`
` 248:                  // *        DELTA( 1 ) = DELTA( 1 ) / TEMP`
` 249:                  // *        DELTA( 2 ) = DELTA( 2 ) / TEMP`
` 250:              }`
` 251:              else`
` 252:              {`
` 253:                  // *`
` 254:                  // *        Now I=2`
` 255:                  // *`
` 256:                  B =  - DELSQ + RHO * (Z[1 + o_z] * Z[1 + o_z] + Z[2 + o_z] * Z[2 + o_z]);`
` 257:                  C = RHO * Z[2 + o_z] * Z[2 + o_z] * DELSQ;`
` 258:                  // *`
` 259:                  // *        The following TAU is DSIGMA * DSIGMA - D( 2 ) * D( 2 )`
` 260:                  // *`
` 261:                  if (B > ZERO)`
` 262:                  {`
` 263:                      TAU = (B + Math.Sqrt(B * B + FOUR * C)) / TWO;`
` 264:                  }`
` 265:                  else`
` 266:                  {`
` 267:                      TAU = TWO * C / ( - B + Math.Sqrt(B * B + FOUR * C));`
` 268:                  }`
` 269:                  // *`
` 270:                  // *        The following TAU is DSIGMA - D( 2 )`
` 271:                  // *`
` 272:                  TAU = TAU / (D[2 + o_d] + Math.Sqrt(D[2 + o_d] * D[2 + o_d] + TAU));`
` 273:                  DSIGMA = D[2 + o_d] + TAU;`
` 274:                  DELTA[1 + o_delta] =  - (DEL + TAU);`
` 275:                  DELTA[2 + o_delta] =  - TAU;`
` 276:                  WORK[1 + o_work] = D[1 + o_d] + TAU + D[2 + o_d];`
` 277:                  WORK[2 + o_work] = TWO * D[2 + o_d] + TAU;`
` 278:                  // *        DELTA( 1 ) = -Z( 1 ) / ( DEL+TAU )`
` 279:                  // *        DELTA( 2 ) = -Z( 2 ) / TAU`
` 280:                  // *        TEMP = SQRT( DELTA( 1 )*DELTA( 1 )+DELTA( 2 )*DELTA( 2 ) )`
` 281:                  // *        DELTA( 1 ) = DELTA( 1 ) / TEMP`
` 282:                  // *        DELTA( 2 ) = DELTA( 2 ) / TEMP`
` 283:              }`
` 284:              return;`
` 285:              // *`
` 286:              // *     End of DLASD5`
` 287:              // *`
` 288:   `
` 289:              #endregion`
` 290:   `
` 291:          }`
` 292:      }`
` 293:  }`