`   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:  }`