Skip Navigation Links
Numerical Libraries
Linear Algebra
Differential Equations
Optimization
Samples
Skip Navigation Links
Linear Algebra
CSLapack
CSBlas
   1:  #region Translated by Jose Antonio De Santiago-Castillo.
   2:   
   3:  //Translated by Jose Antonio De Santiago-Castillo. 
   4:  //E-mail:JAntonioDeSantiago@gmail.com
   5:  //Web: www.DotNumerics.com
   6:  //
   7:  //Fortran to C# Translation.
   8:  //Translated by:
   9:  //F2CSharp Version 0.71 (November 10, 2009)
  10:  //Code Optimizations: None
  11:  //
  12:  #endregion
  13:   
  14:  using System;
  15:  using DotNumerics.FortranLibrary;
  16:   
  17:  namespace DotNumerics.CSLapack
  18:  {
  19:      /// <summary>
  20:      /// -- LAPACK driver routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DGEEV computes for an N-by-N real nonsymmetric matrix A, the
  27:      /// eigenvalues and, optionally, the left and/or right eigenvectors.
  28:      /// 
  29:      /// The right eigenvector v(j) of A satisfies
  30:      /// A * v(j) = lambda(j) * v(j)
  31:      /// where lambda(j) is its eigenvalue.
  32:      /// The left eigenvector u(j) of A satisfies
  33:      /// u(j)**H * A = lambda(j) * u(j)**H
  34:      /// where u(j)**H denotes the conjugate transpose of u(j).
  35:      /// 
  36:      /// The computed eigenvectors are normalized to have Euclidean norm
  37:      /// equal to 1 and largest component real.
  38:      /// 
  39:      ///</summary>
  40:      public class DGEEV
  41:      {
  42:      
  43:   
  44:          #region Dependencies
  45:          
  46:          DGEBAK _dgebak; DGEBAL _dgebal; DGEHRD _dgehrd; DHSEQR _dhseqr; DLABAD _dlabad; DLACPY _dlacpy; DLARTG _dlartg; 
  47:          DLASCL _dlascl;DORGHR _dorghr; DROT _drot; DSCAL _dscal; DTREVC _dtrevc; XERBLA _xerbla; LSAME _lsame; IDAMAX _idamax; 
  48:          ILAENV _ilaenv;DLAMCH _dlamch; DLANGE _dlange; DLAPY2 _dlapy2; DNRM2 _dnrm2; 
  49:   
  50:          #endregion
  51:   
  52:   
  53:          #region Fields
  54:          
  55:          const double ZERO = 0.0E0; const double ONE = 1.0E0; bool LQUERY = false; bool SCALEA = false; bool WANTVL = false; 
  56:          bool WANTVR = false;string SIDE = new string(' ', 1); int HSWORK = 0; int I = 0; int IBAL = 0; int IERR = 0; int IHI = 0; 
  57:          int ILO = 0;int ITAU = 0; int IWRK = 0; int K = 0; int MAXWRK = 0; int MINWRK = 0; int NOUT = 0; double ANRM = 0; 
  58:          double BIGNUM = 0;double CS = 0; double CSCALE = 0; double EPS = 0; double R = 0; double SCL = 0; double SMLNUM = 0; 
  59:          double SN = 0;bool[] SELECT = new bool[1]; int offset_select = 0;
  60:          double[] DUM = new double[1]; int offset_dum = 0;
  61:   
  62:          #endregion
  63:   
  64:          public DGEEV(DGEBAK dgebak, DGEBAL dgebal, DGEHRD dgehrd, DHSEQR dhseqr, DLABAD dlabad, DLACPY dlacpy, DLARTG dlartg, DLASCL dlascl, DORGHR dorghr, DROT drot
  65:                       , DSCAL dscal, DTREVC dtrevc, XERBLA xerbla, LSAME lsame, IDAMAX idamax, ILAENV ilaenv, DLAMCH dlamch, DLANGE dlange, DLAPY2 dlapy2, DNRM2 dnrm2)
  66:          {
  67:      
  68:   
  69:              #region Set Dependencies
  70:              
  71:              this._dgebak = dgebak; this._dgebal = dgebal; this._dgehrd = dgehrd; this._dhseqr = dhseqr; this._dlabad = dlabad; 
  72:              this._dlacpy = dlacpy;this._dlartg = dlartg; this._dlascl = dlascl; this._dorghr = dorghr; this._drot = drot; 
  73:              this._dscal = dscal;this._dtrevc = dtrevc; this._xerbla = xerbla; this._lsame = lsame; this._idamax = idamax; 
  74:              this._ilaenv = ilaenv;this._dlamch = dlamch; this._dlange = dlange; this._dlapy2 = dlapy2; this._dnrm2 = dnrm2; 
  75:   
  76:              #endregion
  77:   
  78:          }
  79:      
  80:          public DGEEV()
  81:          {
  82:      
  83:   
  84:              #region Dependencies (Initialization)
  85:              
  86:              LSAME lsame = new LSAME();
  87:              DSCAL dscal = new DSCAL();
  88:              DSWAP dswap = new DSWAP();
  89:              XERBLA xerbla = new XERBLA();
  90:              IDAMAX idamax = new IDAMAX();
  91:              DLAMC3 dlamc3 = new DLAMC3();
  92:              DAXPY daxpy = new DAXPY();
  93:              DLAPY2 dlapy2 = new DLAPY2();
  94:              DNRM2 dnrm2 = new DNRM2();
  95:              DCOPY dcopy = new DCOPY();
  96:              IEEECK ieeeck = new IEEECK();
  97:              IPARMQ iparmq = new IPARMQ();
  98:              DLABAD dlabad = new DLABAD();
  99:              DROT drot = new DROT();
 100:              DLASSQ dlassq = new DLASSQ();
 101:              DLAQR1 dlaqr1 = new DLAQR1();
 102:              DDOT ddot = new DDOT();
 103:              DLADIV dladiv = new DLADIV();
 104:              DGEBAK dgebak = new DGEBAK(lsame, dscal, dswap, xerbla);
 105:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
 106:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 107:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 108:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 109:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 110:              DGEBAL dgebal = new DGEBAL(lsame, idamax, dlamch, dscal, dswap, xerbla);
 111:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 112:              DGER dger = new DGER(xerbla);
 113:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
 114:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
 115:              DGEHD2 dgehd2 = new DGEHD2(dlarf, dlarfg, xerbla);
 116:              DGEMM dgemm = new DGEMM(lsame, xerbla);
 117:              DLACPY dlacpy = new DLACPY(lsame);
 118:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
 119:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
 120:              DLAHR2 dlahr2 = new DLAHR2(daxpy, dcopy, dgemm, dgemv, dlacpy, dlarfg, dscal, dtrmm, dtrmv);
 121:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
 122:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
 123:              DGEHRD dgehrd = new DGEHRD(daxpy, dgehd2, dgemm, dlahr2, dlarfb, dtrmm, xerbla, ilaenv);
 124:              DLANV2 dlanv2 = new DLANV2(dlamch, dlapy2);
 125:              DLAHQR dlahqr = new DLAHQR(dlamch, dcopy, dlabad, dlanv2, dlarfg, drot);
 126:              DLASET dlaset = new DLASET(lsame);
 127:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
 128:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
 129:              DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
 130:              DORGHR dorghr = new DORGHR(dorgqr, xerbla, ilaenv);
 131:              DLANGE dlange = new DLANGE(dlassq, lsame);
 132:              DLARFX dlarfx = new DLARFX(lsame, dgemv, dger);
 133:              DLARTG dlartg = new DLARTG(dlamch);
 134:              DLASY2 dlasy2 = new DLASY2(idamax, dlamch, dcopy, dswap);
 135:              DLAEXC dlaexc = new DLAEXC(dlamch, dlange, dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2, drot);
 136:              DTREXC dtrexc = new DTREXC(lsame, dlaexc, xerbla);
 137:              DLAQR2 dlaqr2 = new DLAQR2(dlamch, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlarf, dlarfg
 138:                                         , dlaset, dorghr, dtrexc);
 139:              DLAQR5 dlaqr5 = new DLAQR5(dlamch, dgemm, dlabad, dlacpy, dlaqr1, dlarfg, dlaset, dtrmm);
 140:              DLAQR4 dlaqr4 = new DLAQR4(ilaenv, dlacpy, dlahqr, dlanv2, dlaqr2, dlaqr5);
 141:              DLAQR3 dlaqr3 = new DLAQR3(dlamch, ilaenv, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlaqr4
 142:                                         , dlarf, dlarfg, dlaset, dorghr, dtrexc);
 143:              DLAQR0 dlaqr0 = new DLAQR0(ilaenv, dlacpy, dlahqr, dlanv2, dlaqr3, dlaqr4, dlaqr5);
 144:              DHSEQR dhseqr = new DHSEQR(ilaenv, lsame, dlacpy, dlahqr, dlaqr0, dlaset, xerbla);
 145:              DLASCL dlascl = new DLASCL(lsame, dlamch, xerbla);
 146:              DLALN2 dlaln2 = new DLALN2(dlamch, dladiv);
 147:              DTREVC dtrevc = new DTREVC(lsame, idamax, ddot, dlamch, daxpy, dcopy, dgemv, dlaln2, dscal, xerbla
 148:                                         , dlabad);
 149:   
 150:              #endregion
 151:   
 152:   
 153:              #region Set Dependencies
 154:              
 155:              this._dgebak = dgebak; this._dgebal = dgebal; this._dgehrd = dgehrd; this._dhseqr = dhseqr; this._dlabad = dlabad; 
 156:              this._dlacpy = dlacpy;this._dlartg = dlartg; this._dlascl = dlascl; this._dorghr = dorghr; this._drot = drot; 
 157:              this._dscal = dscal;this._dtrevc = dtrevc; this._xerbla = xerbla; this._lsame = lsame; this._idamax = idamax; 
 158:              this._ilaenv = ilaenv;this._dlamch = dlamch; this._dlange = dlange; this._dlapy2 = dlapy2; this._dnrm2 = dnrm2; 
 159:   
 160:              #endregion
 161:   
 162:          }
 163:          /// <summary>
 164:          /// Purpose
 165:          /// =======
 166:          /// 
 167:          /// DGEEV computes for an N-by-N real nonsymmetric matrix A, the
 168:          /// eigenvalues and, optionally, the left and/or right eigenvectors.
 169:          /// 
 170:          /// The right eigenvector v(j) of A satisfies
 171:          /// A * v(j) = lambda(j) * v(j)
 172:          /// where lambda(j) is its eigenvalue.
 173:          /// The left eigenvector u(j) of A satisfies
 174:          /// u(j)**H * A = lambda(j) * u(j)**H
 175:          /// where u(j)**H denotes the conjugate transpose of u(j).
 176:          /// 
 177:          /// The computed eigenvectors are normalized to have Euclidean norm
 178:          /// equal to 1 and largest component real.
 179:          /// 
 180:          ///</summary>
 181:          /// <param name="JOBVL">
 182:          /// (input) CHARACTER*1
 183:          /// = 'N': left eigenvectors of A are not computed;
 184:          /// = 'V': left eigenvectors of A are computed.
 185:          ///</param>
 186:          /// <param name="JOBVR">
 187:          /// (input) CHARACTER*1
 188:          /// = 'N': right eigenvectors of A are not computed;
 189:          /// = 'V': right eigenvectors of A are computed.
 190:          ///</param>
 191:          /// <param name="N">
 192:          /// (input) INTEGER
 193:          /// The order of the matrix A. N .GE. 0.
 194:          ///</param>
 195:          /// <param name="A">
 196:          /// (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 197:          /// On entry, the N-by-N matrix A.
 198:          /// On exit, A has been overwritten.
 199:          ///</param>
 200:          /// <param name="LDA">
 201:          /// (input) INTEGER
 202:          /// The leading dimension of the array A.  LDA .GE. max(1,N).
 203:          ///</param>
 204:          /// <param name="WR">
 205:          /// (output) DOUBLE PRECISION array, dimension (N)
 206:          ///</param>
 207:          /// <param name="WI">
 208:          /// (output) DOUBLE PRECISION array, dimension (N)
 209:          /// WR and WI contain the real and imaginary parts,
 210:          /// respectively, of the computed eigenvalues.  Complex
 211:          /// conjugate pairs of eigenvalues appear consecutively
 212:          /// with the eigenvalue having the positive imaginary part
 213:          /// first.
 214:          ///</param>
 215:          /// <param name="VL">
 216:          /// (output) DOUBLE PRECISION array, dimension (LDVL,N)
 217:          /// If JOBVL = 'V', the left eigenvectors u(j) are stored one
 218:          /// after another in the columns of VL, in the same order
 219:          /// as their eigenvalues.
 220:          /// If JOBVL = 'N', VL is not referenced.
 221:          /// If the j-th eigenvalue is real, then u(j) = VL(:,j),
 222:          /// the j-th column of VL.
 223:          /// If the j-th and (j+1)-st eigenvalues form a complex
 224:          /// conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
 225:          /// u(j+1) = VL(:,j) - i*VL(:,j+1).
 226:          ///</param>
 227:          /// <param name="LDVL">
 228:          /// (input) INTEGER
 229:          /// The leading dimension of the array VL.  LDVL .GE. 1; if
 230:          /// JOBVL = 'V', LDVL .GE. N.
 231:          ///</param>
 232:          /// <param name="VR">
 233:          /// (output) DOUBLE PRECISION array, dimension (LDVR,N)
 234:          /// If JOBVR = 'V', the right eigenvectors v(j) are stored one
 235:          /// after another in the columns of VR, in the same order
 236:          /// as their eigenvalues.
 237:          /// If JOBVR = 'N', VR is not referenced.
 238:          /// If the j-th eigenvalue is real, then v(j) = VR(:,j),
 239:          /// the j-th column of VR.
 240:          /// If the j-th and (j+1)-st eigenvalues form a complex
 241:          /// conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
 242:          /// v(j+1) = VR(:,j) - i*VR(:,j+1).
 243:          ///</param>
 244:          /// <param name="LDVR">
 245:          /// (input) INTEGER
 246:          /// The leading dimension of the array VR.  LDVR .GE. 1; if
 247:          /// JOBVR = 'V', LDVR .GE. N.
 248:          ///</param>
 249:          /// <param name="WORK">
 250:          /// (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 251:          /// On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 252:          ///</param>
 253:          /// <param name="LWORK">
 254:          /// (input) INTEGER
 255:          /// The dimension of the array WORK.  LWORK .GE. max(1,3*N), and
 256:          /// if JOBVL = 'V' or JOBVR = 'V', LWORK .GE. 4*N.  For good
 257:          /// performance, LWORK must generally be larger.
 258:          /// 
 259:          /// If LWORK = -1, then a workspace query is assumed; the routine
 260:          /// only calculates the optimal size of the WORK array, returns
 261:          /// this value as the first entry of the WORK array, and no error
 262:          /// message related to LWORK is issued by XERBLA.
 263:          ///</param>
 264:          /// <param name="INFO">
 265:          /// (output) INTEGER
 266:          /// = 0:  successful exit
 267:          /// .LT. 0:  if INFO = -i, the i-th argument had an illegal value.
 268:          /// .GT. 0:  if INFO = i, the QR algorithm failed to compute all the
 269:          /// eigenvalues, and no eigenvectors have been computed;
 270:          /// elements i+1:N of WR and WI contain eigenvalues which
 271:          /// have converged.
 272:          ///</param>
 273:          public void Run(string JOBVL, string JOBVR, int N, ref double[] A, int offset_a, int LDA, ref double[] WR, int offset_wr
 274:                           , ref double[] WI, int offset_wi, ref double[] VL, int offset_vl, int LDVL, ref double[] VR, int offset_vr, int LDVR, ref double[] WORK, int offset_work
 275:                           , int LWORK, ref int INFO)
 276:          {
 277:   
 278:              #region Array Index Correction
 279:              
 280:               int o_a = -1 - LDA + offset_a;  int o_wr = -1 + offset_wr;  int o_wi = -1 + offset_wi; 
 281:               int o_vl = -1 - LDVL + offset_vl; int o_vr = -1 - LDVR + offset_vr;  int o_work = -1 + offset_work; 
 282:   
 283:              #endregion
 284:   
 285:   
 286:              #region Strings
 287:              
 288:              JOBVL = JOBVL.Substring(0, 1);  JOBVR = JOBVR.Substring(0, 1);  
 289:   
 290:              #endregion
 291:   
 292:   
 293:              #region Prolog
 294:              
 295:              // *
 296:              // *  -- LAPACK driver routine (version 3.1) --
 297:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 298:              // *     November 2006
 299:              // *
 300:              // *     .. Scalar Arguments ..
 301:              // *     ..
 302:              // *     .. Array Arguments ..
 303:              // *     ..
 304:              // *
 305:              // *  Purpose
 306:              // *  =======
 307:              // *
 308:              // *  DGEEV computes for an N-by-N real nonsymmetric matrix A, the
 309:              // *  eigenvalues and, optionally, the left and/or right eigenvectors.
 310:              // *
 311:              // *  The right eigenvector v(j) of A satisfies
 312:              // *                   A * v(j) = lambda(j) * v(j)
 313:              // *  where lambda(j) is its eigenvalue.
 314:              // *  The left eigenvector u(j) of A satisfies
 315:              // *                u(j)**H * A = lambda(j) * u(j)**H
 316:              // *  where u(j)**H denotes the conjugate transpose of u(j).
 317:              // *
 318:              // *  The computed eigenvectors are normalized to have Euclidean norm
 319:              // *  equal to 1 and largest component real.
 320:              // *
 321:              // *  Arguments
 322:              // *  =========
 323:              // *
 324:              // *  JOBVL   (input) CHARACTER*1
 325:              // *          = 'N': left eigenvectors of A are not computed;
 326:              // *          = 'V': left eigenvectors of A are computed.
 327:              // *
 328:              // *  JOBVR   (input) CHARACTER*1
 329:              // *          = 'N': right eigenvectors of A are not computed;
 330:              // *          = 'V': right eigenvectors of A are computed.
 331:              // *
 332:              // *  N       (input) INTEGER
 333:              // *          The order of the matrix A. N >= 0.
 334:              // *
 335:              // *  A       (input/output) DOUBLE PRECISION array, dimension (LDA,N)
 336:              // *          On entry, the N-by-N matrix A.
 337:              // *          On exit, A has been overwritten.
 338:              // *
 339:              // *  LDA     (input) INTEGER
 340:              // *          The leading dimension of the array A.  LDA >= max(1,N).
 341:              // *
 342:              // *  WR      (output) DOUBLE PRECISION array, dimension (N)
 343:              // *  WI      (output) DOUBLE PRECISION array, dimension (N)
 344:              // *          WR and WI contain the real and imaginary parts,
 345:              // *          respectively, of the computed eigenvalues.  Complex
 346:              // *          conjugate pairs of eigenvalues appear consecutively
 347:              // *          with the eigenvalue having the positive imaginary part
 348:              // *          first.
 349:              // *
 350:              // *  VL      (output) DOUBLE PRECISION array, dimension (LDVL,N)
 351:              // *          If JOBVL = 'V', the left eigenvectors u(j) are stored one
 352:              // *          after another in the columns of VL, in the same order
 353:              // *          as their eigenvalues.
 354:              // *          If JOBVL = 'N', VL is not referenced.
 355:              // *          If the j-th eigenvalue is real, then u(j) = VL(:,j),
 356:              // *          the j-th column of VL.
 357:              // *          If the j-th and (j+1)-st eigenvalues form a complex
 358:              // *          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
 359:              // *          u(j+1) = VL(:,j) - i*VL(:,j+1).
 360:              // *
 361:              // *  LDVL    (input) INTEGER
 362:              // *          The leading dimension of the array VL.  LDVL >= 1; if
 363:              // *          JOBVL = 'V', LDVL >= N.
 364:              // *
 365:              // *  VR      (output) DOUBLE PRECISION array, dimension (LDVR,N)
 366:              // *          If JOBVR = 'V', the right eigenvectors v(j) are stored one
 367:              // *          after another in the columns of VR, in the same order
 368:              // *          as their eigenvalues.
 369:              // *          If JOBVR = 'N', VR is not referenced.
 370:              // *          If the j-th eigenvalue is real, then v(j) = VR(:,j),
 371:              // *          the j-th column of VR.
 372:              // *          If the j-th and (j+1)-st eigenvalues form a complex
 373:              // *          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
 374:              // *          v(j+1) = VR(:,j) - i*VR(:,j+1).
 375:              // *
 376:              // *  LDVR    (input) INTEGER
 377:              // *          The leading dimension of the array VR.  LDVR >= 1; if
 378:              // *          JOBVR = 'V', LDVR >= N.
 379:              // *
 380:              // *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 381:              // *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 382:              // *
 383:              // *  LWORK   (input) INTEGER
 384:              // *          The dimension of the array WORK.  LWORK >= max(1,3*N), and
 385:              // *          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
 386:              // *          performance, LWORK must generally be larger.
 387:              // *
 388:              // *          If LWORK = -1, then a workspace query is assumed; the routine
 389:              // *          only calculates the optimal size of the WORK array, returns
 390:              // *          this value as the first entry of the WORK array, and no error
 391:              // *          message related to LWORK is issued by XERBLA.
 392:              // *
 393:              // *  INFO    (output) INTEGER
 394:              // *          = 0:  successful exit
 395:              // *          < 0:  if INFO = -i, the i-th argument had an illegal value.
 396:              // *          > 0:  if INFO = i, the QR algorithm failed to compute all the
 397:              // *                eigenvalues, and no eigenvectors have been computed;
 398:              // *                elements i+1:N of WR and WI contain eigenvalues which
 399:              // *                have converged.
 400:              // *
 401:              // *  =====================================================================
 402:              // *
 403:              // *     .. Parameters ..
 404:              // *     ..
 405:              // *     .. Local Scalars ..
 406:              // *     ..
 407:              // *     .. Local Arrays ..
 408:              // *     ..
 409:              // *     .. External Subroutines ..
 410:              // *     ..
 411:              // *     .. External Functions ..
 412:              // *     ..
 413:              // *     .. Intrinsic Functions ..
 414:              //      INTRINSIC          MAX, SQRT;
 415:              // *     ..
 416:              // *     .. Executable Statements ..
 417:              // *
 418:              // *     Test the input arguments
 419:              // *
 420:   
 421:              #endregion
 422:   
 423:   
 424:              #region Body
 425:              
 426:              INFO = 0;
 427:              LQUERY = (LWORK ==  - 1);
 428:              WANTVL = this._lsame.Run(JOBVL, "V");
 429:              WANTVR = this._lsame.Run(JOBVR, "V");
 430:              if ((!WANTVL) && (!this._lsame.Run(JOBVL, "N")))
 431:              {
 432:                  INFO =  - 1;
 433:              }
 434:              else
 435:              {
 436:                  if ((!WANTVR) && (!this._lsame.Run(JOBVR, "N")))
 437:                  {
 438:                      INFO =  - 2;
 439:                  }
 440:                  else
 441:                  {
 442:                      if (N < 0)
 443:                      {
 444:                          INFO =  - 3;
 445:                      }
 446:                      else
 447:                      {
 448:                          if (LDA < Math.Max(1, N))
 449:                          {
 450:                              INFO =  - 5;
 451:                          }
 452:                          else
 453:                          {
 454:                              if (LDVL < 1 || (WANTVL && LDVL < N))
 455:                              {
 456:                                  INFO =  - 9;
 457:                              }
 458:                              else
 459:                              {
 460:                                  if (LDVR < 1 || (WANTVR && LDVR < N))
 461:                                  {
 462:                                      INFO =  - 11;
 463:                                  }
 464:                              }
 465:                          }
 466:                      }
 467:                  }
 468:              }
 469:              // *
 470:              // *     Compute workspace
 471:              // *      (Note: Comments in the code beginning "Workspace:" describe the
 472:              // *       minimal amount of workspace needed at that point in the code,
 473:              // *       as well as the preferred amount for good performance.
 474:              // *       NB refers to the optimal block size for the immediately
 475:              // *       following subroutine, as returned by ILAENV.
 476:              // *       HSWORK refers to the workspace preferred by DHSEQR, as
 477:              // *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
 478:              // *       the worst case.)
 479:              // *
 480:              if (INFO == 0)
 481:              {
 482:                  if (N == 0)
 483:                  {
 484:                      MINWRK = 1;
 485:                      MAXWRK = 1;
 486:                  }
 487:                  else
 488:                  {
 489:                      MAXWRK = 2 * N + N * this._ilaenv.Run(1, "DGEHRD", " ", N, 1, N, 0);
 490:                      if (WANTVL)
 491:                      {
 492:                          MINWRK = 4 * N;
 493:                          MAXWRK = Math.Max(MAXWRK, 2 * N + (N - 1) * this._ilaenv.Run(1, "DORGHR", " ", N, 1, N,  - 1));
 494:                          this._dhseqr.Run("S", "V", N, 1, N, ref A, offset_a
 495:                                           , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VL, offset_vl, LDVL, ref WORK, offset_work
 496:                                           ,  - 1, ref INFO);
 497:                          HSWORK = (int)WORK[1 + o_work];
 498:                          MAXWRK = Math.Max(MAXWRK, Math.Max(N + 1, N + HSWORK));
 499:                          MAXWRK = Math.Max(MAXWRK, 4 * N);
 500:                      }
 501:                      else
 502:                      {
 503:                          if (WANTVR)
 504:                          {
 505:                              MINWRK = 4 * N;
 506:                              MAXWRK = Math.Max(MAXWRK, 2 * N + (N - 1) * this._ilaenv.Run(1, "DORGHR", " ", N, 1, N,  - 1));
 507:                              this._dhseqr.Run("S", "V", N, 1, N, ref A, offset_a
 508:                                               , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VR, offset_vr, LDVR, ref WORK, offset_work
 509:                                               ,  - 1, ref INFO);
 510:                              HSWORK = (int)WORK[1 + o_work];
 511:                              MAXWRK = Math.Max(MAXWRK, Math.Max(N + 1, N + HSWORK));
 512:                              MAXWRK = Math.Max(MAXWRK, 4 * N);
 513:                          }
 514:                          else
 515:                          {
 516:                              MINWRK = 3 * N;
 517:                              this._dhseqr.Run("E", "N", N, 1, N, ref A, offset_a
 518:                                               , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VR, offset_vr, LDVR, ref WORK, offset_work
 519:                                               ,  - 1, ref INFO);
 520:                              HSWORK = (int)WORK[1 + o_work];
 521:                              MAXWRK = Math.Max(MAXWRK, Math.Max(N + 1, N + HSWORK));
 522:                          }
 523:                      }
 524:                      MAXWRK = Math.Max(MAXWRK, MINWRK);
 525:                  }
 526:                  WORK[1 + o_work] = MAXWRK;
 527:                  // *
 528:                  if (LWORK < MINWRK && !LQUERY)
 529:                  {
 530:                      INFO =  - 13;
 531:                  }
 532:              }
 533:              // *
 534:              if (INFO != 0)
 535:              {
 536:                  this._xerbla.Run("DGEEV ",  - INFO);
 537:                  return;
 538:              }
 539:              else
 540:              {
 541:                  if (LQUERY)
 542:                  {
 543:                      return;
 544:                  }
 545:              }
 546:              // *
 547:              // *     Quick return if possible
 548:              // *
 549:              if (N == 0) return;
 550:              // *
 551:              // *     Get machine constants
 552:              // *
 553:              EPS = this._dlamch.Run("P");
 554:              SMLNUM = this._dlamch.Run("S");
 555:              BIGNUM = ONE / SMLNUM;
 556:              this._dlabad.Run(ref SMLNUM, ref BIGNUM);
 557:              SMLNUM = Math.Sqrt(SMLNUM) / EPS;
 558:              BIGNUM = ONE / SMLNUM;
 559:              // *
 560:              // *     Scale A if max element outside range [SMLNUM,BIGNUM]
 561:              // *
 562:              ANRM = this._dlange.Run("M", N, N, A, offset_a, LDA, ref DUM, offset_dum);
 563:              SCALEA = false;
 564:              if (ANRM > ZERO && ANRM < SMLNUM)
 565:              {
 566:                  SCALEA = true;
 567:                  CSCALE = SMLNUM;
 568:              }
 569:              else
 570:              {
 571:                  if (ANRM > BIGNUM)
 572:                  {
 573:                      SCALEA = true;
 574:                      CSCALE = BIGNUM;
 575:                  }
 576:              }
 577:              if (SCALEA)
 578:              {
 579:                  this._dlascl.Run("G", 0, 0, ANRM, CSCALE, N
 580:                                   , N, ref A, offset_a, LDA, ref IERR);
 581:              }
 582:              // *
 583:              // *     Balance the matrix
 584:              // *     (Workspace: need N)
 585:              // *
 586:              IBAL = 1;
 587:              this._dgebal.Run("B", N, ref A, offset_a, LDA, ref ILO, ref IHI
 588:                               , ref WORK, IBAL + o_work, ref IERR);
 589:              // *
 590:              // *     Reduce to upper Hessenberg form
 591:              // *     (Workspace: need 3*N, prefer 2*N+N*NB)
 592:              // *
 593:              ITAU = IBAL + N;
 594:              IWRK = ITAU + N;
 595:              this._dgehrd.Run(N, ILO, IHI, ref A, offset_a, LDA, ref WORK, ITAU + o_work
 596:                               , ref WORK, IWRK + o_work, LWORK - IWRK + 1, ref IERR);
 597:              // *
 598:              if (WANTVL)
 599:              {
 600:                  // *
 601:                  // *        Want left eigenvectors
 602:                  // *        Copy Householder vectors to VL
 603:                  // *
 604:                  FortranLib.Copy(ref SIDE , "L");
 605:                  this._dlacpy.Run("L", N, N, A, offset_a, LDA, ref VL, offset_vl
 606:                                   , LDVL);
 607:                  // *
 608:                  // *        Generate orthogonal matrix in VL
 609:                  // *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
 610:                  // *
 611:                  this._dorghr.Run(N, ILO, IHI, ref VL, offset_vl, LDVL, WORK, ITAU + o_work
 612:                                   , ref WORK, IWRK + o_work, LWORK - IWRK + 1, ref IERR);
 613:                  // *
 614:                  // *        Perform QR iteration, accumulating Schur vectors in VL
 615:                  // *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
 616:                  // *
 617:                  IWRK = ITAU;
 618:                  this._dhseqr.Run("S", "V", N, ILO, IHI, ref A, offset_a
 619:                                   , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VL, offset_vl, LDVL, ref WORK, IWRK + o_work
 620:                                   , LWORK - IWRK + 1, ref INFO);
 621:                  // *
 622:                  if (WANTVR)
 623:                  {
 624:                      // *
 625:                      // *           Want left and right eigenvectors
 626:                      // *           Copy Schur vectors to VR
 627:                      // *
 628:                      FortranLib.Copy(ref SIDE , "B");
 629:                      this._dlacpy.Run("F", N, N, VL, offset_vl, LDVL, ref VR, offset_vr
 630:                                       , LDVR);
 631:                  }
 632:                  // *
 633:              }
 634:              else
 635:              {
 636:                  if (WANTVR)
 637:                  {
 638:                      // *
 639:                      // *        Want right eigenvectors
 640:                      // *        Copy Householder vectors to VR
 641:                      // *
 642:                      FortranLib.Copy(ref SIDE , "R");
 643:                      this._dlacpy.Run("L", N, N, A, offset_a, LDA, ref VR, offset_vr
 644:                                       , LDVR);
 645:                      // *
 646:                      // *        Generate orthogonal matrix in VR
 647:                      // *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
 648:                      // *
 649:                      this._dorghr.Run(N, ILO, IHI, ref VR, offset_vr, LDVR, WORK, ITAU + o_work
 650:                                       , ref WORK, IWRK + o_work, LWORK - IWRK + 1, ref IERR);
 651:                      // *
 652:                      // *        Perform QR iteration, accumulating Schur vectors in VR
 653:                      // *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
 654:                      // *
 655:                      IWRK = ITAU;
 656:                      this._dhseqr.Run("S", "V", N, ILO, IHI, ref A, offset_a
 657:                                       , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VR, offset_vr, LDVR, ref WORK, IWRK + o_work
 658:                                       , LWORK - IWRK + 1, ref INFO);
 659:                      // *
 660:                  }
 661:                  else
 662:                  {
 663:                      // *
 664:                      // *        Compute eigenvalues only
 665:                      // *        (Workspace: need N+1, prefer N+HSWORK (see comments) )
 666:                      // *
 667:                      IWRK = ITAU;
 668:                      this._dhseqr.Run("E", "N", N, ILO, IHI, ref A, offset_a
 669:                                       , LDA, ref WR, offset_wr, ref WI, offset_wi, ref VR, offset_vr, LDVR, ref WORK, IWRK + o_work
 670:                                       , LWORK - IWRK + 1, ref INFO);
 671:                  }
 672:              }
 673:              // *
 674:              // *     If INFO > 0 from DHSEQR, then quit
 675:              // *
 676:              if (INFO > 0) goto LABEL50;
 677:              // *
 678:              if (WANTVL || WANTVR)
 679:              {
 680:                  // *
 681:                  // *        Compute left and/or right eigenvectors
 682:                  // *        (Workspace: need 4*N)
 683:                  // *
 684:                  this._dtrevc.Run(SIDE, "B", ref SELECT, offset_select, N, A, offset_a, LDA
 685:                                   , ref VL, offset_vl, LDVL, ref VR, offset_vr, LDVR, N, ref NOUT
 686:                                   , ref WORK, IWRK + o_work, ref IERR);
 687:              }
 688:              // *
 689:              if (WANTVL)
 690:              {
 691:                  // *
 692:                  // *        Undo balancing of left eigenvectors
 693:                  // *        (Workspace: need N)
 694:                  // *
 695:                  this._dgebak.Run("B", "L", N, ILO, IHI, WORK, IBAL + o_work
 696:                                   , N, ref VL, offset_vl, LDVL, ref IERR);
 697:                  // *
 698:                  // *        Normalize left eigenvectors and make largest component real
 699:                  // *
 700:                  for (I = 1; I <= N; I++)
 701:                  {
 702:                      if (WI[I + o_wi] == ZERO)
 703:                      {
 704:                          SCL = ONE / this._dnrm2.Run(N, VL, 1+I * LDVL + o_vl, 1);
 705:                          this._dscal.Run(N, SCL, ref VL, 1+I * LDVL + o_vl, 1);
 706:                      }
 707:                      else
 708:                      {
 709:                          if (WI[I + o_wi] > ZERO)
 710:                          {
 711:                              SCL = ONE / this._dlapy2.Run(this._dnrm2.Run(N, VL, 1+I * LDVL + o_vl, 1), this._dnrm2.Run(N, VL, 1+(I + 1) * LDVL + o_vl, 1));
 712:                              this._dscal.Run(N, SCL, ref VL, 1+I * LDVL + o_vl, 1);
 713:                              this._dscal.Run(N, SCL, ref VL, 1+(I + 1) * LDVL + o_vl, 1);
 714:                              for (K = 1; K <= N; K++)
 715:                              {
 716:                                  WORK[IWRK + K - 1 + o_work] = Math.Pow(VL[K+I * LDVL + o_vl],2) + Math.Pow(VL[K+(I + 1) * LDVL + o_vl],2);
 717:                              }
 718:                              K = this._idamax.Run(N, WORK, IWRK + o_work, 1);
 719:                              this._dlartg.Run(VL[K+I * LDVL + o_vl], VL[K+(I + 1) * LDVL + o_vl], ref CS, ref SN, ref R);
 720:                              this._drot.Run(N, ref VL, 1+I * LDVL + o_vl, 1, ref VL, 1+(I + 1) * LDVL + o_vl, 1, CS
 721:                                             , SN);
 722:                              VL[K+(I + 1) * LDVL + o_vl] = ZERO;
 723:                          }
 724:                      }
 725:                  }
 726:              }
 727:              // *
 728:              if (WANTVR)
 729:              {
 730:                  // *
 731:                  // *        Undo balancing of right eigenvectors
 732:                  // *        (Workspace: need N)
 733:                  // *
 734:                  this._dgebak.Run("B", "R", N, ILO, IHI, WORK, IBAL + o_work
 735:                                   , N, ref VR, offset_vr, LDVR, ref IERR);
 736:                  // *
 737:                  // *        Normalize right eigenvectors and make largest component real
 738:                  // *
 739:                  for (I = 1; I <= N; I++)
 740:                  {
 741:                      if (WI[I + o_wi] == ZERO)
 742:                      {
 743:                          SCL = ONE / this._dnrm2.Run(N, VR, 1+I * LDVR + o_vr, 1);
 744:                          this._dscal.Run(N, SCL, ref VR, 1+I * LDVR + o_vr, 1);
 745:                      }
 746:                      else
 747:                      {
 748:                          if (WI[I + o_wi] > ZERO)
 749:                          {
 750:                              SCL = ONE / this._dlapy2.Run(this._dnrm2.Run(N, VR, 1+I * LDVR + o_vr, 1), this._dnrm2.Run(N, VR, 1+(I + 1) * LDVR + o_vr, 1));
 751:                              this._dscal.Run(N, SCL, ref VR, 1+I * LDVR + o_vr, 1);
 752:                              this._dscal.Run(N, SCL, ref VR, 1+(I + 1) * LDVR + o_vr, 1);
 753:                              for (K = 1; K <= N; K++)
 754:                              {
 755:                                  WORK[IWRK + K - 1 + o_work] = Math.Pow(VR[K+I * LDVR + o_vr],2) + Math.Pow(VR[K+(I + 1) * LDVR + o_vr],2);
 756:                              }
 757:                              K = this._idamax.Run(N, WORK, IWRK + o_work, 1);
 758:                              this._dlartg.Run(VR[K+I * LDVR + o_vr], VR[K+(I + 1) * LDVR + o_vr], ref CS, ref SN, ref R);
 759:                              this._drot.Run(N, ref VR, 1+I * LDVR + o_vr, 1, ref VR, 1+(I + 1) * LDVR + o_vr, 1, CS
 760:                                             , SN);
 761:                              VR[K+(I + 1) * LDVR + o_vr] = ZERO;
 762:                          }
 763:                      }
 764:                  }
 765:              }
 766:              // *
 767:              // *     Undo scaling if necessary
 768:              // *
 769:          LABEL50:;
 770:              if (SCALEA)
 771:              {
 772:                  this._dlascl.Run("G", 0, 0, CSCALE, ANRM, N - INFO
 773:                                   , 1, ref WR, INFO + 1 + o_wr, Math.Max(N - INFO, 1), ref IERR);
 774:                  this._dlascl.Run("G", 0, 0, CSCALE, ANRM, N - INFO
 775:                                   , 1, ref WI, INFO + 1 + o_wi, Math.Max(N - INFO, 1), ref IERR);
 776:                  if (INFO > 0)
 777:                  {
 778:                      this._dlascl.Run("G", 0, 0, CSCALE, ANRM, ILO - 1
 779:                                       , 1, ref WR, offset_wr, N, ref IERR);
 780:                      this._dlascl.Run("G", 0, 0, CSCALE, ANRM, ILO - 1
 781:                                       , 1, ref WI, offset_wi, N, ref IERR);
 782:                  }
 783:              }
 784:              // *
 785:              WORK[1 + o_work] = MAXWRK;
 786:              return;
 787:              // *
 788:              // *     End of DGEEV
 789:              // *
 790:   
 791:              #endregion
 792:   
 793:          }
 794:      }
 795:  }