Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   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:      /// DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
  27:      /// symmetric tridiagonal matrix using the implicit QL or QR method.
  28:      /// The eigenvectors of a full or band symmetric matrix can also be found
  29:      /// if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
  30:      /// tridiagonal form.
  31:      /// 
  32:      ///</summary>
  33:      public class DSTEQR
  34:      {
  35:      
  36:   
  37:          #region Dependencies
  38:          
  39:          LSAME _lsame; DLAMCH _dlamch; DLANST _dlanst; DLAPY2 _dlapy2; DLAE2 _dlae2; DLAEV2 _dlaev2; DLARTG _dlartg; 
  40:          DLASCL _dlascl;DLASET _dlaset; DLASR _dlasr; DLASRT _dlasrt; DSWAP _dswap; XERBLA _xerbla; 
  41:   
  42:          #endregion
  43:   
  44:   
  45:          #region Fields
  46:          
  47:          const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; const double THREE = 3.0E0; 
  48:          const int MAXIT = 30;int I = 0; int ICOMPZ = 0; int II = 0; int ISCALE = 0; int J = 0; int JTOT = 0; int K = 0; int L = 0; 
  49:          int L1 = 0;int LEND = 0; int LENDM1 = 0; int LENDP1 = 0; int LENDSV = 0; int LM1 = 0; int LSV = 0; int M = 0; int MM = 0; 
  50:          int MM1 = 0;int NM1 = 0; int NMAXIT = 0; double ANORM = 0; double B = 0; double C = 0; double EPS = 0; double EPS2 = 0; 
  51:          double F = 0;double G = 0; double P = 0; double R = 0; double RT1 = 0; double RT2 = 0; double S = 0; double SAFMAX = 0; 
  52:          double SAFMIN = 0;double SSFMAX = 0; double SSFMIN = 0; double TST = 0; 
  53:   
  54:          #endregion
  55:   
  56:          public DSTEQR(LSAME lsame, DLAMCH dlamch, DLANST dlanst, DLAPY2 dlapy2, DLAE2 dlae2, DLAEV2 dlaev2, DLARTG dlartg, DLASCL dlascl, DLASET dlaset, DLASR dlasr
  57:                        , DLASRT dlasrt, DSWAP dswap, XERBLA xerbla)
  58:          {
  59:      
  60:   
  61:              #region Set Dependencies
  62:              
  63:              this._lsame = lsame; this._dlamch = dlamch; this._dlanst = dlanst; this._dlapy2 = dlapy2; this._dlae2 = dlae2; 
  64:              this._dlaev2 = dlaev2;this._dlartg = dlartg; this._dlascl = dlascl; this._dlaset = dlaset; this._dlasr = dlasr; 
  65:              this._dlasrt = dlasrt;this._dswap = dswap; this._xerbla = xerbla; 
  66:   
  67:              #endregion
  68:   
  69:          }
  70:      
  71:          public DSTEQR()
  72:          {
  73:      
  74:   
  75:              #region Dependencies (Initialization)
  76:              
  77:              LSAME lsame = new LSAME();
  78:              DLAMC3 dlamc3 = new DLAMC3();
  79:              DLASSQ dlassq = new DLASSQ();
  80:              DLAPY2 dlapy2 = new DLAPY2();
  81:              DLAE2 dlae2 = new DLAE2();
  82:              DLAEV2 dlaev2 = new DLAEV2();
  83:              XERBLA xerbla = new XERBLA();
  84:              DSWAP dswap = new DSWAP();
  85:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  86:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  87:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  88:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  89:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  90:              DLANST dlanst = new DLANST(lsame, dlassq);
  91:              DLARTG dlartg = new DLARTG(dlamch);
  92:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
  93:              DLASET dlaset = new DLASET(lsame);
  94:              DLASR dlasr = new DLASR(lsame, xerbla);
  95:              DLASRT dlasrt = new DLASRT(lsame, xerbla);
  96:   
  97:              #endregion
  98:   
  99:   
 100:              #region Set Dependencies
 101:              
 102:              this._lsame = lsame; this._dlamch = dlamch; this._dlanst = dlanst; this._dlapy2 = dlapy2; this._dlae2 = dlae2; 
 103:              this._dlaev2 = dlaev2;this._dlartg = dlartg; this._dlascl = dlascl; this._dlaset = dlaset; this._dlasr = dlasr; 
 104:              this._dlasrt = dlasrt;this._dswap = dswap; this._xerbla = xerbla; 
 105:   
 106:              #endregion
 107:   
 108:          }
 109:          /// <summary>
 110:          /// Purpose
 111:          /// =======
 112:          /// 
 113:          /// DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
 114:          /// symmetric tridiagonal matrix using the implicit QL or QR method.
 115:          /// The eigenvectors of a full or band symmetric matrix can also be found
 116:          /// if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
 117:          /// tridiagonal form.
 118:          /// 
 119:          ///</summary>
 120:          /// <param name="COMPZ">
 121:          /// (input) CHARACTER*1
 122:          /// = 'N':  Compute eigenvalues only.
 123:          /// = 'V':  Compute eigenvalues and eigenvectors of the original
 124:          /// symmetric matrix.  On entry, Z must contain the
 125:          /// orthogonal matrix used to reduce the original matrix
 126:          /// to tridiagonal form.
 127:          /// = 'I':  Compute eigenvalues and eigenvectors of the
 128:          /// tridiagonal matrix.  Z is initialized to the identity
 129:          /// matrix.
 130:          ///</param>
 131:          /// <param name="N">
 132:          /// (input) INTEGER
 133:          /// The order of the matrix.  N .GE. 0.
 134:          ///</param>
 135:          /// <param name="D">
 136:          /// (input/output) DOUBLE PRECISION array, dimension (N)
 137:          /// On entry, the diagonal elements of the tridiagonal matrix.
 138:          /// On exit, if INFO = 0, the eigenvalues in ascending order.
 139:          ///</param>
 140:          /// <param name="E">
 141:          /// (input/output) DOUBLE PRECISION array, dimension (N-1)
 142:          /// On entry, the (n-1) subdiagonal elements of the tridiagonal
 143:          /// matrix.
 144:          /// On exit, E has been destroyed.
 145:          ///</param>
 146:          /// <param name="Z">
 147:          /// (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
 148:          /// On entry, if  COMPZ = 'V', then Z contains the orthogonal
 149:          /// matrix used in the reduction to tridiagonal form.
 150:          /// On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
 151:          /// orthonormal eigenvectors of the original symmetric matrix,
 152:          /// and if COMPZ = 'I', Z contains the orthonormal eigenvectors
 153:          /// of the symmetric tridiagonal matrix.
 154:          /// If COMPZ = 'N', then Z is not referenced.
 155:          ///</param>
 156:          /// <param name="LDZ">
 157:          /// (input) INTEGER
 158:          /// The leading dimension of the array Z.  LDZ .GE. 1, and if
 159:          /// eigenvectors are desired, then  LDZ .GE. max(1,N).
 160:          ///</param>
 161:          /// <param name="WORK">
 162:          /// (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
 163:          /// If COMPZ = 'N', then WORK is not referenced.
 164:          ///</param>
 165:          /// <param name="INFO">
 166:          /// (output) INTEGER
 167:          /// = 0:  successful exit
 168:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value
 169:          /// .GT. 0:  the algorithm has failed to find all the eigenvalues in
 170:          /// a total of 30*N iterations; if INFO = i, then i
 171:          /// elements of E have not converged to zero; on exit, D
 172:          /// and E contain the elements of a symmetric tridiagonal
 173:          /// matrix which is orthogonally similar to the original
 174:          /// matrix.
 175:          ///</param>
 176:          public void Run(string COMPZ, int N, ref double[] D, int offset_d, ref double[] E, int offset_e, ref double[] Z, int offset_z, int LDZ
 177:                           , ref double[] WORK, int offset_work, ref int INFO)
 178:          {
 179:   
 180:              #region Array Index Correction
 181:              
 182:               int o_d = -1 + offset_d;  int o_e = -1 + offset_e;  int o_z = -1 - LDZ + offset_z;  int o_work = -1 + offset_work; 
 183:   
 184:              #endregion
 185:   
 186:   
 187:              #region Strings
 188:              
 189:              COMPZ = COMPZ.Substring(0, 1);  
 190:   
 191:              #endregion
 192:   
 193:   
 194:              #region Prolog
 195:              
 196:              // *
 197:              // *  -- LAPACK routine (version 3.1) --
 198:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 199:              // *     November 2006
 200:              // *
 201:              // *     .. Scalar Arguments ..
 202:              // *     ..
 203:              // *     .. Array Arguments ..
 204:              // *     ..
 205:              // *
 206:              // *  Purpose
 207:              // *  =======
 208:              // *
 209:              // *  DSTEQR computes all eigenvalues and, optionally, eigenvectors of a
 210:              // *  symmetric tridiagonal matrix using the implicit QL or QR method.
 211:              // *  The eigenvectors of a full or band symmetric matrix can also be found
 212:              // *  if DSYTRD or DSPTRD or DSBTRD has been used to reduce this matrix to
 213:              // *  tridiagonal form.
 214:              // *
 215:              // *  Arguments
 216:              // *  =========
 217:              // *
 218:              // *  COMPZ   (input) CHARACTER*1
 219:              // *          = 'N':  Compute eigenvalues only.
 220:              // *          = 'V':  Compute eigenvalues and eigenvectors of the original
 221:              // *                  symmetric matrix.  On entry, Z must contain the
 222:              // *                  orthogonal matrix used to reduce the original matrix
 223:              // *                  to tridiagonal form.
 224:              // *          = 'I':  Compute eigenvalues and eigenvectors of the
 225:              // *                  tridiagonal matrix.  Z is initialized to the identity
 226:              // *                  matrix.
 227:              // *
 228:              // *  N       (input) INTEGER
 229:              // *          The order of the matrix.  N >= 0.
 230:              // *
 231:              // *  D       (input/output) DOUBLE PRECISION array, dimension (N)
 232:              // *          On entry, the diagonal elements of the tridiagonal matrix.
 233:              // *          On exit, if INFO = 0, the eigenvalues in ascending order.
 234:              // *
 235:              // *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
 236:              // *          On entry, the (n-1) subdiagonal elements of the tridiagonal
 237:              // *          matrix.
 238:              // *          On exit, E has been destroyed.
 239:              // *
 240:              // *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
 241:              // *          On entry, if  COMPZ = 'V', then Z contains the orthogonal
 242:              // *          matrix used in the reduction to tridiagonal form.
 243:              // *          On exit, if INFO = 0, then if  COMPZ = 'V', Z contains the
 244:              // *          orthonormal eigenvectors of the original symmetric matrix,
 245:              // *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
 246:              // *          of the symmetric tridiagonal matrix.
 247:              // *          If COMPZ = 'N', then Z is not referenced.
 248:              // *
 249:              // *  LDZ     (input) INTEGER
 250:              // *          The leading dimension of the array Z.  LDZ >= 1, and if
 251:              // *          eigenvectors are desired, then  LDZ >= max(1,N).
 252:              // *
 253:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (max(1,2*N-2))
 254:              // *          If COMPZ = 'N', then WORK is not referenced.
 255:              // *
 256:              // *  INFO    (output) INTEGER
 257:              // *          = 0:  successful exit
 258:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value
 259:              // *          > 0:  the algorithm has failed to find all the eigenvalues in
 260:              // *                a total of 30*N iterations; if INFO = i, then i
 261:              // *                elements of E have not converged to zero; on exit, D
 262:              // *                and E contain the elements of a symmetric tridiagonal
 263:              // *                matrix which is orthogonally similar to the original
 264:              // *                matrix.
 265:              // *
 266:              // *  =====================================================================
 267:              // *
 268:              // *     .. Parameters ..
 269:              // *     ..
 270:              // *     .. Local Scalars ..
 271:              // *     ..
 272:              // *     .. External Functions ..
 273:              // *     ..
 274:              // *     .. External Subroutines ..
 275:              // *     ..
 276:              // *     .. Intrinsic Functions ..
 277:              //      INTRINSIC          ABS, MAX, SIGN, SQRT;
 278:              // *     ..
 279:              // *     .. Executable Statements ..
 280:              // *
 281:              // *     Test the input parameters.
 282:              // *
 283:   
 284:              #endregion
 285:   
 286:   
 287:              #region Body
 288:              
 289:              INFO = 0;
 290:              // *
 291:              if (this._lsame.Run(COMPZ, "N"))
 292:              {
 293:                  ICOMPZ = 0;
 294:              }
 295:              else
 296:              {
 297:                  if (this._lsame.Run(COMPZ, "V"))
 298:                  {
 299:                      ICOMPZ = 1;
 300:                  }
 301:                  else
 302:                  {
 303:                      if (this._lsame.Run(COMPZ, "I"))
 304:                      {
 305:                          ICOMPZ = 2;
 306:                      }
 307:                      else
 308:                      {
 309:                          ICOMPZ =  - 1;
 310:                      }
 311:                  }
 312:              }
 313:              if (ICOMPZ < 0)
 314:              {
 315:                  INFO =  - 1;
 316:              }
 317:              else
 318:              {
 319:                  if (N < 0)
 320:                  {
 321:                      INFO =  - 2;
 322:                  }
 323:                  else
 324:                  {
 325:                      if ((LDZ < 1) || (ICOMPZ > 0 && LDZ < Math.Max(1, N)))
 326:                      {
 327:                          INFO =  - 6;
 328:                      }
 329:                  }
 330:              }
 331:              if (INFO != 0)
 332:              {
 333:                  this._xerbla.Run("DSTEQR",  - INFO);
 334:                  return;
 335:              }
 336:              // *
 337:              // *     Quick return if possible
 338:              // *
 339:              if (N == 0) return;
 340:              // *
 341:              if (N == 1)
 342:              {
 343:                  if (ICOMPZ == 2) Z[1+1 * LDZ + o_z] = ONE;
 344:                  return;
 345:              }
 346:              // *
 347:              // *     Determine the unit roundoff and over/underflow thresholds.
 348:              // *
 349:              EPS = this._dlamch.Run("E");
 350:              EPS2 = Math.Pow(EPS,2);
 351:              SAFMIN = this._dlamch.Run("S");
 352:              SAFMAX = ONE / SAFMIN;
 353:              SSFMAX = Math.Sqrt(SAFMAX) / THREE;
 354:              SSFMIN = Math.Sqrt(SAFMIN) / EPS2;
 355:              // *
 356:              // *     Compute the eigenvalues and eigenvectors of the tridiagonal
 357:              // *     matrix.
 358:              // *
 359:              if (ICOMPZ == 2)
 360:              {
 361:                  this._dlaset.Run("Full", N, N, ZERO, ONE, ref Z, offset_z
 362:                                   , LDZ);
 363:              }
 364:              // *
 365:              NMAXIT = N * MAXIT;
 366:              JTOT = 0;
 367:              // *
 368:              // *     Determine where the matrix splits and choose QL or QR iteration
 369:              // *     for each block, according to whether top or bottom diagonal
 370:              // *     element is smaller.
 371:              // *
 372:              L1 = 1;
 373:              NM1 = N - 1;
 374:              // *
 375:          LABEL10:;
 376:              if (L1 > N) goto LABEL160;
 377:              if (L1 > 1) E[L1 - 1 + o_e] = ZERO;
 378:              if (L1 <= NM1)
 379:              {
 380:                  for (M = L1; M <= NM1; M++)
 381:                  {
 382:                      TST = Math.Abs(E[M + o_e]);
 383:                      if (TST == ZERO) goto LABEL30;
 384:                      if (TST <= (Math.Sqrt(Math.Abs(D[M + o_d])) * Math.Sqrt(Math.Abs(D[M + 1 + o_d]))) * EPS)
 385:                      {
 386:                          E[M + o_e] = ZERO;
 387:                          goto LABEL30;
 388:                      }
 389:                  }
 390:              }
 391:              M = N;
 392:              // *
 393:          LABEL30:;
 394:              L = L1;
 395:              LSV = L;
 396:              LEND = M;
 397:              LENDSV = LEND;
 398:              L1 = M + 1;
 399:              if (LEND == L) goto LABEL10;
 400:              // *
 401:              // *     Scale submatrix in rows and columns L to LEND
 402:              // *
 403:              ANORM = this._dlanst.Run("I", LEND - L + 1, D, L + o_d, E, L + o_e);
 404:              ISCALE = 0;
 405:              if (ANORM == ZERO) goto LABEL10;
 406:              if (ANORM > SSFMAX)
 407:              {
 408:                  ISCALE = 1;
 409:                  this._dlascl.Run("G", 0, 0, ANORM, SSFMAX, LEND - L + 1
 410:                                   , 1, ref D, L + o_d, N, ref INFO);
 411:                  this._dlascl.Run("G", 0, 0, ANORM, SSFMAX, LEND - L
 412:                                   , 1, ref E, L + o_e, N, ref INFO);
 413:              }
 414:              else
 415:              {
 416:                  if (ANORM < SSFMIN)
 417:                  {
 418:                      ISCALE = 2;
 419:                      this._dlascl.Run("G", 0, 0, ANORM, SSFMIN, LEND - L + 1
 420:                                       , 1, ref D, L + o_d, N, ref INFO);
 421:                      this._dlascl.Run("G", 0, 0, ANORM, SSFMIN, LEND - L
 422:                                       , 1, ref E, L + o_e, N, ref INFO);
 423:                  }
 424:              }
 425:              // *
 426:              // *     Choose between QL and QR iteration
 427:              // *
 428:              if (Math.Abs(D[LEND + o_d]) < Math.Abs(D[L + o_d]))
 429:              {
 430:                  LEND = LSV;
 431:                  L = LENDSV;
 432:              }
 433:              // *
 434:              if (LEND > L)
 435:              {
 436:                  // *
 437:                  // *        QL Iteration
 438:                  // *
 439:                  // *        Look for small subdiagonal element.
 440:                  // *
 441:              LABEL40:;
 442:                  if (L != LEND)
 443:                  {
 444:                      LENDM1 = LEND - 1;
 445:                      for (M = L; M <= LENDM1; M++)
 446:                      {
 447:                          TST = Math.Pow(Math.Abs(E[M + o_e]),2);
 448:                          if (TST <= (EPS2 * Math.Abs(D[M + o_d])) * Math.Abs(D[M + 1 + o_d]) + SAFMIN) goto LABEL60;
 449:                      }
 450:                  }
 451:                  // *
 452:                  M = LEND;
 453:                  // *
 454:              LABEL60:;
 455:                  if (M < LEND) E[M + o_e] = ZERO;
 456:                  P = D[L + o_d];
 457:                  if (M == L) goto LABEL80;
 458:                  // *
 459:                  // *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
 460:                  // *        to compute its eigensystem.
 461:                  // *
 462:                  if (M == L + 1)
 463:                  {
 464:                      if (ICOMPZ > 0)
 465:                      {
 466:                          this._dlaev2.Run(D[L + o_d], E[L + o_e], D[L + 1 + o_d], ref RT1, ref RT2, ref C
 467:                                           , ref S);
 468:                          WORK[L + o_work] = C;
 469:                          WORK[N - 1 + L + o_work] = S;
 470:                          this._dlasr.Run("R", "V", "B", N, 2, WORK, L + o_work
 471:                                          , WORK, N - 1 + L + o_work, ref Z, 1+L * LDZ + o_z, LDZ);
 472:                      }
 473:                      else
 474:                      {
 475:                          this._dlae2.Run(D[L + o_d], E[L + o_e], D[L + 1 + o_d], ref RT1, ref RT2);
 476:                      }
 477:                      D[L + o_d] = RT1;
 478:                      D[L + 1 + o_d] = RT2;
 479:                      E[L + o_e] = ZERO;
 480:                      L = L + 2;
 481:                      if (L <= LEND) goto LABEL40;
 482:                      goto LABEL140;
 483:                  }
 484:                  // *
 485:                  if (JTOT == NMAXIT) goto LABEL140;
 486:                  JTOT = JTOT + 1;
 487:                  // *
 488:                  // *        Form shift.
 489:                  // *
 490:                  G = (D[L + 1 + o_d] - P) / (TWO * E[L + o_e]);
 491:                  R = this._dlapy2.Run(G, ONE);
 492:                  G = D[M + o_d] - P + (E[L + o_e] / (G + FortranLib.Sign(R,G)));
 493:                  // *
 494:                  S = ONE;
 495:                  C = ONE;
 496:                  P = ZERO;
 497:                  // *
 498:                  // *        Inner loop
 499:                  // *
 500:                  MM1 = M - 1;
 501:                  for (I = MM1; I >= L; I +=  - 1)
 502:                  {
 503:                      F = S * E[I + o_e];
 504:                      B = C * E[I + o_e];
 505:                      this._dlartg.Run(G, F, ref C, ref S, ref R);
 506:                      if (I != M - 1) E[I + 1 + o_e] = R;
 507:                      G = D[I + 1 + o_d] - P;
 508:                      R = (D[I + o_d] - G) * S + TWO * C * B;
 509:                      P = S * R;
 510:                      D[I + 1 + o_d] = G + P;
 511:                      G = C * R - B;
 512:                      // *
 513:                      // *           If eigenvectors are desired, then save rotations.
 514:                      // *
 515:                      if (ICOMPZ > 0)
 516:                      {
 517:                          WORK[I + o_work] = C;
 518:                          WORK[N - 1 + I + o_work] =  - S;
 519:                      }
 520:                      // *
 521:                  }
 522:                  // *
 523:                  // *        If eigenvectors are desired, then apply saved rotations.
 524:                  // *
 525:                  if (ICOMPZ > 0)
 526:                  {
 527:                      MM = M - L + 1;
 528:                      this._dlasr.Run("R", "V", "B", N, MM, WORK, L + o_work
 529:                                      , WORK, N - 1 + L + o_work, ref Z, 1+L * LDZ + o_z, LDZ);
 530:                  }
 531:                  // *
 532:                  D[L + o_d] = D[L + o_d] - P;
 533:                  E[L + o_e] = G;
 534:                  goto LABEL40;
 535:                  // *
 536:                  // *        Eigenvalue found.
 537:                  // *
 538:              LABEL80:;
 539:                  D[L + o_d] = P;
 540:                  // *
 541:                  L = L + 1;
 542:                  if (L <= LEND) goto LABEL40;
 543:                  goto LABEL140;
 544:                  // *
 545:              }
 546:              else
 547:              {
 548:                  // *
 549:                  // *        QR Iteration
 550:                  // *
 551:                  // *        Look for small superdiagonal element.
 552:                  // *
 553:              LABEL90:;
 554:                  if (L != LEND)
 555:                  {
 556:                      LENDP1 = LEND + 1;
 557:                      for (M = L; M >= LENDP1; M +=  - 1)
 558:                      {
 559:                          TST = Math.Pow(Math.Abs(E[M - 1 + o_e]),2);
 560:                          if (TST <= (EPS2 * Math.Abs(D[M + o_d])) * Math.Abs(D[M - 1 + o_d]) + SAFMIN) goto LABEL110;
 561:                      }
 562:                  }
 563:                  // *
 564:                  M = LEND;
 565:                  // *
 566:              LABEL110:;
 567:                  if (M > LEND) E[M - 1 + o_e] = ZERO;
 568:                  P = D[L + o_d];
 569:                  if (M == L) goto LABEL130;
 570:                  // *
 571:                  // *        If remaining matrix is 2-by-2, use DLAE2 or SLAEV2
 572:                  // *        to compute its eigensystem.
 573:                  // *
 574:                  if (M == L - 1)
 575:                  {
 576:                      if (ICOMPZ > 0)
 577:                      {
 578:                          this._dlaev2.Run(D[L - 1 + o_d], E[L - 1 + o_e], D[L + o_d], ref RT1, ref RT2, ref C
 579:                                           , ref S);
 580:                          WORK[M + o_work] = C;
 581:                          WORK[N - 1 + M + o_work] = S;
 582:                          this._dlasr.Run("R", "V", "F", N, 2, WORK, M + o_work
 583:                                          , WORK, N - 1 + M + o_work, ref Z, 1+(L - 1) * LDZ + o_z, LDZ);
 584:                      }
 585:                      else
 586:                      {
 587:                          this._dlae2.Run(D[L - 1 + o_d], E[L - 1 + o_e], D[L + o_d], ref RT1, ref RT2);
 588:                      }
 589:                      D[L - 1 + o_d] = RT1;
 590:                      D[L + o_d] = RT2;
 591:                      E[L - 1 + o_e] = ZERO;
 592:                      L = L - 2;
 593:                      if (L >= LEND) goto LABEL90;
 594:                      goto LABEL140;
 595:                  }
 596:                  // *
 597:                  if (JTOT == NMAXIT) goto LABEL140;
 598:                  JTOT = JTOT + 1;
 599:                  // *
 600:                  // *        Form shift.
 601:                  // *
 602:                  G = (D[L - 1 + o_d] - P) / (TWO * E[L - 1 + o_e]);
 603:                  R = this._dlapy2.Run(G, ONE);
 604:                  G = D[M + o_d] - P + (E[L - 1 + o_e] / (G + FortranLib.Sign(R,G)));
 605:                  // *
 606:                  S = ONE;
 607:                  C = ONE;
 608:                  P = ZERO;
 609:                  // *
 610:                  // *        Inner loop
 611:                  // *
 612:                  LM1 = L - 1;
 613:                  for (I = M; I <= LM1; I++)
 614:                  {
 615:                      F = S * E[I + o_e];
 616:                      B = C * E[I + o_e];
 617:                      this._dlartg.Run(G, F, ref C, ref S, ref R);
 618:                      if (I != M) E[I - 1 + o_e] = R;
 619:                      G = D[I + o_d] - P;
 620:                      R = (D[I + 1 + o_d] - G) * S + TWO * C * B;
 621:                      P = S * R;
 622:                      D[I + o_d] = G + P;
 623:                      G = C * R - B;
 624:                      // *
 625:                      // *           If eigenvectors are desired, then save rotations.
 626:                      // *
 627:                      if (ICOMPZ > 0)
 628:                      {
 629:                          WORK[I + o_work] = C;
 630:                          WORK[N - 1 + I + o_work] = S;
 631:                      }
 632:                      // *
 633:                  }
 634:                  // *
 635:                  // *        If eigenvectors are desired, then apply saved rotations.
 636:                  // *
 637:                  if (ICOMPZ > 0)
 638:                  {
 639:                      MM = L - M + 1;
 640:                      this._dlasr.Run("R", "V", "F", N, MM, WORK, M + o_work
 641:                                      , WORK, N - 1 + M + o_work, ref Z, 1+M * LDZ + o_z, LDZ);
 642:                  }
 643:                  // *
 644:                  D[L + o_d] = D[L + o_d] - P;
 645:                  E[LM1 + o_e] = G;
 646:                  goto LABEL90;
 647:                  // *
 648:                  // *        Eigenvalue found.
 649:                  // *
 650:              LABEL130:;
 651:                  D[L + o_d] = P;
 652:                  // *
 653:                  L = L - 1;
 654:                  if (L >= LEND) goto LABEL90;
 655:                  goto LABEL140;
 656:                  // *
 657:              }
 658:              // *
 659:              // *     Undo scaling if necessary
 660:              // *
 661:          LABEL140:;
 662:              if (ISCALE == 1)
 663:              {
 664:                  this._dlascl.Run("G", 0, 0, SSFMAX, ANORM, LENDSV - LSV + 1
 665:                                   , 1, ref D, LSV + o_d, N, ref INFO);
 666:                  this._dlascl.Run("G", 0, 0, SSFMAX, ANORM, LENDSV - LSV
 667:                                   , 1, ref E, LSV + o_e, N, ref INFO);
 668:              }
 669:              else
 670:              {
 671:                  if (ISCALE == 2)
 672:                  {
 673:                      this._dlascl.Run("G", 0, 0, SSFMIN, ANORM, LENDSV - LSV + 1
 674:                                       , 1, ref D, LSV + o_d, N, ref INFO);
 675:                      this._dlascl.Run("G", 0, 0, SSFMIN, ANORM, LENDSV - LSV
 676:                                       , 1, ref E, LSV + o_e, N, ref INFO);
 677:                  }
 678:              }
 679:              // *
 680:              // *     Check for no convergence to an eigenvalue after a total
 681:              // *     of N*MAXIT iterations.
 682:              // *
 683:              if (JTOT < NMAXIT) goto LABEL10;
 684:              for (I = 1; I <= N - 1; I++)
 685:              {
 686:                  if (E[I + o_e] != ZERO) INFO = INFO + 1;
 687:              }
 688:              goto LABEL190;
 689:              // *
 690:              // *     Order eigenvalues and eigenvectors.
 691:              // *
 692:          LABEL160:;
 693:              if (ICOMPZ == 0)
 694:              {
 695:                  // *
 696:                  // *        Use Quick Sort
 697:                  // *
 698:                  this._dlasrt.Run("I", N, ref D, offset_d, ref INFO);
 699:                  // *
 700:              }
 701:              else
 702:              {
 703:                  // *
 704:                  // *        Use Selection Sort to minimize swaps of eigenvectors
 705:                  // *
 706:                  for (II = 2; II <= N; II++)
 707:                  {
 708:                      I = II - 1;
 709:                      K = I;
 710:                      P = D[I + o_d];
 711:                      for (J = II; J <= N; J++)
 712:                      {
 713:                          if (D[J + o_d] < P)
 714:                          {
 715:                              K = J;
 716:                              P = D[J + o_d];
 717:                          }
 718:                      }
 719:                      if (K != I)
 720:                      {
 721:                          D[K + o_d] = D[I + o_d];
 722:                          D[I + o_d] = P;
 723:                          this._dswap.Run(N, ref Z, 1+I * LDZ + o_z, 1, ref Z, 1+K * LDZ + o_z, 1);
 724:                      }
 725:                  }
 726:              }
 727:              // *
 728:          LABEL190:;
 729:              return;
 730:              // *
 731:              // *     End of DSTEQR
 732:              // *
 733:   
 734:              #endregion
 735:   
 736:          }
 737:      }
 738:  }