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 driver routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
  27:      /// symmetric tridiagonal matrix using the divide and conquer method.
  28:      /// The eigenvectors of a full or band real symmetric matrix can also be
  29:      /// found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
  30:      /// matrix to tridiagonal form.
  31:      /// 
  32:      /// This code makes very mild assumptions about floating point
  33:      /// arithmetic. It will work on machines with a guard digit in
  34:      /// add/subtract, or on those binary machines without guard digits
  35:      /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
  36:      /// It could conceivably fail on hexadecimal or decimal machines
  37:      /// without guard digits, but we know of none.  See DLAED3 for details.
  38:      /// 
  39:      ///</summary>
  40:      public class DSTEDC
  41:      {
  42:      
  43:   
  44:          #region Dependencies
  45:          
  46:          LSAME _lsame; ILAENV _ilaenv; DLAMCH _dlamch; DLANST _dlanst; DGEMM _dgemm; DLACPY _dlacpy; DLAED0 _dlaed0; 
  47:          DLASCL _dlascl;DLASET _dlaset; DLASRT _dlasrt; DSTEQR _dsteqr; DSTERF _dsterf; DSWAP _dswap; XERBLA _xerbla; 
  48:   
  49:          #endregion
  50:   
  51:   
  52:          #region Fields
  53:          
  54:          const double ZERO = 0.0E0; const double ONE = 1.0E0; const double TWO = 2.0E0; bool LQUERY = false; int FINISH = 0; 
  55:          int I = 0;int ICOMPZ = 0; int II = 0; int J = 0; int K = 0; int LGN = 0; int LIWMIN = 0; int LWMIN = 0; int M = 0; 
  56:          int SMLSIZ = 0;int START = 0; int STOREZ = 0; int STRTRW = 0; double EPS = 0; double ORGNRM = 0; double P = 0; 
  57:          double TINY = 0;
  58:   
  59:          #endregion
  60:   
  61:          public DSTEDC(LSAME lsame, ILAENV ilaenv, DLAMCH dlamch, DLANST dlanst, DGEMM dgemm, DLACPY dlacpy, DLAED0 dlaed0, DLASCL dlascl, DLASET dlaset, DLASRT dlasrt
  62:                        , DSTEQR dsteqr, DSTERF dsterf, DSWAP dswap, XERBLA xerbla)
  63:          {
  64:      
  65:   
  66:              #region Set Dependencies
  67:              
  68:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlanst = dlanst; this._dgemm = dgemm; 
  69:              this._dlacpy = dlacpy;this._dlaed0 = dlaed0; this._dlascl = dlascl; this._dlaset = dlaset; this._dlasrt = dlasrt; 
  70:              this._dsteqr = dsteqr;this._dsterf = dsterf; this._dswap = dswap; this._xerbla = xerbla; 
  71:   
  72:              #endregion
  73:   
  74:          }
  75:      
  76:          public DSTEDC()
  77:          {
  78:      
  79:   
  80:              #region Dependencies (Initialization)
  81:              
  82:              LSAME lsame = new LSAME();
  83:              IEEECK ieeeck = new IEEECK();
  84:              IPARMQ iparmq = new IPARMQ();
  85:              DLAMC3 dlamc3 = new DLAMC3();
  86:              DLASSQ dlassq = new DLASSQ();
  87:              XERBLA xerbla = new XERBLA();
  88:              DCOPY dcopy = new DCOPY();
  89:              IDAMAX idamax = new IDAMAX();
  90:              DLAPY2 dlapy2 = new DLAPY2();
  91:              DLAMRG dlamrg = new DLAMRG();
  92:              DROT drot = new DROT();
  93:              DSCAL dscal = new DSCAL();
  94:              DNRM2 dnrm2 = new DNRM2();
  95:              DLAED5 dlaed5 = new DLAED5();
  96:              DLAE2 dlae2 = new DLAE2();
  97:              DLAEV2 dlaev2 = new DLAEV2();
  98:              DSWAP dswap = new DSWAP();
  99:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
 100:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
 101:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 102:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 103:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 104:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 105:              DLANST dlanst = new DLANST(lsame, dlassq);
 106:              DGEMM dgemm = new DGEMM(lsame, xerbla);
 107:              DLACPY dlacpy = new DLACPY(lsame);
 108:              DLAED2 dlaed2 = new DLAED2(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
 109:              DLAED6 dlaed6 = new DLAED6(dlamch);
 110:              DLAED4 dlaed4 = new DLAED4(dlamch, dlaed5, dlaed6);
 111:              DLASET dlaset = new DLASET(lsame);
 112:              DLAED3 dlaed3 = new DLAED3(dlamc3, dnrm2, dcopy, dgemm, dlacpy, dlaed4, dlaset, xerbla);
 113:              DLAED1 dlaed1 = new DLAED1(dcopy, dlaed2, dlaed3, dlamrg, xerbla);
 114:              DLAED8 dlaed8 = new DLAED8(idamax, dlamch, dlapy2, dcopy, dlacpy, dlamrg, drot, dscal, xerbla);
 115:              DLAED9 dlaed9 = new DLAED9(dlamc3, dnrm2, dcopy, dlaed4, xerbla);
 116:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 117:              DLAEDA dlaeda = new DLAEDA(dcopy, dgemv, drot, xerbla);
 118:              DLAED7 dlaed7 = new DLAED7(dgemm, dlaed8, dlaed9, dlaeda, dlamrg, xerbla);
 119:              DLARTG dlartg = new DLARTG(dlamch);
 120:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
 121:              DLASR dlasr = new DLASR(lsame, xerbla);
 122:              DLASRT dlasrt = new DLASRT(lsame, xerbla);
 123:              DSTEQR dsteqr = new DSTEQR(lsame, dlamch, dlanst, dlapy2, dlae2, dlaev2, dlartg, dlascl, dlaset, dlasr
 124:                                         , dlasrt, dswap, xerbla);
 125:              DLAED0 dlaed0 = new DLAED0(dcopy, dgemm, dlacpy, dlaed1, dlaed7, dsteqr, xerbla, ilaenv);
 126:              DSTERF dsterf = new DSTERF(dlamch, dlanst, dlapy2, dlae2, dlascl, dlasrt, xerbla);
 127:   
 128:              #endregion
 129:   
 130:   
 131:              #region Set Dependencies
 132:              
 133:              this._lsame = lsame; this._ilaenv = ilaenv; this._dlamch = dlamch; this._dlanst = dlanst; this._dgemm = dgemm; 
 134:              this._dlacpy = dlacpy;this._dlaed0 = dlaed0; this._dlascl = dlascl; this._dlaset = dlaset; this._dlasrt = dlasrt; 
 135:              this._dsteqr = dsteqr;this._dsterf = dsterf; this._dswap = dswap; this._xerbla = xerbla; 
 136:   
 137:              #endregion
 138:   
 139:          }
 140:          /// <summary>
 141:          /// Purpose
 142:          /// =======
 143:          /// 
 144:          /// DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
 145:          /// symmetric tridiagonal matrix using the divide and conquer method.
 146:          /// The eigenvectors of a full or band real symmetric matrix can also be
 147:          /// found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
 148:          /// matrix to tridiagonal form.
 149:          /// 
 150:          /// This code makes very mild assumptions about floating point
 151:          /// arithmetic. It will work on machines with a guard digit in
 152:          /// add/subtract, or on those binary machines without guard digits
 153:          /// which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
 154:          /// It could conceivably fail on hexadecimal or decimal machines
 155:          /// without guard digits, but we know of none.  See DLAED3 for details.
 156:          /// 
 157:          ///</summary>
 158:          /// <param name="COMPZ">
 159:          /// (input) CHARACTER*1
 160:          /// = 'N':  Compute eigenvalues only.
 161:          /// = 'I':  Compute eigenvectors of tridiagonal matrix also.
 162:          /// = 'V':  Compute eigenvectors of original dense symmetric
 163:          /// matrix also.  On entry, Z contains the orthogonal
 164:          /// matrix used to reduce the original matrix to
 165:          /// tridiagonal form.
 166:          ///</param>
 167:          /// <param name="N">
 168:          /// (input) INTEGER
 169:          /// The dimension of the symmetric tridiagonal matrix.  N .GE. 0.
 170:          ///</param>
 171:          /// <param name="D">
 172:          /// (input/output) DOUBLE PRECISION array, dimension (N)
 173:          /// On entry, the diagonal elements of the tridiagonal matrix.
 174:          /// On exit, if INFO = 0, the eigenvalues in ascending order.
 175:          ///</param>
 176:          /// <param name="E">
 177:          /// (input/output) DOUBLE PRECISION array, dimension (N-1)
 178:          /// On entry, the subdiagonal elements of the tridiagonal matrix.
 179:          /// On exit, E has been destroyed.
 180:          ///</param>
 181:          /// <param name="Z">
 182:          /// (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
 183:          /// On entry, if COMPZ = 'V', then Z contains the orthogonal
 184:          /// matrix used in the reduction to tridiagonal form.
 185:          /// On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
 186:          /// orthonormal eigenvectors of the original symmetric matrix,
 187:          /// and if COMPZ = 'I', Z contains the orthonormal eigenvectors
 188:          /// of the symmetric tridiagonal matrix.
 189:          /// If  COMPZ = 'N', then Z is not referenced.
 190:          ///</param>
 191:          /// <param name="LDZ">
 192:          /// (input) INTEGER
 193:          /// The leading dimension of the array Z.  LDZ .GE. 1.
 194:          /// If eigenvectors are desired, then LDZ .GE. max(1,N).
 195:          ///</param>
 196:          /// <param name="WORK">
 197:          /// (workspace/output) DOUBLE PRECISION array,
 198:          /// dimension (LWORK)
 199:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 200:          ///</param>
 201:          /// <param name="LWORK">
 202:          /// (input) INTEGER
 203:          /// The dimension of the array WORK.
 204:          /// If COMPZ = 'N' or N .LE. 1 then LWORK must be at least 1.
 205:          /// If COMPZ = 'V' and N .GT. 1 then LWORK must be at least
 206:          /// ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
 207:          /// where lg( N ) = smallest integer k such
 208:          /// that 2**k .GE. N.
 209:          /// If COMPZ = 'I' and N .GT. 1 then LWORK must be at least
 210:          /// ( 1 + 4*N + N**2 ).
 211:          /// Note that for COMPZ = 'I' or 'V', then if N is less than or
 212:          /// equal to the minimum divide size, usually 25, then LWORK need
 213:          /// only be max(1,2*(N-1)).
 214:          /// 
 215:          /// If LWORK = -1, then a workspace query is assumed; the routine
 216:          /// only calculates the optimal size of the WORK array, returns
 217:          /// this value as the first entry of the WORK array, and no error
 218:          /// message related to LWORK is issued by XERBLA.
 219:          ///</param>
 220:          /// <param name="IWORK">
 221:          /// (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
 222:          /// On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 223:          ///</param>
 224:          /// <param name="LIWORK">
 225:          /// (input) INTEGER
 226:          /// The dimension of the array IWORK.
 227:          /// If COMPZ = 'N' or N .LE. 1 then LIWORK must be at least 1.
 228:          /// If COMPZ = 'V' and N .GT. 1 then LIWORK must be at least
 229:          /// ( 6 + 6*N + 5*N*lg N ).
 230:          /// If COMPZ = 'I' and N .GT. 1 then LIWORK must be at least
 231:          /// ( 3 + 5*N ).
 232:          /// Note that for COMPZ = 'I' or 'V', then if N is less than or
 233:          /// equal to the minimum divide size, usually 25, then LIWORK
 234:          /// need only be 1.
 235:          /// 
 236:          /// If LIWORK = -1, then a workspace query is assumed; the
 237:          /// routine only calculates the optimal size of the IWORK array,
 238:          /// returns this value as the first entry of the IWORK array, and
 239:          /// no error message related to LIWORK is issued by XERBLA.
 240:          ///</param>
 241:          /// <param name="INFO">
 242:          /// (output) INTEGER
 243:          /// = 0:  successful exit.
 244:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 245:          /// .GT. 0:  The algorithm failed to compute an eigenvalue while
 246:          /// working on the submatrix lying in rows and columns
 247:          /// INFO/(N+1) through mod(INFO,N+1).
 248:          ///</param>
 249:          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
 250:                           , ref double[] WORK, int offset_work, int LWORK, ref int[] IWORK, int offset_iwork, int LIWORK, ref int INFO)
 251:          {
 252:   
 253:              #region Array Index Correction
 254:              
 255:               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; 
 256:               int o_iwork = -1 + offset_iwork;
 257:   
 258:              #endregion
 259:   
 260:   
 261:              #region Strings
 262:              
 263:              COMPZ = COMPZ.Substring(0, 1);  
 264:   
 265:              #endregion
 266:   
 267:   
 268:              #region Prolog
 269:              
 270:              // *
 271:              // *  -- LAPACK driver routine (version 3.1) --
 272:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 273:              // *     November 2006
 274:              // *
 275:              // *     .. Scalar Arguments ..
 276:              // *     ..
 277:              // *     .. Array Arguments ..
 278:              // *     ..
 279:              // *
 280:              // *  Purpose
 281:              // *  =======
 282:              // *
 283:              // *  DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
 284:              // *  symmetric tridiagonal matrix using the divide and conquer method.
 285:              // *  The eigenvectors of a full or band real symmetric matrix can also be
 286:              // *  found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
 287:              // *  matrix to tridiagonal form.
 288:              // *
 289:              // *  This code makes very mild assumptions about floating point
 290:              // *  arithmetic. It will work on machines with a guard digit in
 291:              // *  add/subtract, or on those binary machines without guard digits
 292:              // *  which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
 293:              // *  It could conceivably fail on hexadecimal or decimal machines
 294:              // *  without guard digits, but we know of none.  See DLAED3 for details.
 295:              // *
 296:              // *  Arguments
 297:              // *  =========
 298:              // *
 299:              // *  COMPZ   (input) CHARACTER*1
 300:              // *          = 'N':  Compute eigenvalues only.
 301:              // *          = 'I':  Compute eigenvectors of tridiagonal matrix also.
 302:              // *          = 'V':  Compute eigenvectors of original dense symmetric
 303:              // *                  matrix also.  On entry, Z contains the orthogonal
 304:              // *                  matrix used to reduce the original matrix to
 305:              // *                  tridiagonal form.
 306:              // *
 307:              // *  N       (input) INTEGER
 308:              // *          The dimension of the symmetric tridiagonal matrix.  N >= 0.
 309:              // *
 310:              // *  D       (input/output) DOUBLE PRECISION array, dimension (N)
 311:              // *          On entry, the diagonal elements of the tridiagonal matrix.
 312:              // *          On exit, if INFO = 0, the eigenvalues in ascending order.
 313:              // *
 314:              // *  E       (input/output) DOUBLE PRECISION array, dimension (N-1)
 315:              // *          On entry, the subdiagonal elements of the tridiagonal matrix.
 316:              // *          On exit, E has been destroyed.
 317:              // *
 318:              // *  Z       (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
 319:              // *          On entry, if COMPZ = 'V', then Z contains the orthogonal
 320:              // *          matrix used in the reduction to tridiagonal form.
 321:              // *          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
 322:              // *          orthonormal eigenvectors of the original symmetric matrix,
 323:              // *          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
 324:              // *          of the symmetric tridiagonal matrix.
 325:              // *          If  COMPZ = 'N', then Z is not referenced.
 326:              // *
 327:              // *  LDZ     (input) INTEGER
 328:              // *          The leading dimension of the array Z.  LDZ >= 1.
 329:              // *          If eigenvectors are desired, then LDZ >= max(1,N).
 330:              // *
 331:              // *  WORK    (workspace/output) DOUBLE PRECISION array,
 332:              // *                                         dimension (LWORK)
 333:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 334:              // *
 335:              // *  LWORK   (input) INTEGER
 336:              // *          The dimension of the array WORK.
 337:              // *          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
 338:              // *          If COMPZ = 'V' and N > 1 then LWORK must be at least
 339:              // *                         ( 1 + 3*N + 2*N*lg N + 3*N**2 ),
 340:              // *                         where lg( N ) = smallest integer k such
 341:              // *                         that 2**k >= N.
 342:              // *          If COMPZ = 'I' and N > 1 then LWORK must be at least
 343:              // *                         ( 1 + 4*N + N**2 ).
 344:              // *          Note that for COMPZ = 'I' or 'V', then if N is less than or
 345:              // *          equal to the minimum divide size, usually 25, then LWORK need
 346:              // *          only be max(1,2*(N-1)).
 347:              // *
 348:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 349:              // *          only calculates the optimal size of the WORK array, returns
 350:              // *          this value as the first entry of the WORK array, and no error
 351:              // *          message related to LWORK is issued by XERBLA.
 352:              // *
 353:              // *  IWORK   (workspace/output) INTEGER array, dimension (MAX(1,LIWORK))
 354:              // *          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 355:              // *
 356:              // *  LIWORK  (input) INTEGER
 357:              // *          The dimension of the array IWORK.
 358:              // *          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
 359:              // *          If COMPZ = 'V' and N > 1 then LIWORK must be at least
 360:              // *                         ( 6 + 6*N + 5*N*lg N ).
 361:              // *          If COMPZ = 'I' and N > 1 then LIWORK must be at least
 362:              // *                         ( 3 + 5*N ).
 363:              // *          Note that for COMPZ = 'I' or 'V', then if N is less than or
 364:              // *          equal to the minimum divide size, usually 25, then LIWORK
 365:              // *          need only be 1.
 366:              // *
 367:              // *          If LIWORK = -1, then a workspace query is assumed; the
 368:              // *          routine only calculates the optimal size of the IWORK array,
 369:              // *          returns this value as the first entry of the IWORK array, and
 370:              // *          no error message related to LIWORK is issued by XERBLA.
 371:              // *
 372:              // *  INFO    (output) INTEGER
 373:              // *          = 0:  successful exit.
 374:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 375:              // *          > 0:  The algorithm failed to compute an eigenvalue while
 376:              // *                working on the submatrix lying in rows and columns
 377:              // *                INFO/(N+1) through mod(INFO,N+1).
 378:              // *
 379:              // *  Further Details
 380:              // *  ===============
 381:              // *
 382:              // *  Based on contributions by
 383:              // *     Jeff Rutter, Computer Science Division, University of California
 384:              // *     at Berkeley, USA
 385:              // *  Modified by Francoise Tisseur, University of Tennessee.
 386:              // *
 387:              // *  =====================================================================
 388:              // *
 389:              // *     .. Parameters ..
 390:              // *     ..
 391:              // *     .. Local Scalars ..
 392:              // *     ..
 393:              // *     .. External Functions ..
 394:              // *     ..
 395:              // *     .. External Subroutines ..
 396:              // *     ..
 397:              // *     .. Intrinsic Functions ..
 398:              //      INTRINSIC          ABS, DBLE, INT, LOG, MAX, MOD, SQRT;
 399:              // *     ..
 400:              // *     .. Executable Statements ..
 401:              // *
 402:              // *     Test the input parameters.
 403:              // *
 404:   
 405:              #endregion
 406:   
 407:   
 408:              #region Body
 409:              
 410:              INFO = 0;
 411:              LQUERY = (LWORK ==  - 1 || LIWORK ==  - 1);
 412:              // *
 413:              if (this._lsame.Run(COMPZ, "N"))
 414:              {
 415:                  ICOMPZ = 0;
 416:              }
 417:              else
 418:              {
 419:                  if (this._lsame.Run(COMPZ, "V"))
 420:                  {
 421:                      ICOMPZ = 1;
 422:                  }
 423:                  else
 424:                  {
 425:                      if (this._lsame.Run(COMPZ, "I"))
 426:                      {
 427:                          ICOMPZ = 2;
 428:                      }
 429:                      else
 430:                      {
 431:                          ICOMPZ =  - 1;
 432:                      }
 433:                  }
 434:              }
 435:              if (ICOMPZ < 0)
 436:              {
 437:                  INFO =  - 1;
 438:              }
 439:              else
 440:              {
 441:                  if (N < 0)
 442:                  {
 443:                      INFO =  - 2;
 444:                  }
 445:                  else
 446:                  {
 447:                      if ((LDZ < 1) || (ICOMPZ > 0 && LDZ < Math.Max(1, N)))
 448:                      {
 449:                          INFO =  - 6;
 450:                      }
 451:                  }
 452:              }
 453:              // *
 454:              if (INFO == 0)
 455:              {
 456:                  // *
 457:                  // *        Compute the workspace requirements
 458:                  // *
 459:                  SMLSIZ = this._ilaenv.Run(9, "DSTEDC", " ", 0, 0, 0, 0);
 460:                  if (N <= 1 || ICOMPZ == 0)
 461:                  {
 462:                      LIWMIN = 1;
 463:                      LWMIN = 1;
 464:                  }
 465:                  else
 466:                  {
 467:                      if (N <= SMLSIZ)
 468:                      {
 469:                          LIWMIN = 1;
 470:                          LWMIN = 2 * (N - 1);
 471:                      }
 472:                      else
 473:                      {
 474:                          LGN = Convert.ToInt32(Math.Truncate(Math.Log(Convert.ToDouble(N)) / Math.Log(TWO)));
 475:                          if (Math.Pow(2,LGN) < N) LGN = LGN + 1;
 476:                          if (Math.Pow(2,LGN) < N) LGN = LGN + 1;
 477:                          if (ICOMPZ == 1)
 478:                          {
 479:                              LWMIN = 1 + 3 * N + 2 * N * LGN + 3 * (int)Math.Pow(N, 2);
 480:                              LIWMIN = 6 + 6 * N + 5 * N * LGN;
 481:                          }
 482:                          else
 483:                          {
 484:                              if (ICOMPZ == 2)
 485:                              {
 486:                                  LWMIN = 1 + 4 * N + (int)Math.Pow(N, 2);
 487:                                  LIWMIN = 3 + 5 * N;
 488:                              }
 489:                          }
 490:                      }
 491:                  }
 492:                  WORK[1 + o_work] = LWMIN;
 493:                  IWORK[1 + o_iwork] = LIWMIN;
 494:                  // *
 495:                  if (LWORK < LWMIN && !LQUERY)
 496:                  {
 497:                      INFO =  - 8;
 498:                  }
 499:                  else
 500:                  {
 501:                      if (LIWORK < LIWMIN && !LQUERY)
 502:                      {
 503:                          INFO =  - 10;
 504:                      }
 505:                  }
 506:              }
 507:              // *
 508:              if (INFO != 0)
 509:              {
 510:                  this._xerbla.Run("DSTEDC",  - INFO);
 511:                  return;
 512:              }
 513:              else
 514:              {
 515:                  if (LQUERY)
 516:                  {
 517:                      return;
 518:                  }
 519:              }
 520:              // *
 521:              // *     Quick return if possible
 522:              // *
 523:              if (N == 0) return;
 524:              if (N == 1)
 525:              {
 526:                  if (ICOMPZ != 0) Z[1+1 * LDZ + o_z] = ONE;
 527:                  return;
 528:              }
 529:              // *
 530:              // *     If the following conditional clause is removed, then the routine
 531:              // *     will use the Divide and Conquer routine to compute only the
 532:              // *     eigenvalues, which requires (3N + 3N**2) real workspace and
 533:              // *     (2 + 5N + 2N lg(N)) integer workspace.
 534:              // *     Since on many architectures DSTERF is much faster than any other
 535:              // *     algorithm for finding eigenvalues only, it is used here
 536:              // *     as the default. If the conditional clause is removed, then
 537:              // *     information on the size of workspace needs to be changed.
 538:              // *
 539:              // *     If COMPZ = 'N', use DSTERF to compute the eigenvalues.
 540:              // *
 541:              if (ICOMPZ == 0)
 542:              {
 543:                  this._dsterf.Run(N, ref D, offset_d, ref E, offset_e, ref INFO);
 544:                  goto LABEL50;
 545:              }
 546:              // *
 547:              // *     If N is smaller than the minimum divide size (SMLSIZ+1), then
 548:              // *     solve the problem with another solver.
 549:              // *
 550:              if (N <= SMLSIZ)
 551:              {
 552:                  // *
 553:                  this._dsteqr.Run(COMPZ, N, ref D, offset_d, ref E, offset_e, ref Z, offset_z, LDZ
 554:                                   , ref WORK, offset_work, ref INFO);
 555:                  // *
 556:              }
 557:              else
 558:              {
 559:                  // *
 560:                  // *        If COMPZ = 'V', the Z matrix must be stored elsewhere for later
 561:                  // *        use.
 562:                  // *
 563:                  if (ICOMPZ == 1)
 564:                  {
 565:                      STOREZ = 1 + N * N;
 566:                  }
 567:                  else
 568:                  {
 569:                      STOREZ = 1;
 570:                  }
 571:                  // *
 572:                  if (ICOMPZ == 2)
 573:                  {
 574:                      this._dlaset.Run("Full", N, N, ZERO, ONE, ref Z, offset_z
 575:                                       , LDZ);
 576:                  }
 577:                  // *
 578:                  // *        Scale.
 579:                  // *
 580:                  ORGNRM = this._dlanst.Run("M", N, D, offset_d, E, offset_e);
 581:                  if (ORGNRM == ZERO) goto LABEL50;
 582:                  // *
 583:                  EPS = this._dlamch.Run("Epsilon");
 584:                  // *
 585:                  START = 1;
 586:                  // *
 587:                  // *        while ( START <= N )
 588:                  // *
 589:              LABEL10:;
 590:                  if (START <= N)
 591:                  {
 592:                      // *
 593:                      // *           Let FINISH be the position of the next subdiagonal entry
 594:                      // *           such that E( FINISH ) <= TINY or FINISH = N if no such
 595:                      // *           subdiagonal exists.  The matrix identified by the elements
 596:                      // *           between START and FINISH constitutes an independent
 597:                      // *           sub-problem.
 598:                      // *
 599:                      FINISH = START;
 600:                  LABEL20:;
 601:                      if (FINISH < N)
 602:                      {
 603:                          TINY = EPS * Math.Sqrt(Math.Abs(D[FINISH + o_d])) * Math.Sqrt(Math.Abs(D[FINISH + 1 + o_d]));
 604:                          if (Math.Abs(E[FINISH + o_e]) > TINY)
 605:                          {
 606:                              FINISH = FINISH + 1;
 607:                              goto LABEL20;
 608:                          }
 609:                      }
 610:                      // *
 611:                      // *           (Sub) Problem determined.  Compute its size and solve it.
 612:                      // *
 613:                      M = FINISH - START + 1;
 614:                      if (M == 1)
 615:                      {
 616:                          START = FINISH + 1;
 617:                          goto LABEL10;
 618:                      }
 619:                      if (M > SMLSIZ)
 620:                      {
 621:                          // *
 622:                          // *              Scale.
 623:                          // *
 624:                          ORGNRM = this._dlanst.Run("M", M, D, START + o_d, E, START + o_e);
 625:                          this._dlascl.Run("G", 0, 0, ORGNRM, ONE, M
 626:                                           , 1, ref D, START + o_d, M, ref INFO);
 627:                          this._dlascl.Run("G", 0, 0, ORGNRM, ONE, M - 1
 628:                                           , 1, ref E, START + o_e, M - 1, ref INFO);
 629:                          // *
 630:                          if (ICOMPZ == 1)
 631:                          {
 632:                              STRTRW = 1;
 633:                          }
 634:                          else
 635:                          {
 636:                              STRTRW = START;
 637:                          }
 638:                          this._dlaed0.Run(ICOMPZ, N, M, ref D, START + o_d, ref E, START + o_e, ref Z, STRTRW+START * LDZ + o_z
 639:                                           , LDZ, ref WORK, 1 + o_work, N, ref WORK, STOREZ + o_work, ref IWORK, offset_iwork, ref INFO);
 640:                          if (INFO != 0)
 641:                          {
 642:                              INFO = (INFO / (M + 1) + START - 1) * (N + 1) + FortranLib.Mod(INFO,(M + 1)) + START - 1;
 643:                              goto LABEL50;
 644:                          }
 645:                          // *
 646:                          // *              Scale back.
 647:                          // *
 648:                          this._dlascl.Run("G", 0, 0, ONE, ORGNRM, M
 649:                                           , 1, ref D, START + o_d, M, ref INFO);
 650:                          // *
 651:                      }
 652:                      else
 653:                      {
 654:                          if (ICOMPZ == 1)
 655:                          {
 656:                              // *
 657:                              // *                 Since QR won't update a Z matrix which is larger than
 658:                              // *                 the length of D, we must solve the sub-problem in a
 659:                              // *                 workspace and then multiply back into Z.
 660:                              // *
 661:                              this._dsteqr.Run("I", M, ref D, START + o_d, ref E, START + o_e, ref WORK, offset_work, M
 662:                                               , ref WORK, M * M + 1 + o_work, ref INFO);
 663:                              this._dlacpy.Run("A", N, M, Z, 1+START * LDZ + o_z, LDZ, ref WORK, STOREZ + o_work
 664:                                               , N);
 665:                              this._dgemm.Run("N", "N", N, M, M, ONE
 666:                                              , WORK, STOREZ + o_work, N, WORK, offset_work, M, ZERO, ref Z, 1+START * LDZ + o_z
 667:                                              , LDZ);
 668:                          }
 669:                          else
 670:                          {
 671:                              if (ICOMPZ == 2)
 672:                              {
 673:                                  this._dsteqr.Run("I", M, ref D, START + o_d, ref E, START + o_e, ref Z, START+START * LDZ + o_z, LDZ
 674:                                                   , ref WORK, offset_work, ref INFO);
 675:                              }
 676:                              else
 677:                              {
 678:                                  this._dsterf.Run(M, ref D, START + o_d, ref E, START + o_e, ref INFO);
 679:                              }
 680:                          }
 681:                          if (INFO != 0)
 682:                          {
 683:                              INFO = START * (N + 1) + FINISH;
 684:                              goto LABEL50;
 685:                          }
 686:                      }
 687:                      // *
 688:                      START = FINISH + 1;
 689:                      goto LABEL10;
 690:                  }
 691:                  // *
 692:                  // *        endwhile
 693:                  // *
 694:                  // *        If the problem split any number of times, then the eigenvalues
 695:                  // *        will not be properly ordered.  Here we permute the eigenvalues
 696:                  // *        (and the associated eigenvectors) into ascending order.
 697:                  // *
 698:                  if (M != N)
 699:                  {
 700:                      if (ICOMPZ == 0)
 701:                      {
 702:                          // *
 703:                          // *              Use Quick Sort
 704:                          // *
 705:                          this._dlasrt.Run("I", N, ref D, offset_d, ref INFO);
 706:                          // *
 707:                      }
 708:                      else
 709:                      {
 710:                          // *
 711:                          // *              Use Selection Sort to minimize swaps of eigenvectors
 712:                          // *
 713:                          for (II = 2; II <= N; II++)
 714:                          {
 715:                              I = II - 1;
 716:                              K = I;
 717:                              P = D[I + o_d];
 718:                              for (J = II; J <= N; J++)
 719:                              {
 720:                                  if (D[J + o_d] < P)
 721:                                  {
 722:                                      K = J;
 723:                                      P = D[J + o_d];
 724:                                  }
 725:                              }
 726:                              if (K != I)
 727:                              {
 728:                                  D[K + o_d] = D[I + o_d];
 729:                                  D[I + o_d] = P;
 730:                                  this._dswap.Run(N, ref Z, 1+I * LDZ + o_z, 1, ref Z, 1+K * LDZ + o_z, 1);
 731:                              }
 732:                          }
 733:                      }
 734:                  }
 735:              }
 736:              // *
 737:          LABEL50:;
 738:              WORK[1 + o_work] = LWMIN;
 739:              IWORK[1 + o_iwork] = LIWMIN;
 740:              // *
 741:              return;
 742:              // *
 743:              // *     End of DSTEDC
 744:              // *
 745:   
 746:              #endregion
 747:   
 748:          }
 749:      }
 750:  }