  `   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:      /// DLAED1 computes the updated eigensystem of a diagonal`
`  27:      /// matrix after modification by a rank-one symmetric matrix.  This`
`  28:      /// routine is used only for the eigenproblem which requires all`
`  29:      /// eigenvalues and eigenvectors of a tridiagonal matrix.  DLAED7 handles`
`  30:      /// the case in which eigenvalues only or eigenvalues and eigenvectors`
`  31:      /// of a full symmetric matrix (which was reduced to tridiagonal form)`
`  32:      /// are desired.`
`  33:      /// `
`  34:      /// T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)`
`  35:      /// `
`  36:      /// where Z = Q'u, u is a vector of length N with ones in the`
`  37:      /// CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.`
`  38:      /// `
`  39:      /// The eigenvectors of the original matrix are stored in Q, and the`
`  40:      /// eigenvalues are in D.  The algorithm consists of three stages:`
`  41:      /// `
`  42:      /// The first stage consists of deflating the size of the problem`
`  43:      /// when there are multiple eigenvalues or if there is a zero in`
`  44:      /// the Z vector.  For each such occurence the dimension of the`
`  45:      /// secular equation problem is reduced by one.  This stage is`
`  46:      /// performed by the routine DLAED2.`
`  47:      /// `
`  48:      /// The second stage consists of calculating the updated`
`  49:      /// eigenvalues. This is done by finding the roots of the secular`
`  50:      /// equation via the routine DLAED4 (as called by DLAED3).`
`  51:      /// This routine also calculates the eigenvectors of the current`
`  52:      /// problem.`
`  53:      /// `
`  54:      /// The final stage consists of computing the updated eigenvectors`
`  55:      /// directly using the updated eigenvalues.  The eigenvectors for`
`  56:      /// the current problem are multiplied with the eigenvectors from`
`  57:      /// the overall problem.`
`  58:      /// `
`  59:      ///</summary>`
`  60:      public class DLAED1`
`  61:      {`
`  62:      `
`  63:   `
`  64:          #region Dependencies`
`  65:          `
`  66:          DCOPY _dcopy; DLAED2 _dlaed2; DLAED3 _dlaed3; DLAMRG _dlamrg; XERBLA _xerbla; `
`  67:   `
`  68:          #endregion`
`  69:   `
`  70:   `
`  71:          #region Fields`
`  72:          `
`  73:          int COLTYP = 0; int I = 0; int IDLMDA = 0; int INDX = 0; int INDXC = 0; int INDXP = 0; int IQ2 = 0; int IS = 0; `
`  74:          int IW = 0;int IZ = 0; int K = 0; int N1 = 0; int N2 = 0; int ZPP1 = 0; `
`  75:   `
`  76:          #endregion`
`  77:   `
`  78:          public DLAED1(DCOPY dcopy, DLAED2 dlaed2, DLAED3 dlaed3, DLAMRG dlamrg, XERBLA xerbla)`
`  79:          {`
`  80:      `
`  81:   `
`  82:              #region Set Dependencies`
`  83:              `
`  84:              this._dcopy = dcopy; this._dlaed2 = dlaed2; this._dlaed3 = dlaed3; this._dlamrg = dlamrg; this._xerbla = xerbla; `
`  85:   `
`  86:              #endregion`
`  87:   `
`  88:          }`
`  89:      `
`  90:          public DLAED1()`
`  91:          {`
`  92:      `
`  93:   `
`  94:              #region Dependencies (Initialization)`
`  95:              `
`  96:              DCOPY dcopy = new DCOPY();`
`  97:              IDAMAX idamax = new IDAMAX();`
`  98:              LSAME lsame = new LSAME();`
`  99:              DLAMC3 dlamc3 = new DLAMC3();`
` 100:              DLAPY2 dlapy2 = new DLAPY2();`
` 101:              DLAMRG dlamrg = new DLAMRG();`
` 102:              DROT drot = new DROT();`
` 103:              DSCAL dscal = new DSCAL();`
` 104:              XERBLA xerbla = new XERBLA();`
` 105:              DNRM2 dnrm2 = new DNRM2();`
` 106:              DLAED5 dlaed5 = new DLAED5();`
` 107:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 108:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 109:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 110:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 111:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 112:              DLACPY dlacpy = new DLACPY(lsame);`
` 113:              DLAED2 dlaed2 = new DLAED2(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);`
` 114:              DGEMM dgemm = new DGEMM(lsame, xerbla);`
` 115:              DLAED6 dlaed6 = new DLAED6(dlamch);`
` 116:              DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);`
` 117:              DLASET dlaset = new DLASET(lsame);`
` 118:              DLAED3 dlaed3 = new DLAED3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlaed4, dlaset, xerbla);`
` 119:   `
` 120:              #endregion`
` 121:   `
` 122:   `
` 123:              #region Set Dependencies`
` 124:              `
` 125:              this._dcopy = dcopy; this._dlaed2 = dlaed2; this._dlaed3 = dlaed3; this._dlamrg = dlamrg; this._xerbla = xerbla; `
` 126:   `
` 127:              #endregion`
` 128:   `
` 129:          }`
` 130:          /// <summary>`
` 131:          /// Purpose`
` 132:          /// =======`
` 133:          /// `
` 134:          /// DLAED1 computes the updated eigensystem of a diagonal`
` 135:          /// matrix after modification by a rank-one symmetric matrix.  This`
` 136:          /// routine is used only for the eigenproblem which requires all`
` 137:          /// eigenvalues and eigenvectors of a tridiagonal matrix.  DLAED7 handles`
` 138:          /// the case in which eigenvalues only or eigenvalues and eigenvectors`
` 139:          /// of a full symmetric matrix (which was reduced to tridiagonal form)`
` 140:          /// are desired.`
` 141:          /// `
` 142:          /// T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)`
` 143:          /// `
` 144:          /// where Z = Q'u, u is a vector of length N with ones in the`
` 145:          /// CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.`
` 146:          /// `
` 147:          /// The eigenvectors of the original matrix are stored in Q, and the`
` 148:          /// eigenvalues are in D.  The algorithm consists of three stages:`
` 149:          /// `
` 150:          /// The first stage consists of deflating the size of the problem`
` 151:          /// when there are multiple eigenvalues or if there is a zero in`
` 152:          /// the Z vector.  For each such occurence the dimension of the`
` 153:          /// secular equation problem is reduced by one.  This stage is`
` 154:          /// performed by the routine DLAED2.`
` 155:          /// `
` 156:          /// The second stage consists of calculating the updated`
` 157:          /// eigenvalues. This is done by finding the roots of the secular`
` 158:          /// equation via the routine DLAED4 (as called by DLAED3).`
` 159:          /// This routine also calculates the eigenvectors of the current`
` 160:          /// problem.`
` 161:          /// `
` 162:          /// The final stage consists of computing the updated eigenvectors`
` 163:          /// directly using the updated eigenvalues.  The eigenvectors for`
` 164:          /// the current problem are multiplied with the eigenvectors from`
` 165:          /// the overall problem.`
` 166:          /// `
` 167:          ///</summary>`
` 168:          /// <param name="N">`
` 169:          /// (input) INTEGER`
` 170:          /// The dimension of the symmetric tridiagonal matrix.  N .GE. 0.`
` 171:          ///</param>`
` 172:          /// <param name="D">`
` 173:          /// (input/output) DOUBLE PRECISION array, dimension (N)`
` 174:          /// On entry, the eigenvalues of the rank-1-perturbed matrix.`
` 175:          /// On exit, the eigenvalues of the repaired matrix.`
` 176:          ///</param>`
` 177:          /// <param name="Q">`
` 178:          /// (input/output) DOUBLE PRECISION array, dimension (LDQ,N)`
` 179:          /// On entry, the eigenvectors of the rank-1-perturbed matrix.`
` 180:          /// On exit, the eigenvectors of the repaired tridiagonal matrix.`
` 181:          ///</param>`
` 182:          /// <param name="LDQ">`
` 183:          /// (input) INTEGER`
` 184:          /// The leading dimension of the array Q.  LDQ .GE. max(1,N).`
` 185:          ///</param>`
` 186:          /// <param name="INDXQ">`
` 187:          /// (input/output) INTEGER array, dimension (N)`
` 188:          /// On entry, the permutation which separately sorts the two`
` 189:          /// subproblems in D into ascending order.`
` 190:          /// On exit, the permutation which will reintegrate the`
` 191:          /// subproblems back into sorted order,`
` 192:          /// i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.`
` 193:          ///</param>`
` 194:          /// <param name="RHO">`
` 195:          /// (input) DOUBLE PRECISION`
` 196:          /// The subdiagonal entry used to create the rank-1 modification.`
` 197:          ///</param>`
` 198:          /// <param name="CUTPNT">`
` 199:          /// (input) INTEGER`
` 200:          /// The location of the last eigenvalue in the leading sub-matrix.`
` 201:          /// min(1,N) .LE. CUTPNT .LE. N/2.`
` 202:          ///</param>`
` 203:          /// <param name="WORK">`
` 204:          /// (workspace) DOUBLE PRECISION array, dimension (4*N + N**2)`
` 205:          ///</param>`
` 206:          /// <param name="IWORK">`
` 207:          /// (workspace) INTEGER array, dimension (4*N)`
` 208:          ///</param>`
` 209:          /// <param name="INFO">`
` 210:          /// (output) INTEGER`
` 211:          /// = 0:  successful exit.`
` 212:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.`
` 213:          /// .GT. 0:  if INFO = 1, an eigenvalue did not converge`
` 214:          ///</param>`
` 215:          public void Run(int N, ref double[] D, int offset_d, ref double[] Q, int offset_q, int LDQ, ref int[] INDXQ, int offset_indxq, ref double RHO`
` 216:                           , int CUTPNT, ref double[] WORK, int offset_work, ref int[] IWORK, int offset_iwork, ref int INFO)`
` 217:          {`
` 218:   `
` 219:              #region Array Index Correction`
` 220:              `
` 221:               int o_d = -1 + offset_d;  int o_q = -1 - LDQ + offset_q;  int o_indxq = -1 + offset_indxq; `
` 222:               int o_work = -1 + offset_work; int o_iwork = -1 + offset_iwork; `
` 223:   `
` 224:              #endregion`
` 225:   `
` 226:   `
` 227:              #region Prolog`
` 228:              `
` 229:              // *`
` 230:              // *  -- LAPACK routine (version 3.1) --`
` 231:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 232:              // *     November 2006`
` 233:              // *`
` 234:              // *     .. Scalar Arguments ..`
` 235:              // *     ..`
` 236:              // *     .. Array Arguments ..`
` 237:              // *     ..`
` 238:              // *`
` 239:              // *  Purpose`
` 240:              // *  =======`
` 241:              // *`
` 242:              // *  DLAED1 computes the updated eigensystem of a diagonal`
` 243:              // *  matrix after modification by a rank-one symmetric matrix.  This`
` 244:              // *  routine is used only for the eigenproblem which requires all`
` 245:              // *  eigenvalues and eigenvectors of a tridiagonal matrix.  DLAED7 handles`
` 246:              // *  the case in which eigenvalues only or eigenvalues and eigenvectors`
` 247:              // *  of a full symmetric matrix (which was reduced to tridiagonal form)`
` 248:              // *  are desired.`
` 249:              // *`
` 250:              // *    T = Q(in) ( D(in) + RHO * Z*Z' ) Q'(in) = Q(out) * D(out) * Q'(out)`
` 251:              // *`
` 252:              // *     where Z = Q'u, u is a vector of length N with ones in the`
` 253:              // *     CUTPNT and CUTPNT + 1 th elements and zeros elsewhere.`
` 254:              // *`
` 255:              // *     The eigenvectors of the original matrix are stored in Q, and the`
` 256:              // *     eigenvalues are in D.  The algorithm consists of three stages:`
` 257:              // *`
` 258:              // *        The first stage consists of deflating the size of the problem`
` 259:              // *        when there are multiple eigenvalues or if there is a zero in`
` 260:              // *        the Z vector.  For each such occurence the dimension of the`
` 261:              // *        secular equation problem is reduced by one.  This stage is`
` 262:              // *        performed by the routine DLAED2.`
` 263:              // *`
` 264:              // *        The second stage consists of calculating the updated`
` 265:              // *        eigenvalues. This is done by finding the roots of the secular`
` 266:              // *        equation via the routine DLAED4 (as called by DLAED3).`
` 267:              // *        This routine also calculates the eigenvectors of the current`
` 268:              // *        problem.`
` 269:              // *`
` 270:              // *        The final stage consists of computing the updated eigenvectors`
` 271:              // *        directly using the updated eigenvalues.  The eigenvectors for`
` 272:              // *        the current problem are multiplied with the eigenvectors from`
` 273:              // *        the overall problem.`
` 274:              // *`
` 275:              // *  Arguments`
` 276:              // *  =========`
` 277:              // *`
` 278:              // *  N      (input) INTEGER`
` 279:              // *         The dimension of the symmetric tridiagonal matrix.  N >= 0.`
` 280:              // *`
` 281:              // *  D      (input/output) DOUBLE PRECISION array, dimension (N)`
` 282:              // *         On entry, the eigenvalues of the rank-1-perturbed matrix.`
` 283:              // *         On exit, the eigenvalues of the repaired matrix.`
` 284:              // *`
` 285:              // *  Q      (input/output) DOUBLE PRECISION array, dimension (LDQ,N)`
` 286:              // *         On entry, the eigenvectors of the rank-1-perturbed matrix.`
` 287:              // *         On exit, the eigenvectors of the repaired tridiagonal matrix.`
` 288:              // *`
` 289:              // *  LDQ    (input) INTEGER`
` 290:              // *         The leading dimension of the array Q.  LDQ >= max(1,N).`
` 291:              // *`
` 292:              // *  INDXQ  (input/output) INTEGER array, dimension (N)`
` 293:              // *         On entry, the permutation which separately sorts the two`
` 294:              // *         subproblems in D into ascending order.`
` 295:              // *         On exit, the permutation which will reintegrate the`
` 296:              // *         subproblems back into sorted order,`
` 297:              // *         i.e. D( INDXQ( I = 1, N ) ) will be in ascending order.`
` 298:              // *`
` 299:              // *  RHO    (input) DOUBLE PRECISION`
` 300:              // *         The subdiagonal entry used to create the rank-1 modification.`
` 301:              // *`
` 302:              // *  CUTPNT (input) INTEGER`
` 303:              // *         The location of the last eigenvalue in the leading sub-matrix.`
` 304:              // *         min(1,N) <= CUTPNT <= N/2.`
` 305:              // *`
` 306:              // *  WORK   (workspace) DOUBLE PRECISION array, dimension (4*N + N**2)`
` 307:              // *`
` 308:              // *  IWORK  (workspace) INTEGER array, dimension (4*N)`
` 309:              // *`
` 310:              // *  INFO   (output) INTEGER`
` 311:              // *          = 0:  successful exit.`
` 312:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.`
` 313:              // *          > 0:  if INFO = 1, an eigenvalue did not converge`
` 314:              // *`
` 315:              // *  Further Details`
` 316:              // *  ===============`
` 317:              // *`
` 318:              // *  Based on contributions by`
` 319:              // *     Jeff Rutter, Computer Science Division, University of California`
` 320:              // *     at Berkeley, USA`
` 321:              // *  Modified by Francoise Tisseur, University of Tennessee.`
` 322:              // *`
` 323:              // *  =====================================================================`
` 324:              // *`
` 325:              // *     .. Local Scalars ..`
` 326:              // *     ..`
` 327:              // *     .. External Subroutines ..`
` 328:              // *     ..`
` 329:              // *     .. Intrinsic Functions ..`
` 330:              //      INTRINSIC          MAX, MIN;`
` 331:              // *     ..`
` 332:              // *     .. Executable Statements ..`
` 333:              // *`
` 334:              // *     Test the input parameters.`
` 335:              // *`
` 336:   `
` 337:              #endregion`
` 338:   `
` 339:   `
` 340:              #region Body`
` 341:              `
` 342:              INFO = 0;`
` 343:              // *`
` 344:              if (N < 0)`
` 345:              {`
` 346:                  INFO =  - 1;`
` 347:              }`
` 348:              else`
` 349:              {`
` 350:                  if (LDQ < Math.Max(1, N))`
` 351:                  {`
` 352:                      INFO =  - 4;`
` 353:                  }`
` 354:                  else`
` 355:                  {`
` 356:                      if (Math.Min(1, N / 2) > CUTPNT || (N / 2) < CUTPNT)`
` 357:                      {`
` 358:                          INFO =  - 7;`
` 359:                      }`
` 360:                  }`
` 361:              }`
` 362:              if (INFO != 0)`
` 363:              {`
` 364:                  this._xerbla.Run("DLAED1",  - INFO);`
` 365:                  return;`
` 366:              }`
` 367:              // *`
` 368:              // *     Quick return if possible`
` 369:              // *`
` 370:              if (N == 0) return;`
` 371:              // *`
` 372:              // *     The following values are integer pointers which indicate`
` 373:              // *     the portion of the workspace`
` 374:              // *     used by a particular array in DLAED2 and DLAED3.`
` 375:              // *`
` 376:              IZ = 1;`
` 377:              IDLMDA = IZ + N;`
` 378:              IW = IDLMDA + N;`
` 379:              IQ2 = IW + N;`
` 380:              // *`
` 381:              INDX = 1;`
` 382:              INDXC = INDX + N;`
` 383:              COLTYP = INDXC + N;`
` 384:              INDXP = COLTYP + N;`
` 385:              // *`
` 386:              // *`
` 387:              // *     Form the z-vector which consists of the last row of Q_1 and the`
` 388:              // *     first row of Q_2.`
` 389:              // *`
` 390:              this._dcopy.Run(CUTPNT, Q, CUTPNT+1 * LDQ + o_q, LDQ, ref WORK, IZ + o_work, 1);`
` 391:              ZPP1 = CUTPNT + 1;`
` 392:              this._dcopy.Run(N - CUTPNT, Q, ZPP1+ZPP1 * LDQ + o_q, LDQ, ref WORK, IZ + CUTPNT + o_work, 1);`
` 393:              // *`
` 394:              // *     Deflate eigenvalues.`
` 395:              // *`
` 396:              this._dlaed2.Run(ref K, N, CUTPNT, ref D, offset_d, ref Q, offset_q, LDQ`
` 397:                               , ref INDXQ, offset_indxq, ref RHO, ref WORK, IZ + o_work, ref WORK, IDLMDA + o_work, ref WORK, IW + o_work, ref WORK, IQ2 + o_work`
` 398:                               , ref IWORK, INDX + o_iwork, ref IWORK, INDXC + o_iwork, ref IWORK, INDXP + o_iwork, ref IWORK, COLTYP + o_iwork, ref INFO);`
` 399:              // *`
` 400:              if (INFO != 0) goto LABEL20;`
` 401:              // *`
` 402:              // *     Solve Secular Equation.`
` 403:              // *`
` 404:              if (K != 0)`
` 405:              {`
` 406:                  IS = (IWORK[COLTYP + o_iwork] + IWORK[COLTYP + 1 + o_iwork]) * CUTPNT + (IWORK[COLTYP + 1 + o_iwork] + IWORK[COLTYP + 2 + o_iwork]) * (N - CUTPNT) + IQ2;`
` 407:                  this._dlaed3.Run(K, N, CUTPNT, ref D, offset_d, ref Q, offset_q, LDQ`
` 408:                                   , RHO, ref WORK, IDLMDA + o_work, WORK, IQ2 + o_work, IWORK, INDXC + o_iwork, IWORK, COLTYP + o_iwork, ref WORK, IW + o_work`
` 409:                                   , ref WORK, IS + o_work, ref INFO);`
` 410:                  if (INFO != 0) goto LABEL20;`
` 411:                  // *`
` 412:                  // *     Prepare the INDXQ sorting permutation.`
` 413:                  // *`
` 414:                  N1 = K;`
` 415:                  N2 = N - K;`
` 416:                  this._dlamrg.Run(N1, N2, D, offset_d, 1,  - 1, ref INDXQ, offset_indxq);`
` 417:              }`
` 418:              else`
` 419:              {`
` 420:                  for (I = 1; I <= N; I++)`
` 421:                  {`
` 422:                      INDXQ[I + o_indxq] = I;`
` 423:                  }`
` 424:              }`
` 425:              // *`
` 426:          LABEL20:;`
` 427:              return;`
` 428:              // *`
` 429:              // *     End of DLAED1`
` 430:              // *`
` 431:   `
` 432:              #endregion`
` 433:   `
` 434:          }`
` 435:      }`
` 436:  }`