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 auxiliary routine (version 3.1) --
  21:      /// Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
  22:      /// November 2006
  23:      /// Purpose
  24:      /// =======
  25:      /// 
  26:      /// DLAQR0 computes the eigenvalues of a Hessenberg matrix H
  27:      /// and, optionally, the matrices T and Z from the Schur decomposition
  28:      /// H = Z T Z**T, where T is an upper quasi-triangular matrix (the
  29:      /// Schur form), and Z is the orthogonal matrix of Schur vectors.
  30:      /// 
  31:      /// Optionally Z may be postmultiplied into an input orthogonal
  32:      /// matrix Q so that this routine can give the Schur factorization
  33:      /// of a matrix A which has been reduced to the Hessenberg form H
  34:      /// by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
  35:      /// 
  36:      ///</summary>
  37:      public class DLAQR0
  38:      {
  39:      
  40:   
  41:          #region Dependencies
  42:          
  43:          ILAENV _ilaenv; DLACPY _dlacpy; DLAHQR _dlahqr; DLANV2 _dlanv2; DLAQR3 _dlaqr3; DLAQR4 _dlaqr4; DLAQR5 _dlaqr5; 
  44:   
  45:          #endregion
  46:   
  47:   
  48:          #region Fields
  49:          
  50:          const int NTINY = 11; const int KEXNW = 5; const int KEXSH = 6; const double WILK1 = 0.75E0; 
  51:          const double WILK2 =  - 0.4375E0;const double ZERO = 0.0E0; const double ONE = 1.0E0; double AA = 0; double BB = 0; 
  52:          double CC = 0;double CS = 0; double DD = 0; double SN = 0; double SS = 0; double SWAP = 0; int I = 0; int INF = 0; 
  53:          int IT = 0;int ITMAX = 0; int K = 0; int KACC22 = 0; int KBOT = 0; int KDU = 0; int KS = 0; int KT = 0; int KTOP = 0; 
  54:          int KU = 0;int KV = 0; int KWH = 0; int KWTOP = 0; int KWV = 0; int LD = 0; int LS = 0; int LWKOPT = 0; int NDFL = 0; 
  55:          int NH = 0;int NHO = 0; int NIBBLE = 0; int NMIN = 0; int NS = 0; int NSMAX = 0; int NSR = 0; int NVE = 0; int NW = 0; 
  56:          int NWMAX = 0;int NWR = 0; bool NWINC = false; bool SORTED = false; string JBCMPZ = new string(' ', 2); 
  57:          double[] ZDUM = new double[1 * 1]; int offset_zdum = 0;
  58:   
  59:          #endregion
  60:   
  61:          public DLAQR0(ILAENV ilaenv, DLACPY dlacpy, DLAHQR dlahqr, DLANV2 dlanv2, DLAQR3 dlaqr3, DLAQR4 dlaqr4, DLAQR5 dlaqr5)
  62:          {
  63:      
  64:   
  65:              #region Set Dependencies
  66:              
  67:              this._ilaenv = ilaenv; this._dlacpy = dlacpy; this._dlahqr = dlahqr; this._dlanv2 = dlanv2; this._dlaqr3 = dlaqr3; 
  68:              this._dlaqr4 = dlaqr4;this._dlaqr5 = dlaqr5; 
  69:   
  70:              #endregion
  71:   
  72:          }
  73:      
  74:          public DLAQR0()
  75:          {
  76:      
  77:   
  78:              #region Dependencies (Initialization)
  79:              
  80:              IEEECK ieeeck = new IEEECK();
  81:              IPARMQ iparmq = new IPARMQ();
  82:              LSAME lsame = new LSAME();
  83:              DLAMC3 dlamc3 = new DLAMC3();
  84:              DCOPY dcopy = new DCOPY();
  85:              DLABAD dlabad = new DLABAD();
  86:              DLAPY2 dlapy2 = new DLAPY2();
  87:              DNRM2 dnrm2 = new DNRM2();
  88:              DSCAL dscal = new DSCAL();
  89:              DROT drot = new DROT();
  90:              DAXPY daxpy = new DAXPY();
  91:              XERBLA xerbla = new XERBLA();
  92:              DLASSQ dlassq = new DLASSQ();
  93:              IDAMAX idamax = new IDAMAX();
  94:              DSWAP dswap = new DSWAP();
  95:              DLAQR1 dlaqr1 = new DLAQR1();
  96:              ILAENV ilaenv = new ILAENV(ieeeck, iparmq);
  97:              DLACPY dlacpy = new DLACPY(lsame);
  98:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  99:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
 100:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
 101:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
 102:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
 103:              DLANV2 dlanv2 = new DLANV2(dlamch, dlapy2);
 104:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
 105:              DLAHQR dlahqr = new DLAHQR(dlamch, dcopy, dlabad, dlanv2, dlarfg, drot);
 106:              DGEMV dgemv = new DGEMV(lsame, xerbla);
 107:              DGER dger = new DGER(xerbla);
 108:              DLARF dlarf = new DLARF(dgemv, dger, lsame);
 109:              DGEHD2 dgehd2 = new DGEHD2(dlarf, dlarfg, xerbla);
 110:              DGEMM dgemm = new DGEMM(lsame, xerbla);
 111:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
 112:              DTRMV dtrmv = new DTRMV(lsame, xerbla);
 113:              DLAHR2 dlahr2 = new DLAHR2(daxpy, dcopy, dgemm, dgemv, dlacpy, dlarfg, dscal, dtrmm, dtrmv);
 114:              DLARFB dlarfb = new DLARFB(lsame, dcopy, dgemm, dtrmm);
 115:              DGEHRD dgehrd = new DGEHRD(daxpy, dgehd2, dgemm, dlahr2, dlarfb, dtrmm, xerbla, ilaenv);
 116:              DLASET dlaset = new DLASET(lsame);
 117:              DLARFT dlarft = new DLARFT(dgemv, dtrmv, lsame);
 118:              DORG2R dorg2r = new DORG2R(dlarf, dscal, xerbla);
 119:              DORGQR dorgqr = new DORGQR(dlarfb, dlarft, dorg2r, xerbla, ilaenv);
 120:              DORGHR dorghr = new DORGHR(dorgqr, xerbla, ilaenv);
 121:              DLANGE dlange = new DLANGE(dlassq, lsame);
 122:              DLARFX dlarfx = new DLARFX(lsame, dgemv, dger);
 123:              DLARTG dlartg = new DLARTG(dlamch);
 124:              DLASY2 dlasy2 = new DLASY2(idamax, dlamch, dcopy, dswap);
 125:              DLAEXC dlaexc = new DLAEXC(dlamch, dlange, dlacpy, dlanv2, dlarfg, dlarfx, dlartg, dlasy2, drot);
 126:              DTREXC dtrexc = new DTREXC(lsame, dlaexc, xerbla);
 127:              DLAQR2 dlaqr2 = new DLAQR2(dlamch, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlarf, dlarfg
 128:                                         , dlaset, dorghr, dtrexc);
 129:              DLAQR5 dlaqr5 = new DLAQR5(dlamch, dgemm, dlabad, dlacpy, dlaqr1, dlarfg, dlaset, dtrmm);
 130:              DLAQR4 dlaqr4 = new DLAQR4(ilaenv, dlacpy, dlahqr, dlanv2, dlaqr2, dlaqr5);
 131:              DLAQR3 dlaqr3 = new DLAQR3(dlamch, ilaenv, dcopy, dgehrd, dgemm, dlabad, dlacpy, dlahqr, dlanv2, dlaqr4
 132:                                         , dlarf, dlarfg, dlaset, dorghr, dtrexc);
 133:   
 134:              #endregion
 135:   
 136:   
 137:              #region Set Dependencies
 138:              
 139:              this._ilaenv = ilaenv; this._dlacpy = dlacpy; this._dlahqr = dlahqr; this._dlanv2 = dlanv2; this._dlaqr3 = dlaqr3; 
 140:              this._dlaqr4 = dlaqr4;this._dlaqr5 = dlaqr5; 
 141:   
 142:              #endregion
 143:   
 144:          }
 145:          /// <summary>
 146:          /// Purpose
 147:          /// =======
 148:          /// 
 149:          /// DLAQR0 computes the eigenvalues of a Hessenberg matrix H
 150:          /// and, optionally, the matrices T and Z from the Schur decomposition
 151:          /// H = Z T Z**T, where T is an upper quasi-triangular matrix (the
 152:          /// Schur form), and Z is the orthogonal matrix of Schur vectors.
 153:          /// 
 154:          /// Optionally Z may be postmultiplied into an input orthogonal
 155:          /// matrix Q so that this routine can give the Schur factorization
 156:          /// of a matrix A which has been reduced to the Hessenberg form H
 157:          /// by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
 158:          /// 
 159:          ///</summary>
 160:          /// <param name="WANTT">
 161:          /// (input) LOGICAL
 162:          /// = .TRUE. : the full Schur form T is required;
 163:          /// = .FALSE.: only eigenvalues are required.
 164:          ///</param>
 165:          /// <param name="WANTZ">
 166:          /// (input) LOGICAL
 167:          /// = .TRUE. : the matrix of Schur vectors Z is required;
 168:          /// = .FALSE.: Schur vectors are not required.
 169:          ///</param>
 170:          /// <param name="N">
 171:          /// (input) INTEGER
 172:          /// The order of the matrix H.  N .GE. 0.
 173:          ///</param>
 174:          /// <param name="ILO">
 175:          /// (input) INTEGER
 176:          ///</param>
 177:          /// <param name="IHI">
 178:          /// (input) INTEGER
 179:          /// It is assumed that H is already upper triangular in rows
 180:          /// and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
 181:          /// H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
 182:          /// previous call to DGEBAL, and then passed to DGEHRD when the
 183:          /// matrix output by DGEBAL is reduced to Hessenberg form.
 184:          /// Otherwise, ILO and IHI should be set to 1 and N,
 185:          /// respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
 186:          /// If N = 0, then ILO = 1 and IHI = 0.
 187:          ///</param>
 188:          /// <param name="H">
 189:          /// (input/output) DOUBLE PRECISION array, dimension (LDH,N)
 190:          /// On entry, the upper Hessenberg matrix H.
 191:          /// On exit, if INFO = 0 and WANTT is .TRUE., then H contains
 192:          /// the upper quasi-triangular matrix T from the Schur
 193:          /// decomposition (the Schur form); 2-by-2 diagonal blocks
 194:          /// (corresponding to complex conjugate pairs of eigenvalues)
 195:          /// are returned in standard form, with H(i,i) = H(i+1,i+1)
 196:          /// and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is
 197:          /// .FALSE., then the contents of H are unspecified on exit.
 198:          /// (The output value of H when INFO.GT.0 is given under the
 199:          /// description of INFO below.)
 200:          /// 
 201:          /// This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
 202:          /// j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
 203:          ///</param>
 204:          /// <param name="LDH">
 205:          /// (input) INTEGER
 206:          /// The leading dimension of the array H. LDH .GE. max(1,N).
 207:          ///</param>
 208:          /// <param name="WR">
 209:          /// (output) DOUBLE PRECISION array, dimension (IHI)
 210:          ///</param>
 211:          /// <param name="WI">
 212:          /// (output) DOUBLE PRECISION array, dimension (IHI)
 213:          /// The real and imaginary parts, respectively, of the computed
 214:          /// eigenvalues of H(ILO:IHI,ILO:IHI) are stored WR(ILO:IHI)
 215:          /// and WI(ILO:IHI). If two eigenvalues are computed as a
 216:          /// complex conjugate pair, they are stored in consecutive
 217:          /// elements of WR and WI, say the i-th and (i+1)th, with
 218:          /// WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then
 219:          /// the eigenvalues are stored in the same order as on the
 220:          /// diagonal of the Schur form returned in H, with
 221:          /// WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
 222:          /// block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
 223:          /// WI(i+1) = -WI(i).
 224:          ///</param>
 225:          /// <param name="ILOZ">
 226:          /// (input) INTEGER
 227:          ///</param>
 228:          /// <param name="IHIZ">
 229:          /// (input) INTEGER
 230:          /// Specify the rows of Z to which transformations must be
 231:          /// applied if WANTZ is .TRUE..
 232:          /// 1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
 233:          ///</param>
 234:          /// <param name="Z">
 235:          /// (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI)
 236:          /// If WANTZ is .FALSE., then Z is not referenced.
 237:          /// If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
 238:          /// replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
 239:          /// orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
 240:          /// (The output value of Z when INFO.GT.0 is given under
 241:          /// the description of INFO below.)
 242:          ///</param>
 243:          /// <param name="LDZ">
 244:          /// (input) INTEGER
 245:          /// The leading dimension of the array Z.  if WANTZ is .TRUE.
 246:          /// then LDZ.GE.MAX(1,IHIZ).  Otherwize, LDZ.GE.1.
 247:          ///</param>
 248:          /// <param name="WORK">
 249:          /// (workspace/output) DOUBLE PRECISION array, dimension LWORK
 250:          /// On exit, if LWORK = -1, WORK(1) returns an estimate of
 251:          /// the optimal value for LWORK.
 252:          ///</param>
 253:          /// <param name="LWORK">
 254:          /// (input) INTEGER
 255:          /// The dimension of the array WORK.  LWORK .GE. max(1,N)
 256:          /// is sufficient, but LWORK typically as large as 6*N may
 257:          /// be required for optimal performance.  A workspace query
 258:          /// to determine the optimal workspace size is recommended.
 259:          /// 
 260:          /// If LWORK = -1, then DLAQR0 does a workspace query.
 261:          /// In this case, DLAQR0 checks the input parameters and
 262:          /// estimates the optimal workspace size for the given
 263:          /// values of N, ILO and IHI.  The estimate is returned
 264:          /// in WORK(1).  No error message related to LWORK is
 265:          /// issued by XERBLA.  Neither H nor Z are accessed.
 266:          /// 
 267:          ///</param>
 268:          /// <param name="INFO">
 269:          /// (output) INTEGER
 270:          /// =  0:  successful exit
 271:          /// .GT. 0:  if INFO = i, DLAQR0 failed to compute all of
 272:          /// the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
 273:          /// and WI contain those eigenvalues which have been
 274:          /// successfully computed.  (Failures are rare.)
 275:          /// 
 276:          /// If INFO .GT. 0 and WANT is .FALSE., then on exit,
 277:          /// the remaining unconverged eigenvalues are the eigen-
 278:          /// values of the upper Hessenberg matrix rows and
 279:          /// columns ILO through INFO of the final, output
 280:          /// value of H.
 281:          /// 
 282:          /// If INFO .GT. 0 and WANTT is .TRUE., then on exit
 283:          /// 
 284:          /// (*)  (initial value of H)*U  = U*(final value of H)
 285:          /// 
 286:          /// where U is an orthogonal matrix.  The final
 287:          /// value of H is upper Hessenberg and quasi-triangular
 288:          /// in rows and columns INFO+1 through IHI.
 289:          /// 
 290:          /// If INFO .GT. 0 and WANTZ is .TRUE., then on exit
 291:          /// 
 292:          /// (final value of Z(ILO:IHI,ILOZ:IHIZ)
 293:          /// =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
 294:          /// 
 295:          /// where U is the orthogonal matrix in (*) (regard-
 296:          /// less of the value of WANTT.)
 297:          /// 
 298:          /// If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
 299:          /// accessed.
 300:          /// 
 301:          ///</param>
 302:          public void Run(bool WANTT, bool WANTZ, int N, int ILO, int IHI, ref double[] H, int offset_h
 303:                           , int LDH, ref double[] WR, int offset_wr, ref double[] WI, int offset_wi, int ILOZ, int IHIZ, ref double[] Z, int offset_z
 304:                           , int LDZ, ref double[] WORK, int offset_work, int LWORK, ref int INFO)
 305:          {
 306:   
 307:              #region Array Index Correction
 308:              
 309:               int o_h = -1 - LDH + offset_h;  int o_wr = -1 + offset_wr;  int o_wi = -1 + offset_wi; 
 310:               int o_z = -1 - LDZ + offset_z; int o_work = -1 + offset_work; 
 311:   
 312:              #endregion
 313:   
 314:   
 315:              #region Prolog
 316:              
 317:              // *
 318:              // *  -- LAPACK auxiliary routine (version 3.1) --
 319:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 320:              // *     November 2006
 321:              // *
 322:              // *     .. Scalar Arguments ..
 323:              // *     ..
 324:              // *     .. Array Arguments ..
 325:              // *     ..
 326:              // *
 327:              // *     Purpose
 328:              // *     =======
 329:              // *
 330:              // *     DLAQR0 computes the eigenvalues of a Hessenberg matrix H
 331:              // *     and, optionally, the matrices T and Z from the Schur decomposition
 332:              // *     H = Z T Z**T, where T is an upper quasi-triangular matrix (the
 333:              // *     Schur form), and Z is the orthogonal matrix of Schur vectors.
 334:              // *
 335:              // *     Optionally Z may be postmultiplied into an input orthogonal
 336:              // *     matrix Q so that this routine can give the Schur factorization
 337:              // *     of a matrix A which has been reduced to the Hessenberg form H
 338:              // *     by the orthogonal matrix Q:  A = Q*H*Q**T = (QZ)*T*(QZ)**T.
 339:              // *
 340:              // *     Arguments
 341:              // *     =========
 342:              // *
 343:              // *     WANTT   (input) LOGICAL
 344:              // *          = .TRUE. : the full Schur form T is required;
 345:              // *          = .FALSE.: only eigenvalues are required.
 346:              // *
 347:              // *     WANTZ   (input) LOGICAL
 348:              // *          = .TRUE. : the matrix of Schur vectors Z is required;
 349:              // *          = .FALSE.: Schur vectors are not required.
 350:              // *
 351:              // *     N     (input) INTEGER
 352:              // *           The order of the matrix H.  N .GE. 0.
 353:              // *
 354:              // *     ILO   (input) INTEGER
 355:              // *     IHI   (input) INTEGER
 356:              // *           It is assumed that H is already upper triangular in rows
 357:              // *           and columns 1:ILO-1 and IHI+1:N and, if ILO.GT.1,
 358:              // *           H(ILO,ILO-1) is zero. ILO and IHI are normally set by a
 359:              // *           previous call to DGEBAL, and then passed to DGEHRD when the
 360:              // *           matrix output by DGEBAL is reduced to Hessenberg form.
 361:              // *           Otherwise, ILO and IHI should be set to 1 and N,
 362:              // *           respectively.  If N.GT.0, then 1.LE.ILO.LE.IHI.LE.N.
 363:              // *           If N = 0, then ILO = 1 and IHI = 0.
 364:              // *
 365:              // *     H     (input/output) DOUBLE PRECISION array, dimension (LDH,N)
 366:              // *           On entry, the upper Hessenberg matrix H.
 367:              // *           On exit, if INFO = 0 and WANTT is .TRUE., then H contains
 368:              // *           the upper quasi-triangular matrix T from the Schur
 369:              // *           decomposition (the Schur form); 2-by-2 diagonal blocks
 370:              // *           (corresponding to complex conjugate pairs of eigenvalues)
 371:              // *           are returned in standard form, with H(i,i) = H(i+1,i+1)
 372:              // *           and H(i+1,i)*H(i,i+1).LT.0. If INFO = 0 and WANTT is
 373:              // *           .FALSE., then the contents of H are unspecified on exit.
 374:              // *           (The output value of H when INFO.GT.0 is given under the
 375:              // *           description of INFO below.)
 376:              // *
 377:              // *           This subroutine may explicitly set H(i,j) = 0 for i.GT.j and
 378:              // *           j = 1, 2, ... ILO-1 or j = IHI+1, IHI+2, ... N.
 379:              // *
 380:              // *     LDH   (input) INTEGER
 381:              // *           The leading dimension of the array H. LDH .GE. max(1,N).
 382:              // *
 383:              // *     WR    (output) DOUBLE PRECISION array, dimension (IHI)
 384:              // *     WI    (output) DOUBLE PRECISION array, dimension (IHI)
 385:              // *           The real and imaginary parts, respectively, of the computed
 386:              // *           eigenvalues of H(ILO:IHI,ILO:IHI) are stored WR(ILO:IHI)
 387:              // *           and WI(ILO:IHI). If two eigenvalues are computed as a
 388:              // *           complex conjugate pair, they are stored in consecutive
 389:              // *           elements of WR and WI, say the i-th and (i+1)th, with
 390:              // *           WI(i) .GT. 0 and WI(i+1) .LT. 0. If WANTT is .TRUE., then
 391:              // *           the eigenvalues are stored in the same order as on the
 392:              // *           diagonal of the Schur form returned in H, with
 393:              // *           WR(i) = H(i,i) and, if H(i:i+1,i:i+1) is a 2-by-2 diagonal
 394:              // *           block, WI(i) = sqrt(-H(i+1,i)*H(i,i+1)) and
 395:              // *           WI(i+1) = -WI(i).
 396:              // *
 397:              // *     ILOZ     (input) INTEGER
 398:              // *     IHIZ     (input) INTEGER
 399:              // *           Specify the rows of Z to which transformations must be
 400:              // *           applied if WANTZ is .TRUE..
 401:              // *           1 .LE. ILOZ .LE. ILO; IHI .LE. IHIZ .LE. N.
 402:              // *
 403:              // *     Z     (input/output) DOUBLE PRECISION array, dimension (LDZ,IHI)
 404:              // *           If WANTZ is .FALSE., then Z is not referenced.
 405:              // *           If WANTZ is .TRUE., then Z(ILO:IHI,ILOZ:IHIZ) is
 406:              // *           replaced by Z(ILO:IHI,ILOZ:IHIZ)*U where U is the
 407:              // *           orthogonal Schur factor of H(ILO:IHI,ILO:IHI).
 408:              // *           (The output value of Z when INFO.GT.0 is given under
 409:              // *           the description of INFO below.)
 410:              // *
 411:              // *     LDZ   (input) INTEGER
 412:              // *           The leading dimension of the array Z.  if WANTZ is .TRUE.
 413:              // *           then LDZ.GE.MAX(1,IHIZ).  Otherwize, LDZ.GE.1.
 414:              // *
 415:              // *     WORK  (workspace/output) DOUBLE PRECISION array, dimension LWORK
 416:              // *           On exit, if LWORK = -1, WORK(1) returns an estimate of
 417:              // *           the optimal value for LWORK.
 418:              // *
 419:              // *     LWORK (input) INTEGER
 420:              // *           The dimension of the array WORK.  LWORK .GE. max(1,N)
 421:              // *           is sufficient, but LWORK typically as large as 6*N may
 422:              // *           be required for optimal performance.  A workspace query
 423:              // *           to determine the optimal workspace size is recommended.
 424:              // *
 425:              // *           If LWORK = -1, then DLAQR0 does a workspace query.
 426:              // *           In this case, DLAQR0 checks the input parameters and
 427:              // *           estimates the optimal workspace size for the given
 428:              // *           values of N, ILO and IHI.  The estimate is returned
 429:              // *           in WORK(1).  No error message related to LWORK is
 430:              // *           issued by XERBLA.  Neither H nor Z are accessed.
 431:              // *
 432:              // *
 433:              // *     INFO  (output) INTEGER
 434:              // *             =  0:  successful exit
 435:              // *           .GT. 0:  if INFO = i, DLAQR0 failed to compute all of
 436:              // *                the eigenvalues.  Elements 1:ilo-1 and i+1:n of WR
 437:              // *                and WI contain those eigenvalues which have been
 438:              // *                successfully computed.  (Failures are rare.)
 439:              // *
 440:              // *                If INFO .GT. 0 and WANT is .FALSE., then on exit,
 441:              // *                the remaining unconverged eigenvalues are the eigen-
 442:              // *                values of the upper Hessenberg matrix rows and
 443:              // *                columns ILO through INFO of the final, output
 444:              // *                value of H.
 445:              // *
 446:              // *                If INFO .GT. 0 and WANTT is .TRUE., then on exit
 447:              // *
 448:              // *           (*)  (initial value of H)*U  = U*(final value of H)
 449:              // *
 450:              // *                where U is an orthogonal matrix.  The final
 451:              // *                value of H is upper Hessenberg and quasi-triangular
 452:              // *                in rows and columns INFO+1 through IHI.
 453:              // *
 454:              // *                If INFO .GT. 0 and WANTZ is .TRUE., then on exit
 455:              // *
 456:              // *                  (final value of Z(ILO:IHI,ILOZ:IHIZ)
 457:              // *                   =  (initial value of Z(ILO:IHI,ILOZ:IHIZ)*U
 458:              // *
 459:              // *                where U is the orthogonal matrix in (*) (regard-
 460:              // *                less of the value of WANTT.)
 461:              // *
 462:              // *                If INFO .GT. 0 and WANTZ is .FALSE., then Z is not
 463:              // *                accessed.
 464:              // *
 465:              // *
 466:              // *     ================================================================
 467:              // *     Based on contributions by
 468:              // *        Karen Braman and Ralph Byers, Department of Mathematics,
 469:              // *        University of Kansas, USA
 470:              // *
 471:              // *     ================================================================
 472:              // *
 473:              // *     References:
 474:              // *       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
 475:              // *       Algorithm Part I: Maintaining Well Focused Shifts, and Level 3
 476:              // *       Performance, SIAM Journal of Matrix Analysis, volume 23, pages
 477:              // *       929--947, 2002.
 478:              // *
 479:              // *       K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
 480:              // *       Algorithm Part II: Aggressive Early Deflation, SIAM Journal
 481:              // *       of Matrix Analysis, volume 23, pages 948--973, 2002.
 482:              // *
 483:              // *     ================================================================
 484:              // *     .. Parameters ..
 485:              // *
 486:              // *     ==== Matrices of order NTINY or smaller must be processed by
 487:              // *     .    DLAHQR because of insufficient subdiagonal scratch space.
 488:              // *     .    (This is a hard limit.) ====
 489:              // *
 490:              // *     ==== Exceptional deflation windows:  try to cure rare
 491:              // *     .    slow convergence by increasing the size of the
 492:              // *     .    deflation window after KEXNW iterations. =====
 493:              // *
 494:              // *     ==== Exceptional shifts: try to cure rare slow convergence
 495:              // *     .    with ad-hoc exceptional shifts every KEXSH iterations.
 496:              // *     .    The constants WILK1 and WILK2 are used to form the
 497:              // *     .    exceptional shifts. ====
 498:              // *
 499:              // *     ..
 500:              // *     .. Local Scalars ..
 501:              // *     ..
 502:              // *     .. External Functions ..
 503:              // *     ..
 504:              // *     .. Local Arrays ..
 505:              // *     ..
 506:              // *     .. External Subroutines ..
 507:              // *     ..
 508:              // *     .. Intrinsic Functions ..
 509:              //      INTRINSIC          ABS, DBLE, INT, MAX, MIN, MOD;
 510:              // *     ..
 511:              // *     .. Executable Statements ..
 512:   
 513:              #endregion
 514:   
 515:   
 516:              #region Body
 517:              
 518:              INFO = 0;
 519:              // *
 520:              // *     ==== Quick return for N = 0: nothing to do. ====
 521:              // *
 522:              if (N == 0)
 523:              {
 524:                  WORK[1 + o_work] = ONE;
 525:                  return;
 526:              }
 527:              // *
 528:              // *     ==== Set up job flags for ILAENV. ====
 529:              // *
 530:              if (WANTT)
 531:              {
 532:                  FortranLib.Copy(ref JBCMPZ, 1, 1, "S");
 533:              }
 534:              else
 535:              {
 536:                  FortranLib.Copy(ref JBCMPZ, 1, 1, "E");
 537:              }
 538:              if (WANTZ)
 539:              {
 540:                  FortranLib.Copy(ref JBCMPZ, 2, 2, "V");
 541:              }
 542:              else
 543:              {
 544:                  FortranLib.Copy(ref JBCMPZ, 2, 2, "N");
 545:              }
 546:              // *
 547:              // *     ==== Tiny matrices must use DLAHQR. ====
 548:              // *
 549:              if (N <= NTINY)
 550:              {
 551:                  // *
 552:                  // *        ==== Estimate optimal workspace. ====
 553:                  // *
 554:                  LWKOPT = 1;
 555:                  if (LWORK !=  - 1)
 556:                  {
 557:                      this._dlahqr.Run(WANTT, WANTZ, N, ILO, IHI, ref H, offset_h
 558:                                       , LDH, ref WR, offset_wr, ref WI, offset_wi, ILOZ, IHIZ, ref Z, offset_z
 559:                                       , LDZ, ref INFO);
 560:                  }
 561:              }
 562:              else
 563:              {
 564:                  // *
 565:                  // *        ==== Use small bulge multi-shift QR with aggressive early
 566:                  // *        .    deflation on larger-than-tiny matrices. ====
 567:                  // *
 568:                  // *        ==== Hope for the best. ====
 569:                  // *
 570:                  INFO = 0;
 571:                  // *
 572:                  // *        ==== NWR = recommended deflation window size.  At this
 573:                  // *        .    point,  N .GT. NTINY = 11, so there is enough
 574:                  // *        .    subdiagonal workspace for NWR.GE.2 as required.
 575:                  // *        .    (In fact, there is enough subdiagonal space for
 576:                  // *        .    NWR.GE.3.) ====
 577:                  // *
 578:                  NWR = this._ilaenv.Run(13, "DLAQR0", JBCMPZ, N, ILO, IHI, LWORK);
 579:                  NWR = Math.Max(2, NWR);
 580:                  NWR = Math.Min(IHI - ILO + 1, Math.Min((N - 1) / 3, NWR));
 581:                  NW = NWR;
 582:                  // *
 583:                  // *        ==== NSR = recommended number of simultaneous shifts.
 584:                  // *        .    At this point N .GT. NTINY = 11, so there is at
 585:                  // *        .    enough subdiagonal workspace for NSR to be even
 586:                  // *        .    and greater than or equal to two as required. ====
 587:                  // *
 588:                  NSR = this._ilaenv.Run(15, "DLAQR0", JBCMPZ, N, ILO, IHI, LWORK);
 589:                  NSR = Math.Min(NSR, Math.Min((N + 6) / 9, IHI - ILO));
 590:                  NSR = Math.Max(2, NSR - FortranLib.Mod(NSR,2));
 591:                  // *
 592:                  // *        ==== Estimate optimal workspace ====
 593:                  // *
 594:                  // *        ==== Workspace query call to DLAQR3 ====
 595:                  // *
 596:                  this._dlaqr3.Run(WANTT, WANTZ, N, ILO, IHI, NWR + 1
 597:                                   , ref H, offset_h, LDH, ILOZ, IHIZ, ref Z, offset_z, LDZ
 598:                                   , ref LS, ref LD, ref WR, offset_wr, ref WI, offset_wi, ref H, offset_h, LDH
 599:                                   , N, ref H, offset_h, LDH, N, ref H, offset_h, LDH
 600:                                   , ref WORK, offset_work,  - 1);
 601:                  // *
 602:                  // *        ==== Optimal workspace = MAX(DLAQR5, DLAQR3) ====
 603:                  // *
 604:                  LWKOPT = Math.Max(3 * NSR / 2, Convert.ToInt32(Math.Truncate(WORK[1 + o_work])));
 605:                  // *
 606:                  // *        ==== Quick return in case of workspace query. ====
 607:                  // *
 608:                  if (LWORK ==  - 1)
 609:                  {
 610:                      WORK[1 + o_work] = Convert.ToDouble(LWKOPT);
 611:                      return;
 612:                  }
 613:                  // *
 614:                  // *        ==== DLAHQR/DLAQR0 crossover point ====
 615:                  // *
 616:                  NMIN = this._ilaenv.Run(12, "DLAQR0", JBCMPZ, N, ILO, IHI, LWORK);
 617:                  NMIN = Math.Max(NTINY, NMIN);
 618:                  // *
 619:                  // *        ==== Nibble crossover point ====
 620:                  // *
 621:                  NIBBLE = this._ilaenv.Run(14, "DLAQR0", JBCMPZ, N, ILO, IHI, LWORK);
 622:                  NIBBLE = Math.Max(0, NIBBLE);
 623:                  // *
 624:                  // *        ==== Accumulate reflections during ttswp?  Use block
 625:                  // *        .    2-by-2 structure during matrix-matrix multiply? ====
 626:                  // *
 627:                  KACC22 = this._ilaenv.Run(16, "DLAQR0", JBCMPZ, N, ILO, IHI, LWORK);
 628:                  KACC22 = Math.Max(0, KACC22);
 629:                  KACC22 = Math.Min(2, KACC22);
 630:                  // *
 631:                  // *        ==== NWMAX = the largest possible deflation window for
 632:                  // *        .    which there is sufficient workspace. ====
 633:                  // *
 634:                  NWMAX = Math.Min((N - 1) / 3, LWORK / 2);
 635:                  // *
 636:                  // *        ==== NSMAX = the Largest number of simultaneous shifts
 637:                  // *        .    for which there is sufficient workspace. ====
 638:                  // *
 639:                  NSMAX = Math.Min((N + 6) / 9, 2 * LWORK / 3);
 640:                  NSMAX = NSMAX - FortranLib.Mod(NSMAX,2);
 641:                  // *
 642:                  // *        ==== NDFL: an iteration count restarted at deflation. ====
 643:                  // *
 644:                  NDFL = 1;
 645:                  // *
 646:                  // *        ==== ITMAX = iteration limit ====
 647:                  // *
 648:                  ITMAX = Math.Max(30, 2 * KEXSH) * Math.Max(10, (IHI - ILO + 1));
 649:                  // *
 650:                  // *        ==== Last row and column in the active block ====
 651:                  // *
 652:                  KBOT = IHI;
 653:                  // *
 654:                  // *        ==== Main Loop ====
 655:                  // *
 656:                  for (IT = 1; IT <= ITMAX; IT++)
 657:                  {
 658:                      // *
 659:                      // *           ==== Done when KBOT falls below ILO ====
 660:                      // *
 661:                      if (KBOT < ILO) goto LABEL90;
 662:                      // *
 663:                      // *           ==== Locate active block ====
 664:                      // *
 665:                      for (K = KBOT; K >= ILO + 1; K +=  - 1)
 666:                      {
 667:                          if (H[K+(K - 1) * LDH + o_h] == ZERO) goto LABEL20;
 668:                      }
 669:                      K = ILO;
 670:                  LABEL20:;
 671:                      KTOP = K;
 672:                      // *
 673:                      // *           ==== Select deflation window size ====
 674:                      // *
 675:                      NH = KBOT - KTOP + 1;
 676:                      if (NDFL < KEXNW || NH < NW)
 677:                      {
 678:                          // *
 679:                          // *              ==== Typical deflation window.  If possible and
 680:                          // *              .    advisable, nibble the entire active block.
 681:                          // *              .    If not, use size NWR or NWR+1 depending upon
 682:                          // *              .    which has the smaller corresponding subdiagonal
 683:                          // *              .    entry (a heuristic). ====
 684:                          // *
 685:                          NWINC = true;
 686:                          if (NH <= Math.Min(NMIN, NWMAX))
 687:                          {
 688:                              NW = NH;
 689:                          }
 690:                          else
 691:                          {
 692:                              NW = Math.Min(NWR, Math.Min(NH, NWMAX));
 693:                              if (NW < NWMAX)
 694:                              {
 695:                                  if (NW >= NH - 1)
 696:                                  {
 697:                                      NW = NH;
 698:                                  }
 699:                                  else
 700:                                  {
 701:                                      KWTOP = KBOT - NW + 1;
 702:                                      if (Math.Abs(H[KWTOP+(KWTOP - 1) * LDH + o_h]) > Math.Abs(H[KWTOP - 1+(KWTOP - 2) * LDH + o_h])) NW = NW + 1;
 703:                                  }
 704:                              }
 705:                          }
 706:                      }
 707:                      else
 708:                      {
 709:                          // *
 710:                          // *              ==== Exceptional deflation window.  If there have
 711:                          // *              .    been no deflations in KEXNW or more iterations,
 712:                          // *              .    then vary the deflation window size.   At first,
 713:                          // *              .    because, larger windows are, in general, more
 714:                          // *              .    powerful than smaller ones, rapidly increase the
 715:                          // *              .    window up to the maximum reasonable and possible.
 716:                          // *              .    Then maybe try a slightly smaller window.  ====
 717:                          // *
 718:                          if (NWINC && NW < Math.Min(NWMAX, NH))
 719:                          {
 720:                              NW = Math.Min(NWMAX, Math.Min(NH, 2 * NW));
 721:                          }
 722:                          else
 723:                          {
 724:                              NWINC = false;
 725:                              if (NW == NH && NH > 2) NW = NH - 1;
 726:                          }
 727:                      }
 728:                      // *
 729:                      // *           ==== Aggressive early deflation:
 730:                      // *           .    split workspace under the subdiagonal into
 731:                      // *           .      - an nw-by-nw work array V in the lower
 732:                      // *           .        left-hand-corner,
 733:                      // *           .      - an NW-by-at-least-NW-but-more-is-better
 734:                      // *           .        (NW-by-NHO) horizontal work array along
 735:                      // *           .        the bottom edge,
 736:                      // *           .      - an at-least-NW-but-more-is-better (NHV-by-NW)
 737:                      // *           .        vertical work array along the left-hand-edge.
 738:                      // *           .        ====
 739:                      // *
 740:                      KV = N - NW + 1;
 741:                      KT = NW + 1;
 742:                      NHO = (N - NW - 1) - KT + 1;
 743:                      KWV = NW + 2;
 744:                      NVE = (N - NW) - KWV + 1;
 745:                      // *
 746:                      // *           ==== Aggressive early deflation ====
 747:                      // *
 748:                      this._dlaqr3.Run(WANTT, WANTZ, N, KTOP, KBOT, NW
 749:                                       , ref H, offset_h, LDH, ILOZ, IHIZ, ref Z, offset_z, LDZ
 750:                                       , ref LS, ref LD, ref WR, offset_wr, ref WI, offset_wi, ref H, KV+1 * LDH + o_h, LDH
 751:                                       , NHO, ref H, KV+KT * LDH + o_h, LDH, NVE, ref H, KWV+1 * LDH + o_h, LDH
 752:                                       , ref WORK, offset_work, LWORK);
 753:                      // *
 754:                      // *           ==== Adjust KBOT accounting for new deflations. ====
 755:                      // *
 756:                      KBOT = KBOT - LD;
 757:                      // *
 758:                      // *           ==== KS points to the shifts. ====
 759:                      // *
 760:                      KS = KBOT - LS + 1;
 761:                      // *
 762:                      // *           ==== Skip an expensive QR sweep if there is a (partly
 763:                      // *           .    heuristic) reason to expect that many eigenvalues
 764:                      // *           .    will deflate without it.  Here, the QR sweep is
 765:                      // *           .    skipped if many eigenvalues have just been deflated
 766:                      // *           .    or if the remaining active block is small.
 767:                      // *
 768:                      if ((LD == 0) || ((100 * LD <= NW * NIBBLE) && (KBOT - KTOP + 1 > Math.Min(NMIN, NWMAX))))
 769:                      {
 770:                          // *
 771:                          // *              ==== NS = nominal number of simultaneous shifts.
 772:                          // *              .    This may be lowered (slightly) if DLAQR3
 773:                          // *              .    did not provide that many shifts. ====
 774:                          // *
 775:                          NS = Math.Min(NSMAX, Math.Min(NSR, Math.Max(2, KBOT - KTOP)));
 776:                          NS = NS - FortranLib.Mod(NS,2);
 777:                          // *
 778:                          // *              ==== If there have been no deflations
 779:                          // *              .    in a multiple of KEXSH iterations,
 780:                          // *              .    then try exceptional shifts.
 781:                          // *              .    Otherwise use shifts provided by
 782:                          // *              .    DLAQR3 above or from the eigenvalues
 783:                          // *              .    of a trailing principal submatrix. ====
 784:                          // *
 785:                          if (FortranLib.Mod(NDFL,KEXSH) == 0)
 786:                          {
 787:                              KS = KBOT - NS + 1;
 788:                              for (I = KBOT; I >= Math.Max(KS + 1, KTOP + 2); I +=  - 2)
 789:                              {
 790:                                  SS = Math.Abs(H[I+(I - 1) * LDH + o_h]) + Math.Abs(H[I - 1+(I - 2) * LDH + o_h]);
 791:                                  AA = WILK1 * SS + H[I+I * LDH + o_h];
 792:                                  BB = SS;
 793:                                  CC = WILK2 * SS;
 794:                                  DD = AA;
 795:                                  this._dlanv2.Run(ref AA, ref BB, ref CC, ref DD, ref WR[I - 1 + o_wr], ref WI[I - 1 + o_wi]
 796:                                                   , ref WR[I + o_wr], ref WI[I + o_wi], ref CS, ref SN);
 797:                              }
 798:                              if (KS == KTOP)
 799:                              {
 800:                                  WR[KS + 1 + o_wr] = H[KS + 1+(KS + 1) * LDH + o_h];
 801:                                  WI[KS + 1 + o_wi] = ZERO;
 802:                                  WR[KS + o_wr] = WR[KS + 1 + o_wr];
 803:                                  WI[KS + o_wi] = WI[KS + 1 + o_wi];
 804:                              }
 805:                          }
 806:                          else
 807:                          {
 808:                              // *
 809:                              // *                 ==== Got NS/2 or fewer shifts? Use DLAQR4 or
 810:                              // *                 .    DLAHQR on a trailing principal submatrix to
 811:                              // *                 .    get more. (Since NS.LE.NSMAX.LE.(N+6)/9,
 812:                              // *                 .    there is enough space below the subdiagonal
 813:                              // *                 .    to fit an NS-by-NS scratch array.) ====
 814:                              // *
 815:                              if (KBOT - KS + 1 <= NS / 2)
 816:                              {
 817:                                  KS = KBOT - NS + 1;
 818:                                  KT = N - NS + 1;
 819:                                  this._dlacpy.Run("A", NS, NS, H, KS+KS * LDH + o_h, LDH, ref H, KT+1 * LDH + o_h
 820:                                                   , LDH);
 821:                                  if (NS > NMIN)
 822:                                  {
 823:                                      this._dlaqr4.Run(false, false, NS, 1, NS, ref H, KT+1 * LDH + o_h
 824:                                                       , LDH, ref WR, KS + o_wr, ref WI, KS + o_wi, 1, 1, ref ZDUM, offset_zdum
 825:                                                       , 1, ref WORK, offset_work, LWORK, ref INF);
 826:                                  }
 827:                                  else
 828:                                  {
 829:                                      this._dlahqr.Run(false, false, NS, 1, NS, ref H, KT+1 * LDH + o_h
 830:                                                       , LDH, ref WR, KS + o_wr, ref WI, KS + o_wi, 1, 1, ref ZDUM, offset_zdum
 831:                                                       , 1, ref INF);
 832:                                  }
 833:                                  KS = KS + INF;
 834:                                  // *
 835:                                  // *                    ==== In case of a rare QR failure use
 836:                                  // *                    .    eigenvalues of the trailing 2-by-2
 837:                                  // *                    .    principal submatrix.  ====
 838:                                  // *
 839:                                  if (KS >= KBOT)
 840:                                  {
 841:                                      AA = H[KBOT - 1+(KBOT - 1) * LDH + o_h];
 842:                                      CC = H[KBOT+(KBOT - 1) * LDH + o_h];
 843:                                      BB = H[KBOT - 1+KBOT * LDH + o_h];
 844:                                      DD = H[KBOT+KBOT * LDH + o_h];
 845:                                      this._dlanv2.Run(ref AA, ref BB, ref CC, ref DD, ref WR[KBOT - 1 + o_wr], ref WI[KBOT - 1 + o_wi]
 846:                                                       , ref WR[KBOT + o_wr], ref WI[KBOT + o_wi], ref CS, ref SN);
 847:                                      KS = KBOT - 1;
 848:                                  }
 849:                              }
 850:                              // *
 851:                              if (KBOT - KS + 1 > NS)
 852:                              {
 853:                                  // *
 854:                                  // *                    ==== Sort the shifts (Helps a little)
 855:                                  // *                    .    Bubble sort keeps complex conjugate
 856:                                  // *                    .    pairs together. ====
 857:                                  // *
 858:                                  SORTED = false;
 859:                                  for (K = KBOT; K >= KS + 1; K +=  - 1)
 860:                                  {
 861:                                      if (SORTED) goto LABEL60;
 862:                                      SORTED = true;
 863:                                      for (I = KS; I <= K - 1; I++)
 864:                                      {
 865:                                          if (Math.Abs(WR[I + o_wr]) + Math.Abs(WI[I + o_wi]) < Math.Abs(WR[I + 1 + o_wr]) + Math.Abs(WI[I + 1 + o_wi]))
 866:                                          {
 867:                                              SORTED = false;
 868:                                              // *
 869:                                              SWAP = WR[I + o_wr];
 870:                                              WR[I + o_wr] = WR[I + 1 + o_wr];
 871:                                              WR[I + 1 + o_wr] = SWAP;
 872:                                              // *
 873:                                              SWAP = WI[I + o_wi];
 874:                                              WI[I + o_wi] = WI[I + 1 + o_wi];
 875:                                              WI[I + 1 + o_wi] = SWAP;
 876:                                          }
 877:                                      }
 878:                                  }
 879:                              LABEL60:;
 880:                              }
 881:                              // *
 882:                              // *                 ==== Shuffle shifts into pairs of real shifts
 883:                              // *                 .    and pairs of complex conjugate shifts
 884:                              // *                 .    assuming complex conjugate shifts are
 885:                              // *                 .    already adjacent to one another. (Yes,
 886:                              // *                 .    they are.)  ====
 887:                              // *
 888:                              for (I = KBOT; I >= KS + 2; I +=  - 2)
 889:                              {
 890:                                  if (WI[I + o_wi] !=  - WI[I - 1 + o_wi])
 891:                                  {
 892:                                      // *
 893:                                      SWAP = WR[I + o_wr];
 894:                                      WR[I + o_wr] = WR[I - 1 + o_wr];
 895:                                      WR[I - 1 + o_wr] = WR[I - 2 + o_wr];
 896:                                      WR[I - 2 + o_wr] = SWAP;
 897:                                      // *
 898:                                      SWAP = WI[I + o_wi];
 899:                                      WI[I + o_wi] = WI[I - 1 + o_wi];
 900:                                      WI[I - 1 + o_wi] = WI[I - 2 + o_wi];
 901:                                      WI[I - 2 + o_wi] = SWAP;
 902:                                  }
 903:                              }
 904:                          }
 905:                          // *
 906:                          // *              ==== If there are only two shifts and both are
 907:                          // *              .    real, then use only one.  ====
 908:                          // *
 909:                          if (KBOT - KS + 1 == 2)
 910:                          {
 911:                              if (WI[KBOT + o_wi] == ZERO)
 912:                              {
 913:                                  if (Math.Abs(WR[KBOT + o_wr] - H[KBOT+KBOT * LDH + o_h]) < Math.Abs(WR[KBOT - 1 + o_wr] - H[KBOT+KBOT * LDH + o_h]))
 914:                                  {
 915:                                      WR[KBOT - 1 + o_wr] = WR[KBOT + o_wr];
 916:                                  }
 917:                                  else
 918:                                  {
 919:                                      WR[KBOT + o_wr] = WR[KBOT - 1 + o_wr];
 920:                                  }
 921:                              }
 922:                          }
 923:                          // *
 924:                          // *              ==== Use up to NS of the the smallest magnatiude
 925:                          // *              .    shifts.  If there aren't NS shifts available,
 926:                          // *              .    then use them all, possibly dropping one to
 927:                          // *              .    make the number of shifts even. ====
 928:                          // *
 929:                          NS = Math.Min(NS, KBOT - KS + 1);
 930:                          NS = NS - FortranLib.Mod(NS,2);
 931:                          KS = KBOT - NS + 1;
 932:                          // *
 933:                          // *              ==== Small-bulge multi-shift QR sweep:
 934:                          // *              .    split workspace under the subdiagonal into
 935:                          // *              .    - a KDU-by-KDU work array U in the lower
 936:                          // *              .      left-hand-corner,
 937:                          // *              .    - a KDU-by-at-least-KDU-but-more-is-better
 938:                          // *              .      (KDU-by-NHo) horizontal work array WH along
 939:                          // *              .      the bottom edge,
 940:                          // *              .    - and an at-least-KDU-but-more-is-better-by-KDU
 941:                          // *              .      (NVE-by-KDU) vertical work WV arrow along
 942:                          // *              .      the left-hand-edge. ====
 943:                          // *
 944:                          KDU = 3 * NS - 3;
 945:                          KU = N - KDU + 1;
 946:                          KWH = KDU + 1;
 947:                          NHO = (N - KDU + 1 - 4) - (KDU + 1) + 1;
 948:                          KWV = KDU + 4;
 949:                          NVE = N - KDU - KWV + 1;
 950:                          // *
 951:                          // *              ==== Small-bulge multi-shift QR sweep ====
 952:                          // *
 953:                          this._dlaqr5.Run(WANTT, WANTZ, KACC22, N, KTOP, KBOT
 954:                                           , NS, ref WR, KS + o_wr, ref WI, KS + o_wi, ref H, offset_h, LDH, ILOZ
 955:                                           , IHIZ, ref Z, offset_z, LDZ, ref WORK, offset_work, 3, ref H, KU+1 * LDH + o_h
 956:                                           , LDH, NVE, ref H, KWV+1 * LDH + o_h, LDH, NHO, ref H, KU+KWH * LDH + o_h
 957:                                           , LDH);
 958:                      }
 959:                      // *
 960:                      // *           ==== Note progress (or the lack of it). ====
 961:                      // *
 962:                      if (LD > 0)
 963:                      {
 964:                          NDFL = 1;
 965:                      }
 966:                      else
 967:                      {
 968:                          NDFL = NDFL + 1;
 969:                      }
 970:                      // *
 971:                      // *           ==== End of main loop ====
 972:                  }
 973:                  // *
 974:                  // *        ==== Iteration limit exceeded.  Set INFO to show where
 975:                  // *        .    the problem occurred and exit. ====
 976:                  // *
 977:                  INFO = KBOT;
 978:              LABEL90:;
 979:              }
 980:              // *
 981:              // *     ==== Return the optimal value of LWORK. ====
 982:              // *
 983:              WORK[1 + o_work] = Convert.ToDouble(LWKOPT);
 984:              // *
 985:              // *     ==== End of DLAQR0 ====
 986:              // *
 987:   
 988:              #endregion
 989:   
 990:          }
 991:      }
 992:  }