Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   1:  #region Translated by Jose Antonio De Santiago-Castillo.
   2:   
   3:  //Translated by Jose Antonio De Santiago-Castillo. 
   4:  //E-mail:JAntonioDeSantiago@gmail.com
   5:  //Web: www.DotNumerics.com
   6:  //
   7:  //Fortran to C# Translation.
   8:  //Translated by:
   9:  //F2CSharp Version 0.71 (November 10, 2009)
  10:  //Code Optimizations: None
  11:  //
  12:  #endregion
  13:   
  14:  using System;
  15:  using DotNumerics.FortranLibrary;
  16:   
  17:  namespace DotNumerics.CSLapack
  18:  {
  19:      /// <summary>
  20:      /// -- LAPACK routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// 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:  }