`   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:      /// DTREVC computes some or all of the right and/or left eigenvectors of`
`  27:      /// a real upper quasi-triangular matrix T.`
`  28:      /// Matrices of this type are produced by the Schur factorization of`
`  29:      /// a real general matrix:  A = Q*T*Q**T, as computed by DHSEQR.`
`  30:      /// `
`  31:      /// The right eigenvector x and the left eigenvector y of T corresponding`
`  32:      /// to an eigenvalue w are defined by:`
`  33:      /// `
`  34:      /// T*x = w*x,     (y**H)*T = w*(y**H)`
`  35:      /// `
`  36:      /// where y**H denotes the conjugate transpose of y.`
`  37:      /// The eigenvalues are not input to this routine, but are read directly`
`  38:      /// from the diagonal blocks of T.`
`  39:      /// `
`  40:      /// This routine returns the matrices X and/or Y of right and left`
`  41:      /// eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an`
`  42:      /// input matrix.  If Q is the orthogonal factor that reduces a matrix`
`  43:      /// A to Schur form T, then Q*X and Q*Y are the matrices of right and`
`  44:      /// left eigenvectors of A.`
`  45:      /// `
`  46:      ///</summary>`
`  47:      public class DTREVC`
`  48:      {`
`  49:      `
`  50:   `
`  51:          #region Dependencies`
`  52:          `
`  53:          LSAME _lsame; IDAMAX _idamax; DDOT _ddot; DLAMCH _dlamch; DAXPY _daxpy; DCOPY _dcopy; DGEMV _dgemv; DLALN2 _dlaln2; `
`  54:          DSCAL _dscal;XERBLA _xerbla; DLABAD _dlabad; `
`  55:   `
`  56:          #endregion`
`  57:   `
`  58:   `
`  59:          #region Fields`
`  60:          `
`  61:          const double ZERO = 0.0E+0; const double ONE = 1.0E+0; bool ALLV = false; bool BOTHV = false; bool LEFTV = false; `
`  62:          bool OVER = false;bool PAIR = false; bool RIGHTV = false; bool SOMEV = false; int I = 0; int IERR = 0; int II = 0; `
`  63:          int IP = 0;int IS = 0; int J = 0; int J1 = 0; int J2 = 0; int JNXT = 0; int K = 0; int KI = 0; int N2 = 0; `
`  64:          double BETA = 0;double BIGNUM = 0; double EMAX = 0; double OVFL = 0; double REC = 0; double REMAX = 0; double SCALE = 0; `
`  65:          double SMIN = 0;double SMLNUM = 0; double ULP = 0; double UNFL = 0; double VCRIT = 0; double VMAX = 0; double WI = 0; `
`  66:          double WR = 0;double XNORM = 0; double[] X = new double[2 * 2]; int offset_x = 0; int o_x = -3; `
`  67:   `
`  68:          #endregion`
`  69:   `
`  70:          public DTREVC(LSAME lsame, IDAMAX idamax, DDOT ddot, DLAMCH dlamch, DAXPY daxpy, DCOPY dcopy, DGEMV dgemv, DLALN2 dlaln2, DSCAL dscal, XERBLA xerbla`
`  71:                        , DLABAD dlabad)`
`  72:          {`
`  73:      `
`  74:   `
`  75:              #region Set Dependencies`
`  76:              `
`  77:              this._lsame = lsame; this._idamax = idamax; this._ddot = ddot; this._dlamch = dlamch; this._daxpy = daxpy; `
`  78:              this._dcopy = dcopy;this._dgemv = dgemv; this._dlaln2 = dlaln2; this._dscal = dscal; this._xerbla = xerbla; `
`  79:              this._dlabad = dlabad;`
`  80:   `
`  81:              #endregion`
`  82:   `
`  83:          }`
`  84:      `
`  85:          public DTREVC()`
`  86:          {`
`  87:      `
`  88:   `
`  89:              #region Dependencies (Initialization)`
`  90:              `
`  91:              LSAME lsame = new LSAME();`
`  92:              IDAMAX idamax = new IDAMAX();`
`  93:              DDOT ddot = new DDOT();`
`  94:              DLAMC3 dlamc3 = new DLAMC3();`
`  95:              DAXPY daxpy = new DAXPY();`
`  96:              DCOPY dcopy = new DCOPY();`
`  97:              XERBLA xerbla = new XERBLA();`
`  98:              DLADIV dladiv = new DLADIV();`
`  99:              DSCAL dscal = new DSCAL();`
` 100:              DLABAD dlabad = new DLABAD();`
` 101:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);`
` 102:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);`
` 103:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);`
` 104:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);`
` 105:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);`
` 106:              DGEMV dgemv = new DGEMV(lsame, xerbla);`
` 107:              DLALN2 dlaln2 = new DLALN2(dlamch, dladiv);`
` 108:   `
` 109:              #endregion`
` 110:   `
` 111:   `
` 112:              #region Set Dependencies`
` 113:              `
` 114:              this._lsame = lsame; this._idamax = idamax; this._ddot = ddot; this._dlamch = dlamch; this._daxpy = daxpy; `
` 115:              this._dcopy = dcopy;this._dgemv = dgemv; this._dlaln2 = dlaln2; this._dscal = dscal; this._xerbla = xerbla; `
` 116:              this._dlabad = dlabad;`
` 117:   `
` 118:              #endregion`
` 119:   `
` 120:          }`
` 121:          /// <summary>`
` 122:          /// Purpose`
` 123:          /// =======`
` 124:          /// `
` 125:          /// DTREVC computes some or all of the right and/or left eigenvectors of`
` 126:          /// a real upper quasi-triangular matrix T.`
` 127:          /// Matrices of this type are produced by the Schur factorization of`
` 128:          /// a real general matrix:  A = Q*T*Q**T, as computed by DHSEQR.`
` 129:          /// `
` 130:          /// The right eigenvector x and the left eigenvector y of T corresponding`
` 131:          /// to an eigenvalue w are defined by:`
` 132:          /// `
` 133:          /// T*x = w*x,     (y**H)*T = w*(y**H)`
` 134:          /// `
` 135:          /// where y**H denotes the conjugate transpose of y.`
` 136:          /// The eigenvalues are not input to this routine, but are read directly`
` 137:          /// from the diagonal blocks of T.`
` 138:          /// `
` 139:          /// This routine returns the matrices X and/or Y of right and left`
` 140:          /// eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an`
` 141:          /// input matrix.  If Q is the orthogonal factor that reduces a matrix`
` 142:          /// A to Schur form T, then Q*X and Q*Y are the matrices of right and`
` 143:          /// left eigenvectors of A.`
` 144:          /// `
` 145:          ///</summary>`
` 146:          /// <param name="SIDE">`
` 147:          /// (input) CHARACTER*1`
` 148:          /// = 'R':  compute right eigenvectors only;`
` 149:          /// = 'L':  compute left eigenvectors only;`
` 150:          /// = 'B':  compute both right and left eigenvectors.`
` 151:          ///</param>`
` 152:          /// <param name="HOWMNY">`
` 153:          /// (input) CHARACTER*1`
` 154:          /// = 'A':  compute all right and/or left eigenvectors;`
` 155:          /// = 'B':  compute all right and/or left eigenvectors,`
` 156:          /// backtransformed by the matrices in VR and/or VL;`
` 157:          /// = 'S':  compute selected right and/or left eigenvectors,`
` 158:          /// as indicated by the logical array SELECT.`
` 159:          ///</param>`
` 160:          /// <param name="SELECT">`
` 161:          /// (input/output) LOGICAL array, dimension (N)`
` 162:          /// If HOWMNY = 'S', SELECT specifies the eigenvectors to be`
` 163:          /// computed.`
` 164:          /// If w(j) is a real eigenvalue, the corresponding real`
` 165:          /// eigenvector is computed if SELECT(j) is .TRUE..`
` 166:          /// If w(j) and w(j+1) are the real and imaginary parts of a`
` 167:          /// complex eigenvalue, the corresponding complex eigenvector is`
` 168:          /// computed if either SELECT(j) or SELECT(j+1) is .TRUE., and`
` 169:          /// on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to`
` 170:          /// .FALSE..`
` 171:          /// Not referenced if HOWMNY = 'A' or 'B'.`
` 172:          ///</param>`
` 173:          /// <param name="N">`
` 174:          /// (input) INTEGER`
` 175:          /// The order of the matrix T. N .GE. 0.`
` 176:          ///</param>`
` 177:          /// <param name="T">`
` 178:          /// (input) DOUBLE PRECISION array, dimension (LDT,N)`
` 179:          /// The upper quasi-triangular matrix T in Schur canonical form.`
` 180:          ///</param>`
` 181:          /// <param name="LDT">`
` 182:          /// (input) INTEGER`
` 183:          /// The leading dimension of the array T. LDT .GE. max(1,N).`
` 184:          ///</param>`
` 185:          /// <param name="VL">`
` 186:          /// (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)`
` 187:          /// On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must`
` 188:          /// contain an N-by-N matrix Q (usually the orthogonal matrix Q`
` 189:          /// of Schur vectors returned by DHSEQR).`
` 190:          /// On exit, if SIDE = 'L' or 'B', VL contains:`
` 191:          /// if HOWMNY = 'A', the matrix Y of left eigenvectors of T;`
` 192:          /// if HOWMNY = 'B', the matrix Q*Y;`
` 193:          /// if HOWMNY = 'S', the left eigenvectors of T specified by`
` 194:          /// SELECT, stored consecutively in the columns`
` 195:          /// of VL, in the same order as their`
` 196:          /// eigenvalues.`
` 197:          /// A complex eigenvector corresponding to a complex eigenvalue`
` 198:          /// is stored in two consecutive columns, the first holding the`
` 199:          /// real part, and the second the imaginary part.`
` 200:          /// Not referenced if SIDE = 'R'.`
` 201:          ///</param>`
` 202:          /// <param name="LDVL">`
` 203:          /// (input) INTEGER`
` 204:          /// The leading dimension of the array VL.  LDVL .GE. 1, and if`
` 205:          /// SIDE = 'L' or 'B', LDVL .GE. N.`
` 206:          ///</param>`
` 207:          /// <param name="VR">`
` 208:          /// (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)`
` 209:          /// On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must`
` 210:          /// contain an N-by-N matrix Q (usually the orthogonal matrix Q`
` 211:          /// of Schur vectors returned by DHSEQR).`
` 212:          /// On exit, if SIDE = 'R' or 'B', VR contains:`
` 213:          /// if HOWMNY = 'A', the matrix X of right eigenvectors of T;`
` 214:          /// if HOWMNY = 'B', the matrix Q*X;`
` 215:          /// if HOWMNY = 'S', the right eigenvectors of T specified by`
` 216:          /// SELECT, stored consecutively in the columns`
` 217:          /// of VR, in the same order as their`
` 218:          /// eigenvalues.`
` 219:          /// A complex eigenvector corresponding to a complex eigenvalue`
` 220:          /// is stored in two consecutive columns, the first holding the`
` 221:          /// real part and the second the imaginary part.`
` 222:          /// Not referenced if SIDE = 'L'.`
` 223:          ///</param>`
` 224:          /// <param name="LDVR">`
` 225:          /// (input) INTEGER`
` 226:          /// The leading dimension of the array VR.  LDVR .GE. 1, and if`
` 227:          /// SIDE = 'R' or 'B', LDVR .GE. N.`
` 228:          ///</param>`
` 229:          /// <param name="MM">`
` 230:          /// (input) INTEGER`
` 231:          /// The number of columns in the arrays VL and/or VR. MM .GE. M.`
` 232:          ///</param>`
` 233:          /// <param name="M">`
` 234:          /// (output) INTEGER`
` 235:          /// The number of columns in the arrays VL and/or VR actually`
` 236:          /// used to store the eigenvectors.`
` 237:          /// If HOWMNY = 'A' or 'B', M is set to N.`
` 238:          /// Each selected real eigenvector occupies one column and each`
` 239:          /// selected complex eigenvector occupies two columns.`
` 240:          ///</param>`
` 241:          /// <param name="WORK">`
` 242:          /// (workspace) DOUBLE PRECISION array, dimension (3*N)`
` 243:          ///</param>`
` 244:          /// <param name="INFO">`
` 245:          /// (output) INTEGER`
` 246:          /// = 0:  successful exit`
` 247:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value`
` 248:          ///</param>`
` 249:          public void Run(string SIDE, string HOWMNY, ref bool[] SELECT, int offset_select, int N, double[] T, int offset_t, int LDT`
` 250:                           , ref double[] VL, int offset_vl, int LDVL, ref double[] VR, int offset_vr, int LDVR, int MM, ref int M`
` 251:                           , ref double[] WORK, int offset_work, ref int INFO)`
` 252:          {`
` 253:   `
` 254:              #region Array Index Correction`
` 255:              `
` 256:               int o_select = -1 + offset_select;  int o_t = -1 - LDT + offset_t;  int o_vl = -1 - LDVL + offset_vl; `
` 257:               int o_vr = -1 - LDVR + offset_vr; int o_work = -1 + offset_work; `
` 258:   `
` 259:              #endregion`
` 260:   `
` 261:   `
` 262:              #region Strings`
` 263:              `
` 264:              SIDE = SIDE.Substring(0, 1);  HOWMNY = HOWMNY.Substring(0, 1);  `
` 265:   `
` 266:              #endregion`
` 267:   `
` 268:   `
` 269:              #region Prolog`
` 270:              `
` 271:              // *`
` 272:              // *  -- LAPACK routine (version 3.1) --`
` 273:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..`
` 274:              // *     November 2006`
` 275:              // *`
` 276:              // *     .. Scalar Arguments ..`
` 277:              // *     ..`
` 278:              // *     .. Array Arguments ..`
` 279:              // *     ..`
` 280:              // *`
` 281:              // *  Purpose`
` 282:              // *  =======`
` 283:              // *`
` 284:              // *  DTREVC computes some or all of the right and/or left eigenvectors of`
` 285:              // *  a real upper quasi-triangular matrix T.`
` 286:              // *  Matrices of this type are produced by the Schur factorization of`
` 287:              // *  a real general matrix:  A = Q*T*Q**T, as computed by DHSEQR.`
` 288:              // *  `
` 289:              // *  The right eigenvector x and the left eigenvector y of T corresponding`
` 290:              // *  to an eigenvalue w are defined by:`
` 291:              // *  `
` 292:              // *     T*x = w*x,     (y**H)*T = w*(y**H)`
` 293:              // *  `
` 294:              // *  where y**H denotes the conjugate transpose of y.`
` 295:              // *  The eigenvalues are not input to this routine, but are read directly`
` 296:              // *  from the diagonal blocks of T.`
` 297:              // *  `
` 298:              // *  This routine returns the matrices X and/or Y of right and left`
` 299:              // *  eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an`
` 300:              // *  input matrix.  If Q is the orthogonal factor that reduces a matrix`
` 301:              // *  A to Schur form T, then Q*X and Q*Y are the matrices of right and`
` 302:              // *  left eigenvectors of A.`
` 303:              // *`
` 304:              // *  Arguments`
` 305:              // *  =========`
` 306:              // *`
` 307:              // *  SIDE    (input) CHARACTER*1`
` 308:              // *          = 'R':  compute right eigenvectors only;`
` 309:              // *          = 'L':  compute left eigenvectors only;`
` 310:              // *          = 'B':  compute both right and left eigenvectors.`
` 311:              // *`
` 312:              // *  HOWMNY  (input) CHARACTER*1`
` 313:              // *          = 'A':  compute all right and/or left eigenvectors;`
` 314:              // *          = 'B':  compute all right and/or left eigenvectors,`
` 315:              // *                  backtransformed by the matrices in VR and/or VL;`
` 316:              // *          = 'S':  compute selected right and/or left eigenvectors,`
` 317:              // *                  as indicated by the logical array SELECT.`
` 318:              // *`
` 319:              // *  SELECT  (input/output) LOGICAL array, dimension (N)`
` 320:              // *          If HOWMNY = 'S', SELECT specifies the eigenvectors to be`
` 321:              // *          computed.`
` 322:              // *          If w(j) is a real eigenvalue, the corresponding real`
` 323:              // *          eigenvector is computed if SELECT(j) is .TRUE..`
` 324:              // *          If w(j) and w(j+1) are the real and imaginary parts of a`
` 325:              // *          complex eigenvalue, the corresponding complex eigenvector is`
` 326:              // *          computed if either SELECT(j) or SELECT(j+1) is .TRUE., and`
` 327:              // *          on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to`
` 328:              // *          .FALSE..`
` 329:              // *          Not referenced if HOWMNY = 'A' or 'B'.`
` 330:              // *`
` 331:              // *  N       (input) INTEGER`
` 332:              // *          The order of the matrix T. N >= 0.`
` 333:              // *`
` 334:              // *  T       (input) DOUBLE PRECISION array, dimension (LDT,N)`
` 335:              // *          The upper quasi-triangular matrix T in Schur canonical form.`
` 336:              // *`
` 337:              // *  LDT     (input) INTEGER`
` 338:              // *          The leading dimension of the array T. LDT >= max(1,N).`
` 339:              // *`
` 340:              // *  VL      (input/output) DOUBLE PRECISION array, dimension (LDVL,MM)`
` 341:              // *          On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must`
` 342:              // *          contain an N-by-N matrix Q (usually the orthogonal matrix Q`
` 343:              // *          of Schur vectors returned by DHSEQR).`
` 344:              // *          On exit, if SIDE = 'L' or 'B', VL contains:`
` 345:              // *          if HOWMNY = 'A', the matrix Y of left eigenvectors of T;`
` 346:              // *          if HOWMNY = 'B', the matrix Q*Y;`
` 347:              // *          if HOWMNY = 'S', the left eigenvectors of T specified by`
` 348:              // *                           SELECT, stored consecutively in the columns`
` 349:              // *                           of VL, in the same order as their`
` 350:              // *                           eigenvalues.`
` 351:              // *          A complex eigenvector corresponding to a complex eigenvalue`
` 352:              // *          is stored in two consecutive columns, the first holding the`
` 353:              // *          real part, and the second the imaginary part.`
` 354:              // *          Not referenced if SIDE = 'R'.`
` 355:              // *`
` 356:              // *  LDVL    (input) INTEGER`
` 357:              // *          The leading dimension of the array VL.  LDVL >= 1, and if`
` 358:              // *          SIDE = 'L' or 'B', LDVL >= N.`
` 359:              // *`
` 360:              // *  VR      (input/output) DOUBLE PRECISION array, dimension (LDVR,MM)`
` 361:              // *          On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must`
` 362:              // *          contain an N-by-N matrix Q (usually the orthogonal matrix Q`
` 363:              // *          of Schur vectors returned by DHSEQR).`
` 364:              // *          On exit, if SIDE = 'R' or 'B', VR contains:`
` 365:              // *          if HOWMNY = 'A', the matrix X of right eigenvectors of T;`
` 366:              // *          if HOWMNY = 'B', the matrix Q*X;`
` 367:              // *          if HOWMNY = 'S', the right eigenvectors of T specified by`
` 368:              // *                           SELECT, stored consecutively in the columns`
` 369:              // *                           of VR, in the same order as their`
` 370:              // *                           eigenvalues.`
` 371:              // *          A complex eigenvector corresponding to a complex eigenvalue`
` 372:              // *          is stored in two consecutive columns, the first holding the`
` 373:              // *          real part and the second the imaginary part.`
` 374:              // *          Not referenced if SIDE = 'L'.`
` 375:              // *`
` 376:              // *  LDVR    (input) INTEGER`
` 377:              // *          The leading dimension of the array VR.  LDVR >= 1, and if`
` 378:              // *          SIDE = 'R' or 'B', LDVR >= N.`
` 379:              // *`
` 380:              // *  MM      (input) INTEGER`
` 381:              // *          The number of columns in the arrays VL and/or VR. MM >= M.`
` 382:              // *`
` 383:              // *  M       (output) INTEGER`
` 384:              // *          The number of columns in the arrays VL and/or VR actually`
` 385:              // *          used to store the eigenvectors.`
` 386:              // *          If HOWMNY = 'A' or 'B', M is set to N.`
` 387:              // *          Each selected real eigenvector occupies one column and each`
` 388:              // *          selected complex eigenvector occupies two columns.`
` 389:              // *`
` 390:              // *  WORK    (workspace) DOUBLE PRECISION array, dimension (3*N)`
` 391:              // *`
` 392:              // *  INFO    (output) INTEGER`
` 393:              // *          = 0:  successful exit`
` 394:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value`
` 395:              // *`
` 396:              // *  Further Details`
` 397:              // *  ===============`
` 398:              // *`
` 399:              // *  The algorithm used in this program is basically backward (forward)`
` 400:              // *  substitution, with scaling to make the the code robust against`
` 401:              // *  possible overflow.`
` 402:              // *`
` 403:              // *  Each eigenvector is normalized so that the element of largest`
` 404:              // *  magnitude has magnitude 1; here the magnitude of a complex number`
` 405:              // *  (x,y) is taken to be |x| + |y|.`
` 406:              // *`
` 407:              // *  =====================================================================`
` 408:              // *`
` 409:              // *     .. Parameters ..`
` 410:              // *     ..`
` 411:              // *     .. Local Scalars ..`
` 412:              // *     ..`
` 413:              // *     .. External Functions ..`
` 414:              // *     ..`
` 415:              // *     .. External Subroutines ..`
` 416:              // *     ..`
` 417:              // *     .. Intrinsic Functions ..`
` 418:              //      INTRINSIC          ABS, MAX, SQRT;`
` 419:              // *     ..`
` 420:              // *     .. Local Arrays ..`
` 421:              // *     ..`
` 422:              // *     .. Executable Statements ..`
` 423:              // *`
` 424:              // *     Decode and test the input parameters`
` 425:              // *`
` 426:   `
` 427:              #endregion`
` 428:   `
` 429:   `
` 430:              #region Body`
` 431:              `
` 432:              BOTHV = this._lsame.Run(SIDE, "B");`
` 433:              RIGHTV = this._lsame.Run(SIDE, "R") || BOTHV;`
` 434:              LEFTV = this._lsame.Run(SIDE, "L") || BOTHV;`
` 435:              // *`
` 436:              ALLV = this._lsame.Run(HOWMNY, "A");`
` 437:              OVER = this._lsame.Run(HOWMNY, "B");`
` 438:              SOMEV = this._lsame.Run(HOWMNY, "S");`
` 439:              // *`
` 440:              INFO = 0;`
` 441:              if (!RIGHTV && !LEFTV)`
` 442:              {`
` 443:                  INFO =  - 1;`
` 444:              }`
` 445:              else`
` 446:              {`
` 447:                  if (!ALLV && !OVER && !SOMEV)`
` 448:                  {`
` 449:                      INFO =  - 2;`
` 450:                  }`
` 451:                  else`
` 452:                  {`
` 453:                      if (N < 0)`
` 454:                      {`
` 455:                          INFO =  - 4;`
` 456:                      }`
` 457:                      else`
` 458:                      {`
` 459:                          if (LDT < Math.Max(1, N))`
` 460:                          {`
` 461:                              INFO =  - 6;`
` 462:                          }`
` 463:                          else`
` 464:                          {`
` 465:                              if (LDVL < 1 || (LEFTV && LDVL < N))`
` 466:                              {`
` 467:                                  INFO =  - 8;`
` 468:                              }`
` 469:                              else`
` 470:                              {`
` 471:                                  if (LDVR < 1 || (RIGHTV && LDVR < N))`
` 472:                                  {`
` 473:                                      INFO =  - 10;`
` 474:                                  }`
` 475:                                  else`
` 476:                                  {`
` 477:                                      // *`
` 478:                                      // *        Set M to the number of columns required to store the selected`
` 479:                                      // *        eigenvectors, standardize the array SELECT if necessary, and`
` 480:                                      // *        test MM.`
` 481:                                      // *`
` 482:                                      if (SOMEV)`
` 483:                                      {`
` 484:                                          M = 0;`
` 485:                                          PAIR = false;`
` 486:                                          for (J = 1; J <= N; J++)`
` 487:                                          {`
` 488:                                              if (PAIR)`
` 489:                                              {`
` 490:                                                  PAIR = false;`
` 491:                                                  SELECT[J + o_select] = false;`
` 492:                                              }`
` 493:                                              else`
` 494:                                              {`
` 495:                                                  if (J < N)`
` 496:                                                  {`
` 497:                                                      if (T[J + 1+J * LDT + o_t] == ZERO)`
` 498:                                                      {`
` 499:                                                          if (SELECT[J + o_select]) M = M + 1;`
` 500:                                                      }`
` 501:                                                      else`
` 502:                                                      {`
` 503:                                                          PAIR = true;`
` 504:                                                          if (SELECT[J + o_select] || SELECT[J + 1 + o_select])`
` 505:                                                          {`
` 506:                                                              SELECT[J + o_select] = true;`
` 507:                                                              M = M + 2;`
` 508:                                                          }`
` 509:                                                      }`
` 510:                                                  }`
` 511:                                                  else`
` 512:                                                  {`
` 513:                                                      if (SELECT[N + o_select]) M = M + 1;`
` 514:                                                  }`
` 515:                                              }`
` 516:                                          }`
` 517:                                      }`
` 518:                                      else`
` 519:                                      {`
` 520:                                          M = N;`
` 521:                                      }`
` 522:                                      // *`
` 523:                                      if (MM < M)`
` 524:                                      {`
` 525:                                          INFO =  - 11;`
` 526:                                      }`
` 527:                                  }`
` 528:                              }`
` 529:                          }`
` 530:                      }`
` 531:                  }`
` 532:              }`
` 533:              if (INFO != 0)`
` 534:              {`
` 535:                  this._xerbla.Run("DTREVC",  - INFO);`
` 536:                  return;`
` 537:              }`
` 538:              // *`
` 539:              // *     Quick return if possible.`
` 540:              // *`
` 541:              if (N == 0) return;`
` 542:              // *`
` 543:              // *     Set the constants to control overflow.`
` 544:              // *`
` 545:              UNFL = this._dlamch.Run("Safe minimum");`
` 546:              OVFL = ONE / UNFL;`
` 547:              this._dlabad.Run(ref UNFL, ref OVFL);`
` 548:              ULP = this._dlamch.Run("Precision");`
` 549:              SMLNUM = UNFL * (N / ULP);`
` 550:              BIGNUM = (ONE - ULP) / SMLNUM;`
` 551:              // *`
` 552:              // *     Compute 1-norm of each column of strictly upper triangular`
` 553:              // *     part of T to control overflow in triangular solver.`
` 554:              // *`
` 555:              WORK[1 + o_work] = ZERO;`
` 556:              for (J = 2; J <= N; J++)`
` 557:              {`
` 558:                  WORK[J + o_work] = ZERO;`
` 559:                  for (I = 1; I <= J - 1; I++)`
` 560:                  {`
` 561:                      WORK[J + o_work] = WORK[J + o_work] + Math.Abs(T[I+J * LDT + o_t]);`
` 562:                  }`
` 563:              }`
` 564:              // *`
` 565:              // *     Index IP is used to specify the real or complex eigenvalue:`
` 566:              // *       IP = 0, real eigenvalue,`
` 567:              // *            1, first of conjugate complex pair: (wr,wi)`
` 568:              // *           -1, second of conjugate complex pair: (wr,wi)`
` 569:              // *`
` 570:              N2 = 2 * N;`
` 571:              // *`
` 572:              if (RIGHTV)`
` 573:              {`
` 574:                  // *`
` 575:                  // *        Compute right eigenvectors.`
` 576:                  // *`
` 577:                  IP = 0;`
` 578:                  IS = M;`
` 579:                  for (KI = N; KI >= 1; KI +=  - 1)`
` 580:                  {`
` 581:                      // *`
` 582:                      if (IP == 1) goto LABEL130;`
` 583:                      if (KI == 1) goto LABEL40;`
` 584:                      if (T[KI+(KI - 1) * LDT + o_t] == ZERO) goto LABEL40;`
` 585:                      IP =  - 1;`
` 586:                      // *`
` 587:                  LABEL40:;`
` 588:                      if (SOMEV)`
` 589:                      {`
` 590:                          if (IP == 0)`
` 591:                          {`
` 592:                              if (!SELECT[KI + o_select]) goto LABEL130;`
` 593:                          }`
` 594:                          else`
` 595:                          {`
` 596:                              if (!SELECT[KI - 1 + o_select]) goto LABEL130;`
` 597:                          }`
` 598:                      }`
` 599:                      // *`
` 600:                      // *           Compute the KI-th eigenvalue (WR,WI).`
` 601:                      // *`
` 602:                      WR = T[KI+KI * LDT + o_t];`
` 603:                      WI = ZERO;`
` 604:                      if (IP != 0) WI = Math.Sqrt(Math.Abs(T[KI+(KI - 1) * LDT + o_t])) * Math.Sqrt(Math.Abs(T[KI - 1+KI * LDT + o_t]));`
` 605:                      SMIN = Math.Max(ULP * (Math.Abs(WR) + Math.Abs(WI)), SMLNUM);`
` 606:                      // *`
` 607:                      if (IP == 0)`
` 608:                      {`
` 609:                          // *`
` 610:                          // *              Real right eigenvector`
` 611:                          // *`
` 612:                          WORK[KI + N + o_work] = ONE;`
` 613:                          // *`
` 614:                          // *              Form right-hand side`
` 615:                          // *`
` 616:                          for (K = 1; K <= KI - 1; K++)`
` 617:                          {`
` 618:                              WORK[K + N + o_work] =  - T[K+KI * LDT + o_t];`
` 619:                          }`
` 620:                          // *`
` 621:                          // *              Solve the upper quasi-triangular system:`
` 622:                          // *                 (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.`
` 623:                          // *`
` 624:                          JNXT = KI - 1;`
` 625:                          for (J = KI - 1; J >= 1; J +=  - 1)`
` 626:                          {`
` 627:                              if (J > JNXT) goto LABEL60;`
` 628:                              J1 = J;`
` 629:                              J2 = J;`
` 630:                              JNXT = J - 1;`
` 631:                              if (J > 1)`
` 632:                              {`
` 633:                                  if (T[J+(J - 1) * LDT + o_t] != ZERO)`
` 634:                                  {`
` 635:                                      J1 = J - 1;`
` 636:                                      JNXT = J - 2;`
` 637:                                  }`
` 638:                              }`
` 639:                              // *`
` 640:                              if (J1 == J2)`
` 641:                              {`
` 642:                                  // *`
` 643:                                  // *                    1-by-1 diagonal block`
` 644:                                  // *`
` 645:                                  this._dlaln2.Run(false, 1, 1, SMIN, ONE, T, J+J * LDT + o_t`
` 646:                                                   , LDT, ONE, ONE, WORK, J + N + o_work, N, WR`
` 647:                                                   , ZERO, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);`
` 648:                                  // *`
` 649:                                  // *                    Scale X(1,1) to avoid overflow when updating`
` 650:                                  // *                    the right-hand side.`
` 651:                                  // *`
` 652:                                  if (XNORM > ONE)`
` 653:                                  {`
` 654:                                      if (WORK[J + o_work] > BIGNUM / XNORM)`
` 655:                                      {`
` 656:                                          X[1+1 * 2 + o_x] = X[1+1 * 2 + o_x] / XNORM;`
` 657:                                          SCALE = SCALE / XNORM;`
` 658:                                      }`
` 659:                                  }`
` 660:                                  // *`
` 661:                                  // *                    Scale if necessary`
` 662:                                  // *`
` 663:                                  if (SCALE != ONE) this._dscal.Run(KI, SCALE, ref WORK, 1 + N + o_work, 1);`
` 664:                                  WORK[J + N + o_work] = X[1+1 * 2 + o_x];`
` 665:                                  // *`
` 666:                                  // *                    Update right-hand side`
` 667:                                  // *`
` 668:                                  this._daxpy.Run(J - 1,  - X[1+1 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);`
` 669:                                  // *`
` 670:                              }`
` 671:                              else`
` 672:                              {`
` 673:                                  // *`
` 674:                                  // *                    2-by-2 diagonal block`
` 675:                                  // *`
` 676:                                  this._dlaln2.Run(false, 2, 1, SMIN, ONE, T, J - 1+(J - 1) * LDT + o_t`
` 677:                                                   , LDT, ONE, ONE, WORK, J - 1 + N + o_work, N, WR`
` 678:                                                   , ZERO, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);`
` 679:                                  // *`
` 680:                                  // *                    Scale X(1,1) and X(2,1) to avoid overflow when`
` 681:                                  // *                    updating the right-hand side.`
` 682:                                  // *`
` 683:                                  if (XNORM > ONE)`
` 684:                                  {`
` 685:                                      BETA = Math.Max(WORK[J - 1 + o_work], WORK[J + o_work]);`
` 686:                                      if (BETA > BIGNUM / XNORM)`
` 687:                                      {`
` 688:                                          X[1+1 * 2 + o_x] = X[1+1 * 2 + o_x] / XNORM;`
` 689:                                          X[2+1 * 2 + o_x] = X[2+1 * 2 + o_x] / XNORM;`
` 690:                                          SCALE = SCALE / XNORM;`
` 691:                                      }`
` 692:                                  }`
` 693:                                  // *`
` 694:                                  // *                    Scale if necessary`
` 695:                                  // *`
` 696:                                  if (SCALE != ONE) this._dscal.Run(KI, SCALE, ref WORK, 1 + N + o_work, 1);`
` 697:                                  WORK[J - 1 + N + o_work] = X[1+1 * 2 + o_x];`
` 698:                                  WORK[J + N + o_work] = X[2+1 * 2 + o_x];`
` 699:                                  // *`
` 700:                                  // *                    Update right-hand side`
` 701:                                  // *`
` 702:                                  this._daxpy.Run(J - 2,  - X[1+1 * 2 + o_x], T, 1+(J - 1) * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);`
` 703:                                  this._daxpy.Run(J - 2,  - X[2+1 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);`
` 704:                              }`
` 705:                          LABEL60:;`
` 706:                          }`
` 707:                          // *`
` 708:                          // *              Copy the vector x or Q*x to VR and normalize.`
` 709:                          // *`
` 710:                          if (!OVER)`
` 711:                          {`
` 712:                              this._dcopy.Run(KI, WORK, 1 + N + o_work, 1, ref VR, 1+IS * LDVR + o_vr, 1);`
` 713:                              // *`
` 714:                              II = this._idamax.Run(KI, VR, 1+IS * LDVR + o_vr, 1);`
` 715:                              REMAX = ONE / Math.Abs(VR[II+IS * LDVR + o_vr]);`
` 716:                              this._dscal.Run(KI, REMAX, ref VR, 1+IS * LDVR + o_vr, 1);`
` 717:                              // *`
` 718:                              for (K = KI + 1; K <= N; K++)`
` 719:                              {`
` 720:                                  VR[K+IS * LDVR + o_vr] = ZERO;`
` 721:                              }`
` 722:                          }`
` 723:                          else`
` 724:                          {`
` 725:                              if (KI > 1)`
` 726:                              {`
` 727:                                  this._dgemv.Run("N", N, KI - 1, ONE, VR, offset_vr, LDVR`
` 728:                                                  , WORK, 1 + N + o_work, 1, WORK[KI + N + o_work], ref VR, 1+KI * LDVR + o_vr, 1);`
` 729:                              }`
` 730:                              // *`
` 731:                              II = this._idamax.Run(N, VR, 1+KI * LDVR + o_vr, 1);`
` 732:                              REMAX = ONE / Math.Abs(VR[II+KI * LDVR + o_vr]);`
` 733:                              this._dscal.Run(N, REMAX, ref VR, 1+KI * LDVR + o_vr, 1);`
` 734:                          }`
` 735:                          // *`
` 736:                      }`
` 737:                      else`
` 738:                      {`
` 739:                          // *`
` 740:                          // *              Complex right eigenvector.`
` 741:                          // *`
` 742:                          // *              Initial solve`
` 743:                          // *                [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.`
` 744:                          // *                [ (T(KI,KI-1)   T(KI,KI)   )               ]`
` 745:                          // *`
` 746:                          if (Math.Abs(T[KI - 1+KI * LDT + o_t]) >= Math.Abs(T[KI+(KI - 1) * LDT + o_t]))`
` 747:                          {`
` 748:                              WORK[KI - 1 + N + o_work] = ONE;`
` 749:                              WORK[KI + N2 + o_work] = WI / T[KI - 1+KI * LDT + o_t];`
` 750:                          }`
` 751:                          else`
` 752:                          {`
` 753:                              WORK[KI - 1 + N + o_work] =  - WI / T[KI+(KI - 1) * LDT + o_t];`
` 754:                              WORK[KI + N2 + o_work] = ONE;`
` 755:                          }`
` 756:                          WORK[KI + N + o_work] = ZERO;`
` 757:                          WORK[KI - 1 + N2 + o_work] = ZERO;`
` 758:                          // *`
` 759:                          // *              Form right-hand side`
` 760:                          // *`
` 761:                          for (K = 1; K <= KI - 2; K++)`
` 762:                          {`
` 763:                              WORK[K + N + o_work] =  - WORK[KI - 1 + N + o_work] * T[K+(KI - 1) * LDT + o_t];`
` 764:                              WORK[K + N2 + o_work] =  - WORK[KI + N2 + o_work] * T[K+KI * LDT + o_t];`
` 765:                          }`
` 766:                          // *`
` 767:                          // *              Solve upper quasi-triangular system:`
` 768:                          // *              (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)`
` 769:                          // *`
` 770:                          JNXT = KI - 2;`
` 771:                          for (J = KI - 2; J >= 1; J +=  - 1)`
` 772:                          {`
` 773:                              if (J > JNXT) goto LABEL90;`
` 774:                              J1 = J;`
` 775:                              J2 = J;`
` 776:                              JNXT = J - 1;`
` 777:                              if (J > 1)`
` 778:                              {`
` 779:                                  if (T[J+(J - 1) * LDT + o_t] != ZERO)`
` 780:                                  {`
` 781:                                      J1 = J - 1;`
` 782:                                      JNXT = J - 2;`
` 783:                                  }`
` 784:                              }`
` 785:                              // *`
` 786:                              if (J1 == J2)`
` 787:                              {`
` 788:                                  // *`
` 789:                                  // *                    1-by-1 diagonal block`
` 790:                                  // *`
` 791:                                  this._dlaln2.Run(false, 1, 2, SMIN, ONE, T, J+J * LDT + o_t`
` 792:                                                   , LDT, ONE, ONE, WORK, J + N + o_work, N, WR`
` 793:                                                   , WI, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);`
` 794:                                  // *`
` 795:                                  // *                    Scale X(1,1) and X(1,2) to avoid overflow when`
` 796:                                  // *                    updating the right-hand side.`
` 797:                                  // *`
` 798:                                  if (XNORM > ONE)`
` 799:                                  {`
` 800:                                      if (WORK[J + o_work] > BIGNUM / XNORM)`
` 801:                                      {`
` 802:                                          X[1+1 * 2 + o_x] = X[1+1 * 2 + o_x] / XNORM;`
` 803:                                          X[1+2 * 2 + o_x] = X[1+2 * 2 + o_x] / XNORM;`
` 804:                                          SCALE = SCALE / XNORM;`
` 805:                                      }`
` 806:                                  }`
` 807:                                  // *`
` 808:                                  // *                    Scale if necessary`
` 809:                                  // *`
` 810:                                  if (SCALE != ONE)`
` 811:                                  {`
` 812:                                      this._dscal.Run(KI, SCALE, ref WORK, 1 + N + o_work, 1);`
` 813:                                      this._dscal.Run(KI, SCALE, ref WORK, 1 + N2 + o_work, 1);`
` 814:                                  }`
` 815:                                  WORK[J + N + o_work] = X[1+1 * 2 + o_x];`
` 816:                                  WORK[J + N2 + o_work] = X[1+2 * 2 + o_x];`
` 817:                                  // *`
` 818:                                  // *                    Update the right-hand side`
` 819:                                  // *`
` 820:                                  this._daxpy.Run(J - 1,  - X[1+1 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);`
` 821:                                  this._daxpy.Run(J - 1,  - X[1+2 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N2 + o_work, 1);`
` 822:                                  // *`
` 823:                              }`
` 824:                              else`
` 825:                              {`
` 826:                                  // *`
` 827:                                  // *                    2-by-2 diagonal block`
` 828:                                  // *`
` 829:                                  this._dlaln2.Run(false, 2, 2, SMIN, ONE, T, J - 1+(J - 1) * LDT + o_t`
` 830:                                                   , LDT, ONE, ONE, WORK, J - 1 + N + o_work, N, WR`
` 831:                                                   , WI, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);`
` 832:                                  // *`
` 833:                                  // *                    Scale X to avoid overflow when updating`
` 834:                                  // *                    the right-hand side.`
` 835:                                  // *`
` 836:                                  if (XNORM > ONE)`
` 837:                                  {`
` 838:                                      BETA = Math.Max(WORK[J - 1 + o_work], WORK[J + o_work]);`
` 839:                                      if (BETA > BIGNUM / XNORM)`
` 840:                                      {`
` 841:                                          REC = ONE / XNORM;`
` 842:                                          X[1+1 * 2 + o_x] = X[1+1 * 2 + o_x] * REC;`
` 843:                                          X[1+2 * 2 + o_x] = X[1+2 * 2 + o_x] * REC;`
` 844:                                          X[2+1 * 2 + o_x] = X[2+1 * 2 + o_x] * REC;`
` 845:                                          X[2+2 * 2 + o_x] = X[2+2 * 2 + o_x] * REC;`
` 846:                                          SCALE = SCALE * REC;`
` 847:                                      }`
` 848:                                  }`
` 849:                                  // *`
` 850:                                  // *                    Scale if necessary`
` 851:                                  // *`
` 852:                                  if (SCALE != ONE)`
` 853:                                  {`
` 854:                                      this._dscal.Run(KI, SCALE, ref WORK, 1 + N + o_work, 1);`
` 855:                                      this._dscal.Run(KI, SCALE, ref WORK, 1 + N2 + o_work, 1);`
` 856:                                  }`
` 857:                                  WORK[J - 1 + N + o_work] = X[1+1 * 2 + o_x];`
` 858:                                  WORK[J + N + o_work] = X[2+1 * 2 + o_x];`
` 859:                                  WORK[J - 1 + N2 + o_work] = X[1+2 * 2 + o_x];`
` 860:                                  WORK[J + N2 + o_work] = X[2+2 * 2 + o_x];`
` 861:                                  // *`
` 862:                                  // *                    Update the right-hand side`
` 863:                                  // *`
` 864:                                  this._daxpy.Run(J - 2,  - X[1+1 * 2 + o_x], T, 1+(J - 1) * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);`
` 865:                                  this._daxpy.Run(J - 2,  - X[2+1 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N + o_work, 1);`
` 866:                                  this._daxpy.Run(J - 2,  - X[1+2 * 2 + o_x], T, 1+(J - 1) * LDT + o_t, 1, ref WORK, 1 + N2 + o_work, 1);`
` 867:                                  this._daxpy.Run(J - 2,  - X[2+2 * 2 + o_x], T, 1+J * LDT + o_t, 1, ref WORK, 1 + N2 + o_work, 1);`
` 868:                              }`
` 869:                          LABEL90:;`
` 870:                          }`
` 871:                          // *`
` 872:                          // *              Copy the vector x or Q*x to VR and normalize.`
` 873:                          // *`
` 874:                          if (!OVER)`
` 875:                          {`
` 876:                              this._dcopy.Run(KI, WORK, 1 + N + o_work, 1, ref VR, 1+(IS - 1) * LDVR + o_vr, 1);`
` 877:                              this._dcopy.Run(KI, WORK, 1 + N2 + o_work, 1, ref VR, 1+IS * LDVR + o_vr, 1);`
` 878:                              // *`
` 879:                              EMAX = ZERO;`
` 880:                              for (K = 1; K <= KI; K++)`
` 881:                              {`
` 882:                                  EMAX = Math.Max(EMAX, Math.Abs(VR[K+(IS - 1) * LDVR + o_vr]) + Math.Abs(VR[K+IS * LDVR + o_vr]));`
` 883:                              }`
` 884:                              // *`
` 885:                              REMAX = ONE / EMAX;`
` 886:                              this._dscal.Run(KI, REMAX, ref VR, 1+(IS - 1) * LDVR + o_vr, 1);`
` 887:                              this._dscal.Run(KI, REMAX, ref VR, 1+IS * LDVR + o_vr, 1);`
` 888:                              // *`
` 889:                              for (K = KI + 1; K <= N; K++)`
` 890:                              {`
` 891:                                  VR[K+(IS - 1) * LDVR + o_vr] = ZERO;`
` 892:                                  VR[K+IS * LDVR + o_vr] = ZERO;`
` 893:                              }`
` 894:                              // *`
` 895:                          }`
` 896:                          else`
` 897:                          {`
` 898:                              // *`
` 899:                              if (KI > 2)`
` 900:                              {`
` 901:                                  this._dgemv.Run("N", N, KI - 2, ONE, VR, offset_vr, LDVR`
` 902:                                                  , WORK, 1 + N + o_work, 1, WORK[KI - 1 + N + o_work], ref VR, 1+(KI - 1) * LDVR + o_vr, 1);`
` 903:                                  this._dgemv.Run("N", N, KI - 2, ONE, VR, offset_vr, LDVR`
` 904:                                                  , WORK, 1 + N2 + o_work, 1, WORK[KI + N2 + o_work], ref VR, 1+KI * LDVR + o_vr, 1);`
` 905:                              }`
` 906:                              else`
` 907:                              {`
` 908:                                  this._dscal.Run(N, WORK[KI - 1 + N + o_work], ref VR, 1+(KI - 1) * LDVR + o_vr, 1);`
` 909:                                  this._dscal.Run(N, WORK[KI + N2 + o_work], ref VR, 1+KI * LDVR + o_vr, 1);`
` 910:                              }`
` 911:                              // *`
` 912:                              EMAX = ZERO;`
` 913:                              for (K = 1; K <= N; K++)`
` 914:                              {`
` 915:                                  EMAX = Math.Max(EMAX, Math.Abs(VR[K+(KI - 1) * LDVR + o_vr]) + Math.Abs(VR[K+KI * LDVR + o_vr]));`
` 916:                              }`
` 917:                              REMAX = ONE / EMAX;`
` 918:                              this._dscal.Run(N, REMAX, ref VR, 1+(KI - 1) * LDVR + o_vr, 1);`
` 919:                              this._dscal.Run(N, REMAX, ref VR, 1+KI * LDVR + o_vr, 1);`
` 920:                          }`
` 921:                      }`
` 922:                      // *`
` 923:                      IS = IS - 1;`
` 924:                      if (IP != 0) IS = IS - 1;`
` 925:                  LABEL130:;`
` 926:                      if (IP == 1) IP = 0;`
` 927:                      if (IP ==  - 1) IP = 1;`
` 928:                  }`
` 929:              }`
` 930:              // *`
` 931:              if (LEFTV)`
` 932:              {`
` 933:                  // *`
` 934:                  // *        Compute left eigenvectors.`
` 935:                  // *`
` 936:                  IP = 0;`
` 937:                  IS = 1;`
` 938:                  for (KI = 1; KI <= N; KI++)`
` 939:                  {`
` 940:                      // *`
` 941:                      if (IP ==  - 1) goto LABEL250;`
` 942:                      if (KI == N) goto LABEL150;`
` 943:                      if (T[KI + 1+KI * LDT + o_t] == ZERO) goto LABEL150;`
` 944:                      IP = 1;`
` 945:                      // *`
` 946:                  LABEL150:;`
` 947:                      if (SOMEV)`
` 948:                      {`
` 949:                          if (!SELECT[KI + o_select]) goto LABEL250;`
` 950:                      }`
` 951:                      // *`
` 952:                      // *           Compute the KI-th eigenvalue (WR,WI).`
` 953:                      // *`
` 954:                      WR = T[KI+KI * LDT + o_t];`
` 955:                      WI = ZERO;`
` 956:                      if (IP != 0) WI = Math.Sqrt(Math.Abs(T[KI+(KI + 1) * LDT + o_t])) * Math.Sqrt(Math.Abs(T[KI + 1+KI * LDT + o_t]));`
` 957:                      SMIN = Math.Max(ULP * (Math.Abs(WR) + Math.Abs(WI)), SMLNUM);`
` 958:                      // *`
` 959:                      if (IP == 0)`
` 960:                      {`
` 961:                          // *`
` 962:                          // *              Real left eigenvector.`
` 963:                          // *`
` 964:                          WORK[KI + N + o_work] = ONE;`
` 965:                          // *`
` 966:                          // *              Form right-hand side`
` 967:                          // *`
` 968:                          for (K = KI + 1; K <= N; K++)`
` 969:                          {`
` 970:                              WORK[K + N + o_work] =  - T[KI+K * LDT + o_t];`
` 971:                          }`
` 972:                          // *`
` 973:                          // *              Solve the quasi-triangular system:`
` 974:                          // *                 (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK`
` 975:                          // *`
` 976:                          VMAX = ONE;`
` 977:                          VCRIT = BIGNUM;`
` 978:                          // *`
` 979:                          JNXT = KI + 1;`
` 980:                          for (J = KI + 1; J <= N; J++)`
` 981:                          {`
` 982:                              if (J < JNXT) goto LABEL170;`
` 983:                              J1 = J;`
` 984:                              J2 = J;`
` 985:                              JNXT = J + 1;`
` 986:                              if (J < N)`
` 987:                              {`
` 988:                                  if (T[J + 1+J * LDT + o_t] != ZERO)`
` 989:                                  {`
` 990:                                      J2 = J + 1;`
` 991:                                      JNXT = J + 2;`
` 992:                                  }`
` 993:                              }`
` 994:                              // *`
` 995:                              if (J1 == J2)`
` 996:                              {`
` 997:                                  // *`
` 998:                                  // *                    1-by-1 diagonal block`
` 999:                                  // *`
`1000:                                  // *                    Scale if necessary to avoid overflow when forming`
`1001:                                  // *                    the right-hand side.`
`1002:                                  // *`
`1003:                                  if (WORK[J + o_work] > VCRIT)`
`1004:                                  {`
`1005:                                      REC = ONE / VMAX;`
`1006:                                      this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N + o_work, 1);`
`1007:                                      VMAX = ONE;`
`1008:                                      VCRIT = BIGNUM;`
`1009:                                  }`
`1010:                                  // *`
`1011:                                  WORK[J + N + o_work] = WORK[J + N + o_work] - this._ddot.Run(J - KI - 1, T, KI + 1+J * LDT + o_t, 1, WORK, KI + 1 + N + o_work, 1);`
`1012:                                  // *`
`1013:                                  // *                    Solve (T(J,J)-WR)'*X = WORK`
`1014:                                  // *`
`1015:                                  this._dlaln2.Run(false, 1, 1, SMIN, ONE, T, J+J * LDT + o_t`
`1016:                                                   , LDT, ONE, ONE, WORK, J + N + o_work, N, WR`
`1017:                                                   , ZERO, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);`
`1018:                                  // *`
`1019:                                  // *                    Scale if necessary`
`1020:                                  // *`
`1021:                                  if (SCALE != ONE) this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N + o_work, 1);`
`1022:                                  WORK[J + N + o_work] = X[1+1 * 2 + o_x];`
`1023:                                  VMAX = Math.Max(Math.Abs(WORK[J + N + o_work]), VMAX);`
`1024:                                  VCRIT = BIGNUM / VMAX;`
`1025:                                  // *`
`1026:                              }`
`1027:                              else`
`1028:                              {`
`1029:                                  // *`
`1030:                                  // *                    2-by-2 diagonal block`
`1031:                                  // *`
`1032:                                  // *                    Scale if necessary to avoid overflow when forming`
`1033:                                  // *                    the right-hand side.`
`1034:                                  // *`
`1035:                                  BETA = Math.Max(WORK[J + o_work], WORK[J + 1 + o_work]);`
`1036:                                  if (BETA > VCRIT)`
`1037:                                  {`
`1038:                                      REC = ONE / VMAX;`
`1039:                                      this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N + o_work, 1);`
`1040:                                      VMAX = ONE;`
`1041:                                      VCRIT = BIGNUM;`
`1042:                                  }`
`1043:                                  // *`
`1044:                                  WORK[J + N + o_work] = WORK[J + N + o_work] - this._ddot.Run(J - KI - 1, T, KI + 1+J * LDT + o_t, 1, WORK, KI + 1 + N + o_work, 1);`
`1045:                                  // *`
`1046:                                  WORK[J + 1 + N + o_work] = WORK[J + 1 + N + o_work] - this._ddot.Run(J - KI - 1, T, KI + 1+(J + 1) * LDT + o_t, 1, WORK, KI + 1 + N + o_work, 1);`
`1047:                                  // *`
`1048:                                  // *                    Solve`
`1049:                                  // *                      [T(J,J)-WR   T(J,J+1)     ]'* X = SCALE*( WORK1 )`
`1050:                                  // *                      [T(J+1,J)    T(J+1,J+1)-WR]             ( WORK2 )`
`1051:                                  // *`
`1052:                                  this._dlaln2.Run(true, 2, 1, SMIN, ONE, T, J+J * LDT + o_t`
`1053:                                                   , LDT, ONE, ONE, WORK, J + N + o_work, N, WR`
`1054:                                                   , ZERO, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);`
`1055:                                  // *`
`1056:                                  // *                    Scale if necessary`
`1057:                                  // *`
`1058:                                  if (SCALE != ONE) this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N + o_work, 1);`
`1059:                                  WORK[J + N + o_work] = X[1+1 * 2 + o_x];`
`1060:                                  WORK[J + 1 + N + o_work] = X[2+1 * 2 + o_x];`
`1061:                                  // *`
`1062:                                  VMAX = Math.Max(Math.Abs(WORK[J + N + o_work]), Math.Max(Math.Abs(WORK[J + 1 + N + o_work]), VMAX));`
`1063:                                  VCRIT = BIGNUM / VMAX;`
`1064:                                  // *`
`1065:                              }`
`1066:                          LABEL170:;`
`1067:                          }`
`1068:                          // *`
`1069:                          // *              Copy the vector x or Q*x to VL and normalize.`
`1070:                          // *`
`1071:                          if (!OVER)`
`1072:                          {`
`1073:                              this._dcopy.Run(N - KI + 1, WORK, KI + N + o_work, 1, ref VL, KI+IS * LDVL + o_vl, 1);`
`1074:                              // *`
`1075:                              II = this._idamax.Run(N - KI + 1, VL, KI+IS * LDVL + o_vl, 1) + KI - 1;`
`1076:                              REMAX = ONE / Math.Abs(VL[II+IS * LDVL + o_vl]);`
`1077:                              this._dscal.Run(N - KI + 1, REMAX, ref VL, KI+IS * LDVL + o_vl, 1);`
`1078:                              // *`
`1079:                              for (K = 1; K <= KI - 1; K++)`
`1080:                              {`
`1081:                                  VL[K+IS * LDVL + o_vl] = ZERO;`
`1082:                              }`
`1083:                              // *`
`1084:                          }`
`1085:                          else`
`1086:                          {`
`1087:                              // *`
`1088:                              if (KI < N)`
`1089:                              {`
`1090:                                  this._dgemv.Run("N", N, N - KI, ONE, VL, 1+(KI + 1) * LDVL + o_vl, LDVL`
`1091:                                                  , WORK, KI + 1 + N + o_work, 1, WORK[KI + N + o_work], ref VL, 1+KI * LDVL + o_vl, 1);`
`1092:                              }`
`1093:                              // *`
`1094:                              II = this._idamax.Run(N, VL, 1+KI * LDVL + o_vl, 1);`
`1095:                              REMAX = ONE / Math.Abs(VL[II+KI * LDVL + o_vl]);`
`1096:                              this._dscal.Run(N, REMAX, ref VL, 1+KI * LDVL + o_vl, 1);`
`1097:                              // *`
`1098:                          }`
`1099:                          // *`
`1100:                      }`
`1101:                      else`
`1102:                      {`
`1103:                          // *`
`1104:                          // *              Complex left eigenvector.`
`1105:                          // *`
`1106:                          // *               Initial solve:`
`1107:                          // *                 ((T(KI,KI)    T(KI,KI+1) )' - (WR - I* WI))*X = 0.`
`1108:                          // *                 ((T(KI+1,KI) T(KI+1,KI+1))                )`
`1109:                          // *`
`1110:                          if (Math.Abs(T[KI+(KI + 1) * LDT + o_t]) >= Math.Abs(T[KI + 1+KI * LDT + o_t]))`
`1111:                          {`
`1112:                              WORK[KI + N + o_work] = WI / T[KI+(KI + 1) * LDT + o_t];`
`1113:                              WORK[KI + 1 + N2 + o_work] = ONE;`
`1114:                          }`
`1115:                          else`
`1116:                          {`
`1117:                              WORK[KI + N + o_work] = ONE;`
`1118:                              WORK[KI + 1 + N2 + o_work] =  - WI / T[KI + 1+KI * LDT + o_t];`
`1119:                          }`
`1120:                          WORK[KI + 1 + N + o_work] = ZERO;`
`1121:                          WORK[KI + N2 + o_work] = ZERO;`
`1122:                          // *`
`1123:                          // *              Form right-hand side`
`1124:                          // *`
`1125:                          for (K = KI + 2; K <= N; K++)`
`1126:                          {`
`1127:                              WORK[K + N + o_work] =  - WORK[KI + N + o_work] * T[KI+K * LDT + o_t];`
`1128:                              WORK[K + N2 + o_work] =  - WORK[KI + 1 + N2 + o_work] * T[KI + 1+K * LDT + o_t];`
`1129:                          }`
`1130:                          // *`
`1131:                          // *              Solve complex quasi-triangular system:`
`1132:                          // *              ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2`
`1133:                          // *`
`1134:                          VMAX = ONE;`
`1135:                          VCRIT = BIGNUM;`
`1136:                          // *`
`1137:                          JNXT = KI + 2;`
`1138:                          for (J = KI + 2; J <= N; J++)`
`1139:                          {`
`1140:                              if (J < JNXT) goto LABEL200;`
`1141:                              J1 = J;`
`1142:                              J2 = J;`
`1143:                              JNXT = J + 1;`
`1144:                              if (J < N)`
`1145:                              {`
`1146:                                  if (T[J + 1+J * LDT + o_t] != ZERO)`
`1147:                                  {`
`1148:                                      J2 = J + 1;`
`1149:                                      JNXT = J + 2;`
`1150:                                  }`
`1151:                              }`
`1152:                              // *`
`1153:                              if (J1 == J2)`
`1154:                              {`
`1155:                                  // *`
`1156:                                  // *                    1-by-1 diagonal block`
`1157:                                  // *`
`1158:                                  // *                    Scale if necessary to avoid overflow when`
`1159:                                  // *                    forming the right-hand side elements.`
`1160:                                  // *`
`1161:                                  if (WORK[J + o_work] > VCRIT)`
`1162:                                  {`
`1163:                                      REC = ONE / VMAX;`
`1164:                                      this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N + o_work, 1);`
`1165:                                      this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N2 + o_work, 1);`
`1166:                                      VMAX = ONE;`
`1167:                                      VCRIT = BIGNUM;`
`1168:                                  }`
`1169:                                  // *`
`1170:                                  WORK[J + N + o_work] = WORK[J + N + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+J * LDT + o_t, 1, WORK, KI + 2 + N + o_work, 1);`
`1171:                                  WORK[J + N2 + o_work] = WORK[J + N2 + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+J * LDT + o_t, 1, WORK, KI + 2 + N2 + o_work, 1);`
`1172:                                  // *`
`1173:                                  // *                    Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2`
`1174:                                  // *`
`1175:                                  this._dlaln2.Run(false, 1, 2, SMIN, ONE, T, J+J * LDT + o_t`
`1176:                                                   , LDT, ONE, ONE, WORK, J + N + o_work, N, WR`
`1177:                                                   ,  - WI, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);`
`1178:                                  // *`
`1179:                                  // *                    Scale if necessary`
`1180:                                  // *`
`1181:                                  if (SCALE != ONE)`
`1182:                                  {`
`1183:                                      this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N + o_work, 1);`
`1184:                                      this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N2 + o_work, 1);`
`1185:                                  }`
`1186:                                  WORK[J + N + o_work] = X[1+1 * 2 + o_x];`
`1187:                                  WORK[J + N2 + o_work] = X[1+2 * 2 + o_x];`
`1188:                                  VMAX = Math.Max(Math.Abs(WORK[J + N + o_work]), Math.Max(Math.Abs(WORK[J + N2 + o_work]), VMAX));`
`1189:                                  VCRIT = BIGNUM / VMAX;`
`1190:                                  // *`
`1191:                              }`
`1192:                              else`
`1193:                              {`
`1194:                                  // *`
`1195:                                  // *                    2-by-2 diagonal block`
`1196:                                  // *`
`1197:                                  // *                    Scale if necessary to avoid overflow when forming`
`1198:                                  // *                    the right-hand side elements.`
`1199:                                  // *`
`1200:                                  BETA = Math.Max(WORK[J + o_work], WORK[J + 1 + o_work]);`
`1201:                                  if (BETA > VCRIT)`
`1202:                                  {`
`1203:                                      REC = ONE / VMAX;`
`1204:                                      this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N + o_work, 1);`
`1205:                                      this._dscal.Run(N - KI + 1, REC, ref WORK, KI + N2 + o_work, 1);`
`1206:                                      VMAX = ONE;`
`1207:                                      VCRIT = BIGNUM;`
`1208:                                  }`
`1209:                                  // *`
`1210:                                  WORK[J + N + o_work] = WORK[J + N + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+J * LDT + o_t, 1, WORK, KI + 2 + N + o_work, 1);`
`1211:                                  // *`
`1212:                                  WORK[J + N2 + o_work] = WORK[J + N2 + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+J * LDT + o_t, 1, WORK, KI + 2 + N2 + o_work, 1);`
`1213:                                  // *`
`1214:                                  WORK[J + 1 + N + o_work] = WORK[J + 1 + N + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+(J + 1) * LDT + o_t, 1, WORK, KI + 2 + N + o_work, 1);`
`1215:                                  // *`
`1216:                                  WORK[J + 1 + N2 + o_work] = WORK[J + 1 + N2 + o_work] - this._ddot.Run(J - KI - 2, T, KI + 2+(J + 1) * LDT + o_t, 1, WORK, KI + 2 + N2 + o_work, 1);`
`1217:                                  // *`
`1218:                                  // *                    Solve 2-by-2 complex linear equation`
`1219:                                  // *                      ([T(j,j)   T(j,j+1)  ]'-(wr-i*wi)*I)*X = SCALE*B`
`1220:                                  // *                      ([T(j+1,j) T(j+1,j+1)]             )`
`1221:                                  // *`
`1222:                                  this._dlaln2.Run(true, 2, 2, SMIN, ONE, T, J+J * LDT + o_t`
`1223:                                                   , LDT, ONE, ONE, WORK, J + N + o_work, N, WR`
`1224:                                                   ,  - WI, ref X, offset_x, 2, ref SCALE, ref XNORM, ref IERR);`
`1225:                                  // *`
`1226:                                  // *                    Scale if necessary`
`1227:                                  // *`
`1228:                                  if (SCALE != ONE)`
`1229:                                  {`
`1230:                                      this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N + o_work, 1);`
`1231:                                      this._dscal.Run(N - KI + 1, SCALE, ref WORK, KI + N2 + o_work, 1);`
`1232:                                  }`
`1233:                                  WORK[J + N + o_work] = X[1+1 * 2 + o_x];`
`1234:                                  WORK[J + N2 + o_work] = X[1+2 * 2 + o_x];`
`1235:                                  WORK[J + 1 + N + o_work] = X[2+1 * 2 + o_x];`
`1236:                                  WORK[J + 1 + N2 + o_work] = X[2+2 * 2 + o_x];`
`1237:                                  VMAX = Math.Max(Math.Abs(X[1+1 * 2 + o_x]), Math.Max(Math.Abs(X[1+2 * 2 + o_x]), Math.Max(Math.Abs(X[2+1 * 2 + o_x]), Math.Max(Math.Abs(X[2+2 * 2 + o_x]), VMAX))));`
`1238:                                  VCRIT = BIGNUM / VMAX;`
`1239:                                  // *`
`1240:                              }`
`1241:                          LABEL200:;`
`1242:                          }`
`1243:                          // *`
`1244:                          // *              Copy the vector x or Q*x to VL and normalize.`
`1245:                          // *`
`1246:                          if (!OVER)`
`1247:                          {`
`1248:                              this._dcopy.Run(N - KI + 1, WORK, KI + N + o_work, 1, ref VL, KI+IS * LDVL + o_vl, 1);`
`1249:                              this._dcopy.Run(N - KI + 1, WORK, KI + N2 + o_work, 1, ref VL, KI+(IS + 1) * LDVL + o_vl, 1);`
`1250:                              // *`
`1251:                              EMAX = ZERO;`
`1252:                              for (K = KI; K <= N; K++)`
`1253:                              {`
`1254:                                  EMAX = Math.Max(EMAX, Math.Abs(VL[K+IS * LDVL + o_vl]) + Math.Abs(VL[K+(IS + 1) * LDVL + o_vl]));`
`1255:                              }`
`1256:                              REMAX = ONE / EMAX;`
`1257:                              this._dscal.Run(N - KI + 1, REMAX, ref VL, KI+IS * LDVL + o_vl, 1);`
`1258:                              this._dscal.Run(N - KI + 1, REMAX, ref VL, KI+(IS + 1) * LDVL + o_vl, 1);`
`1259:                              // *`
`1260:                              for (K = 1; K <= KI - 1; K++)`
`1261:                              {`
`1262:                                  VL[K+IS * LDVL + o_vl] = ZERO;`
`1263:                                  VL[K+(IS + 1) * LDVL + o_vl] = ZERO;`
`1264:                              }`
`1265:                          }`
`1266:                          else`
`1267:                          {`
`1268:                              if (KI < N - 1)`
`1269:                              {`
`1270:                                  this._dgemv.Run("N", N, N - KI - 1, ONE, VL, 1+(KI + 2) * LDVL + o_vl, LDVL`
`1271:                                                  , WORK, KI + 2 + N + o_work, 1, WORK[KI + N + o_work], ref VL, 1+KI * LDVL + o_vl, 1);`
`1272:                                  this._dgemv.Run("N", N, N - KI - 1, ONE, VL, 1+(KI + 2) * LDVL + o_vl, LDVL`
`1273:                                                  , WORK, KI + 2 + N2 + o_work, 1, WORK[KI + 1 + N2 + o_work], ref VL, 1+(KI + 1) * LDVL + o_vl, 1);`
`1274:                              }`
`1275:                              else`
`1276:                              {`
`1277:                                  this._dscal.Run(N, WORK[KI + N + o_work], ref VL, 1+KI * LDVL + o_vl, 1);`
`1278:                                  this._dscal.Run(N, WORK[KI + 1 + N2 + o_work], ref VL, 1+(KI + 1) * LDVL + o_vl, 1);`
`1279:                              }`
`1280:                              // *`
`1281:                              EMAX = ZERO;`
`1282:                              for (K = 1; K <= N; K++)`
`1283:                              {`
`1284:                                  EMAX = Math.Max(EMAX, Math.Abs(VL[K+KI * LDVL + o_vl]) + Math.Abs(VL[K+(KI + 1) * LDVL + o_vl]));`
`1285:                              }`
`1286:                              REMAX = ONE / EMAX;`
`1287:                              this._dscal.Run(N, REMAX, ref VL, 1+KI * LDVL + o_vl, 1);`
`1288:                              this._dscal.Run(N, REMAX, ref VL, 1+(KI + 1) * LDVL + o_vl, 1);`
`1289:                              // *`
`1290:                          }`
`1291:                          // *`
`1292:                      }`
`1293:                      // *`
`1294:                      IS = IS + 1;`
`1295:                      if (IP != 0) IS = IS + 1;`
`1296:                  LABEL250:;`
`1297:                      if (IP ==  - 1) IP = 0;`
`1298:                      if (IP == 1) IP =  - 1;`
`1299:                      // *`
`1300:                  }`
`1301:                  // *`
`1302:              }`
`1303:              // *`
`1304:              return;`
`1305:              // *`
`1306:              // *     End of DTREVC`
`1307:              // *`
`1308:   `
`1309:              #endregion`
`1310:   `
`1311:          }`
`1312:      }`
`1313:  }`