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:      ///</summary>
  24:      public class DLAQR5
  25:      {
  26:      
  27:   
  28:          #region Dependencies
  29:          
  30:          DLAMCH _dlamch; DGEMM _dgemm; DLABAD _dlabad; DLACPY _dlacpy; DLAQR1 _dlaqr1; DLARFG _dlarfg; DLASET _dlaset; 
  31:          DTRMM _dtrmm;
  32:   
  33:          #endregion
  34:   
  35:   
  36:          #region Fields
  37:          
  38:          const double ZERO = 0.0E0; const double ONE = 1.0E0; double ALPHA = 0; double BETA = 0; double H11 = 0; double H12 = 0; 
  39:          double H21 = 0;double H22 = 0; double REFSUM = 0; double SAFMAX = 0; double SAFMIN = 0; double SCL = 0; double SMLNUM = 0; 
  40:          double SWAP = 0;double TST1 = 0; double TST2 = 0; double ULP = 0; int I = 0; int I2 = 0; int I4 = 0; int INCOL = 0; 
  41:          int J = 0;int J2 = 0; int J4 = 0; int JBOT = 0; int JCOL = 0; int JLEN = 0; int JROW = 0; int JTOP = 0; int K = 0; 
  42:          int K1 = 0;int KDU = 0; int KMS = 0; int KNZ = 0; int KRCOL = 0; int KZS = 0; int M = 0; int M22 = 0; int MBOT = 0; 
  43:          int MEND = 0;int MSTART = 0; int MTOP = 0; int NBMPS = 0; int NDCOL = 0; int NS = 0; int NU = 0; bool ACCUM = false; 
  44:          bool BLK22 = false;bool BMP22 = false; double[] VT = new double[3]; int offset_vt = 0; int o_vt = -1; 
  45:   
  46:          #endregion
  47:   
  48:          public DLAQR5(DLAMCH dlamch, DGEMM dgemm, DLABAD dlabad, DLACPY dlacpy, DLAQR1 dlaqr1, DLARFG dlarfg, DLASET dlaset, DTRMM dtrmm)
  49:          {
  50:      
  51:   
  52:              #region Set Dependencies
  53:              
  54:              this._dlamch = dlamch; this._dgemm = dgemm; this._dlabad = dlabad; this._dlacpy = dlacpy; this._dlaqr1 = dlaqr1; 
  55:              this._dlarfg = dlarfg;this._dlaset = dlaset; this._dtrmm = dtrmm; 
  56:   
  57:              #endregion
  58:   
  59:          }
  60:      
  61:          public DLAQR5()
  62:          {
  63:      
  64:   
  65:              #region Dependencies (Initialization)
  66:              
  67:              LSAME lsame = new LSAME();
  68:              DLAMC3 dlamc3 = new DLAMC3();
  69:              XERBLA xerbla = new XERBLA();
  70:              DLABAD dlabad = new DLABAD();
  71:              DLAQR1 dlaqr1 = new DLAQR1();
  72:              DLAPY2 dlapy2 = new DLAPY2();
  73:              DNRM2 dnrm2 = new DNRM2();
  74:              DSCAL dscal = new DSCAL();
  75:              DLAMC1 dlamc1 = new DLAMC1(dlamc3);
  76:              DLAMC4 dlamc4 = new DLAMC4(dlamc3);
  77:              DLAMC5 dlamc5 = new DLAMC5(dlamc3);
  78:              DLAMC2 dlamc2 = new DLAMC2(dlamc3, dlamc1, dlamc4, dlamc5);
  79:              DLAMCH dlamch = new DLAMCH(lsame, dlamc2);
  80:              DGEMM dgemm = new DGEMM(lsame, xerbla);
  81:              DLACPY dlacpy = new DLACPY(lsame);
  82:              DLARFG dlarfg = new DLARFG(dlamch, dlapy2, dnrm2, dscal);
  83:              DLASET dlaset = new DLASET(lsame);
  84:              DTRMM dtrmm = new DTRMM(lsame, xerbla);
  85:   
  86:              #endregion
  87:   
  88:   
  89:              #region Set Dependencies
  90:              
  91:              this._dlamch = dlamch; this._dgemm = dgemm; this._dlabad = dlabad; this._dlacpy = dlacpy; this._dlaqr1 = dlaqr1; 
  92:              this._dlarfg = dlarfg;this._dlaset = dlaset; this._dtrmm = dtrmm; 
  93:   
  94:              #endregion
  95:   
  96:          }
  97:          /// <param name="WANTT">
  98:          /// (input) logical scalar
  99:          /// WANTT = .true. if the quasi-triangular Schur factor
 100:          /// is being computed.  WANTT is set to .false. otherwise.
 101:          ///</param>
 102:          /// <param name="WANTZ">
 103:          /// (input) logical scalar
 104:          /// WANTZ = .true. if the orthogonal Schur factor is being
 105:          /// computed.  WANTZ is set to .false. otherwise.
 106:          ///</param>
 107:          /// <param name="KACC22">
 108:          /// (input) integer with value 0, 1, or 2.
 109:          /// Specifies the computation mode of far-from-diagonal
 110:          /// orthogonal updates.
 111:          /// = 0: DLAQR5 does not accumulate reflections and does not
 112:          /// use matrix-matrix multiply to update far-from-diagonal
 113:          /// matrix entries.
 114:          /// = 1: DLAQR5 accumulates reflections and uses matrix-matrix
 115:          /// multiply to update the far-from-diagonal matrix entries.
 116:          /// = 2: DLAQR5 accumulates reflections, uses matrix-matrix
 117:          /// multiply to update the far-from-diagonal matrix entries,
 118:          /// and takes advantage of 2-by-2 block structure during
 119:          /// matrix multiplies.
 120:          ///</param>
 121:          /// <param name="N">
 122:          /// (input) integer scalar
 123:          /// N is the order of the Hessenberg matrix H upon which this
 124:          /// subroutine operates.
 125:          ///</param>
 126:          /// <param name="KTOP">
 127:          /// (input) integer scalar
 128:          ///</param>
 129:          /// <param name="KBOT">
 130:          /// (input) integer scalar
 131:          /// These are the first and last rows and columns of an
 132:          /// isolated diagonal block upon which the QR sweep is to be
 133:          /// applied. It is assumed without a check that
 134:          /// either KTOP = 1  or   H(KTOP,KTOP-1) = 0
 135:          /// and
 136:          /// either KBOT = N  or   H(KBOT+1,KBOT) = 0.
 137:          ///</param>
 138:          /// <param name="NSHFTS">
 139:          /// (input) integer scalar
 140:          /// NSHFTS gives the number of simultaneous shifts.  NSHFTS
 141:          /// must be positive and even.
 142:          ///</param>
 143:          /// <param name="SR">
 144:          /// (input) DOUBLE PRECISION array of size (NSHFTS)
 145:          ///</param>
 146:          /// <param name="SI">
 147:          /// (input) DOUBLE PRECISION array of size (NSHFTS)
 148:          /// SR contains the real parts and SI contains the imaginary
 149:          /// parts of the NSHFTS shifts of origin that define the
 150:          /// multi-shift QR sweep.
 151:          ///</param>
 152:          /// <param name="H">
 153:          /// (input/output) DOUBLE PRECISION array of size (LDH,N)
 154:          /// On input H contains a Hessenberg matrix.  On output a
 155:          /// multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
 156:          /// to the isolated diagonal block in rows and columns KTOP
 157:          /// through KBOT.
 158:          ///</param>
 159:          /// <param name="LDH">
 160:          /// (input) integer scalar
 161:          /// LDH is the leading dimension of H just as declared in the
 162:          /// calling procedure.  LDH.GE.MAX(1,N).
 163:          ///</param>
 164:          /// <param name="ILOZ">
 165:          /// (input) INTEGER
 166:          ///</param>
 167:          /// <param name="IHIZ">
 168:          /// (input) INTEGER
 169:          /// Specify the rows of Z to which transformations must be
 170:          /// applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
 171:          ///</param>
 172:          /// <param name="Z">
 173:          /// (input/output) DOUBLE PRECISION array of size (LDZ,IHI)
 174:          /// If WANTZ = .TRUE., then the QR Sweep orthogonal
 175:          /// similarity transformation is accumulated into
 176:          /// Z(ILOZ:IHIZ,ILO:IHI) from the right.
 177:          /// If WANTZ = .FALSE., then Z is unreferenced.
 178:          ///</param>
 179:          /// <param name="LDZ">
 180:          /// (input) integer scalar
 181:          /// LDA is the leading dimension of Z just as declared in
 182:          /// the calling procedure. LDZ.GE.N.
 183:          ///</param>
 184:          /// <param name="V">
 185:          /// (workspace) DOUBLE PRECISION array of size (LDV,NSHFTS/2)
 186:          ///</param>
 187:          /// <param name="LDV">
 188:          /// (input) integer scalar
 189:          /// LDV is the leading dimension of V as declared in the
 190:          /// calling procedure.  LDV.GE.3.
 191:          ///</param>
 192:          /// <param name="U">
 193:          /// (workspace) DOUBLE PRECISION array of size
 194:          /// (LDU,3*NSHFTS-3)
 195:          ///</param>
 196:          /// <param name="LDU">
 197:          /// (input) integer scalar
 198:          /// LDU is the leading dimension of U just as declared in the
 199:          /// in the calling subroutine.  LDU.GE.3*NSHFTS-3.
 200:          ///</param>
 201:          /// <param name="NV">
 202:          /// (input) integer scalar
 203:          /// NV is the number of rows in WV agailable for workspace.
 204:          /// NV.GE.1.
 205:          ///</param>
 206:          /// <param name="WV">
 207:          /// (workspace) DOUBLE PRECISION array of size
 208:          /// (LDWV,3*NSHFTS-3)
 209:          ///</param>
 210:          /// <param name="LDWV">
 211:          /// (input) integer scalar
 212:          /// LDWV is the leading dimension of WV as declared in the
 213:          /// in the calling subroutine.  LDWV.GE.NV.
 214:          /// 
 215:          ///</param>
 216:          /// <param name="NH">
 217:          /// (input) integer scalar
 218:          /// NH is the number of columns in array WH available for
 219:          /// workspace. NH.GE.1.
 220:          ///</param>
 221:          /// <param name="WH">
 222:          /// (workspace) DOUBLE PRECISION array of size (LDWH,NH)
 223:          ///</param>
 224:          /// <param name="LDWH">
 225:          /// (input) integer scalar
 226:          /// Leading dimension of WH just as declared in the
 227:          /// calling procedure.  LDWH.GE.3*NSHFTS-3.
 228:          ///</param>
 229:          public void Run(bool WANTT, bool WANTZ, int KACC22, int N, int KTOP, int KBOT
 230:                           , int NSHFTS, ref double[] SR, int offset_sr, ref double[] SI, int offset_si, ref double[] H, int offset_h, int LDH, int ILOZ
 231:                           , int IHIZ, ref double[] Z, int offset_z, int LDZ, ref double[] V, int offset_v, int LDV, ref double[] U, int offset_u
 232:                           , int LDU, int NV, ref double[] WV, int offset_wv, int LDWV, int NH, ref double[] WH, int offset_wh
 233:                           , int LDWH)
 234:          {
 235:   
 236:              #region Array Index Correction
 237:              
 238:               int o_sr = -1 + offset_sr;  int o_si = -1 + offset_si;  int o_h = -1 - LDH + offset_h; 
 239:               int o_z = -1 - LDZ + offset_z; int o_v = -1 - LDV + offset_v;  int o_u = -1 - LDU + offset_u; 
 240:               int o_wv = -1 - LDWV + offset_wv; int o_wh = -1 - LDWH + offset_wh; 
 241:   
 242:              #endregion
 243:   
 244:   
 245:              #region Prolog
 246:              
 247:              // *
 248:              // *  -- LAPACK auxiliary routine (version 3.1) --
 249:              // *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 250:              // *     November 2006
 251:              // *
 252:              // *     .. Scalar Arguments ..
 253:              // *     ..
 254:              // *     .. Array Arguments ..
 255:              // *     ..
 256:              // *
 257:              // *     This auxiliary subroutine called by DLAQR0 performs a
 258:              // *     single small-bulge multi-shift QR sweep.
 259:              // *
 260:              // *      WANTT  (input) logical scalar
 261:              // *             WANTT = .true. if the quasi-triangular Schur factor
 262:              // *             is being computed.  WANTT is set to .false. otherwise.
 263:              // *
 264:              // *      WANTZ  (input) logical scalar
 265:              // *             WANTZ = .true. if the orthogonal Schur factor is being
 266:              // *             computed.  WANTZ is set to .false. otherwise.
 267:              // *
 268:              // *      KACC22 (input) integer with value 0, 1, or 2.
 269:              // *             Specifies the computation mode of far-from-diagonal
 270:              // *             orthogonal updates.
 271:              // *        = 0: DLAQR5 does not accumulate reflections and does not
 272:              // *             use matrix-matrix multiply to update far-from-diagonal
 273:              // *             matrix entries.
 274:              // *        = 1: DLAQR5 accumulates reflections and uses matrix-matrix
 275:              // *             multiply to update the far-from-diagonal matrix entries.
 276:              // *        = 2: DLAQR5 accumulates reflections, uses matrix-matrix
 277:              // *             multiply to update the far-from-diagonal matrix entries,
 278:              // *             and takes advantage of 2-by-2 block structure during
 279:              // *             matrix multiplies.
 280:              // *
 281:              // *      N      (input) integer scalar
 282:              // *             N is the order of the Hessenberg matrix H upon which this
 283:              // *             subroutine operates.
 284:              // *
 285:              // *      KTOP   (input) integer scalar
 286:              // *      KBOT   (input) integer scalar
 287:              // *             These are the first and last rows and columns of an
 288:              // *             isolated diagonal block upon which the QR sweep is to be
 289:              // *             applied. It is assumed without a check that
 290:              // *                       either KTOP = 1  or   H(KTOP,KTOP-1) = 0
 291:              // *             and
 292:              // *                       either KBOT = N  or   H(KBOT+1,KBOT) = 0.
 293:              // *
 294:              // *      NSHFTS (input) integer scalar
 295:              // *             NSHFTS gives the number of simultaneous shifts.  NSHFTS
 296:              // *             must be positive and even.
 297:              // *
 298:              // *      SR     (input) DOUBLE PRECISION array of size (NSHFTS)
 299:              // *      SI     (input) DOUBLE PRECISION array of size (NSHFTS)
 300:              // *             SR contains the real parts and SI contains the imaginary
 301:              // *             parts of the NSHFTS shifts of origin that define the
 302:              // *             multi-shift QR sweep.
 303:              // *
 304:              // *      H      (input/output) DOUBLE PRECISION array of size (LDH,N)
 305:              // *             On input H contains a Hessenberg matrix.  On output a
 306:              // *             multi-shift QR sweep with shifts SR(J)+i*SI(J) is applied
 307:              // *             to the isolated diagonal block in rows and columns KTOP
 308:              // *             through KBOT.
 309:              // *
 310:              // *      LDH    (input) integer scalar
 311:              // *             LDH is the leading dimension of H just as declared in the
 312:              // *             calling procedure.  LDH.GE.MAX(1,N).
 313:              // *
 314:              // *      ILOZ   (input) INTEGER
 315:              // *      IHIZ   (input) INTEGER
 316:              // *             Specify the rows of Z to which transformations must be
 317:              // *             applied if WANTZ is .TRUE.. 1 .LE. ILOZ .LE. IHIZ .LE. N
 318:              // *
 319:              // *      Z      (input/output) DOUBLE PRECISION array of size (LDZ,IHI)
 320:              // *             If WANTZ = .TRUE., then the QR Sweep orthogonal
 321:              // *             similarity transformation is accumulated into
 322:              // *             Z(ILOZ:IHIZ,ILO:IHI) from the right.
 323:              // *             If WANTZ = .FALSE., then Z is unreferenced.
 324:              // *
 325:              // *      LDZ    (input) integer scalar
 326:              // *             LDA is the leading dimension of Z just as declared in
 327:              // *             the calling procedure. LDZ.GE.N.
 328:              // *
 329:              // *      V      (workspace) DOUBLE PRECISION array of size (LDV,NSHFTS/2)
 330:              // *
 331:              // *      LDV    (input) integer scalar
 332:              // *             LDV is the leading dimension of V as declared in the
 333:              // *             calling procedure.  LDV.GE.3.
 334:              // *
 335:              // *      U      (workspace) DOUBLE PRECISION array of size
 336:              // *             (LDU,3*NSHFTS-3)
 337:              // *
 338:              // *      LDU    (input) integer scalar
 339:              // *             LDU is the leading dimension of U just as declared in the
 340:              // *             in the calling subroutine.  LDU.GE.3*NSHFTS-3.
 341:              // *
 342:              // *      NH     (input) integer scalar
 343:              // *             NH is the number of columns in array WH available for
 344:              // *             workspace. NH.GE.1.
 345:              // *
 346:              // *      WH     (workspace) DOUBLE PRECISION array of size (LDWH,NH)
 347:              // *
 348:              // *      LDWH   (input) integer scalar
 349:              // *             Leading dimension of WH just as declared in the
 350:              // *             calling procedure.  LDWH.GE.3*NSHFTS-3.
 351:              // *
 352:              // *      NV     (input) integer scalar
 353:              // *             NV is the number of rows in WV agailable for workspace.
 354:              // *             NV.GE.1.
 355:              // *
 356:              // *      WV     (workspace) DOUBLE PRECISION array of size
 357:              // *             (LDWV,3*NSHFTS-3)
 358:              // *
 359:              // *      LDWV   (input) integer scalar
 360:              // *             LDWV is the leading dimension of WV as declared in the
 361:              // *             in the calling subroutine.  LDWV.GE.NV.
 362:              // *
 363:              // *
 364:              // *     ================================================================
 365:              // *     Based on contributions by
 366:              // *        Karen Braman and Ralph Byers, Department of Mathematics,
 367:              // *        University of Kansas, USA
 368:              // *
 369:              // *     ============================================================
 370:              // *     Reference:
 371:              // *
 372:              // *     K. Braman, R. Byers and R. Mathias, The Multi-Shift QR
 373:              // *     Algorithm Part I: Maintaining Well Focused Shifts, and
 374:              // *     Level 3 Performance, SIAM Journal of Matrix Analysis,
 375:              // *     volume 23, pages 929--947, 2002.
 376:              // *
 377:              // *     ============================================================
 378:              // *     .. Parameters ..
 379:              // *     ..
 380:              // *     .. Local Scalars ..
 381:              // *     ..
 382:              // *     .. External Functions ..
 383:              // *     ..
 384:              // *     .. Intrinsic Functions ..
 385:              // *
 386:              //      INTRINSIC          ABS, DBLE, MAX, MIN, MOD;
 387:              // *     ..
 388:              // *     .. Local Arrays ..
 389:              // *     ..
 390:              // *     .. External Subroutines ..
 391:              // *     ..
 392:              // *     .. Executable Statements ..
 393:              // *
 394:              // *     ==== If there are no shifts, then there is nothing to do. ====
 395:              // *
 396:   
 397:              #endregion
 398:   
 399:   
 400:              #region Body
 401:              
 402:              if (NSHFTS < 2) return;
 403:              // *
 404:              // *     ==== If the active block is empty or 1-by-1, then there
 405:              // *     .    is nothing to do. ====
 406:              // *
 407:              if (KTOP >= KBOT) return;
 408:              // *
 409:              // *     ==== Shuffle shifts into pairs of real shifts and pairs
 410:              // *     .    of complex conjugate shifts assuming complex
 411:              // *     .    conjugate shifts are already adjacent to one
 412:              // *     .    another. ====
 413:              // *
 414:              for (I = 1; I <= NSHFTS - 2; I += 2)
 415:              {
 416:                  if (SI[I + o_si] !=  - SI[I + 1 + o_si])
 417:                  {
 418:                      // *
 419:                      SWAP = SR[I + o_sr];
 420:                      SR[I + o_sr] = SR[I + 1 + o_sr];
 421:                      SR[I + 1 + o_sr] = SR[I + 2 + o_sr];
 422:                      SR[I + 2 + o_sr] = SWAP;
 423:                      // *
 424:                      SWAP = SI[I + o_si];
 425:                      SI[I + o_si] = SI[I + 1 + o_si];
 426:                      SI[I + 1 + o_si] = SI[I + 2 + o_si];
 427:                      SI[I + 2 + o_si] = SWAP;
 428:                  }
 429:              }
 430:              // *
 431:              // *     ==== NSHFTS is supposed to be even, but if is odd,
 432:              // *     .    then simply reduce it by one.  The shuffle above
 433:              // *     .    ensures that the dropped shift is real and that
 434:              // *     .    the remaining shifts are paired. ====
 435:              // *
 436:              NS = NSHFTS - FortranLib.Mod(NSHFTS,2);
 437:              // *
 438:              // *     ==== Machine constants for deflation ====
 439:              // *
 440:              SAFMIN = this._dlamch.Run("SAFE MINIMUM");
 441:              SAFMAX = ONE / SAFMIN;
 442:              this._dlabad.Run(ref SAFMIN, ref SAFMAX);
 443:              ULP = this._dlamch.Run("PRECISION");
 444:              SMLNUM = SAFMIN * (Convert.ToDouble(N) / ULP);
 445:              // *
 446:              // *     ==== Use accumulated reflections to update far-from-diagonal
 447:              // *     .    entries ? ====
 448:              // *
 449:              ACCUM = (KACC22 == 1) || (KACC22 == 2);
 450:              // *
 451:              // *     ==== If so, exploit the 2-by-2 block structure? ====
 452:              // *
 453:              BLK22 = (NS > 2) && (KACC22 == 2);
 454:              // *
 455:              // *     ==== clear trash ====
 456:              // *
 457:              if (KTOP + 2 <= KBOT) H[KTOP + 2+KTOP * LDH + o_h] = ZERO;
 458:              // *
 459:              // *     ==== NBMPS = number of 2-shift bulges in the chain ====
 460:              // *
 461:              NBMPS = NS / 2;
 462:              // *
 463:              // *     ==== KDU = width of slab ====
 464:              // *
 465:              KDU = 6 * NBMPS - 3;
 466:              // *
 467:              // *     ==== Create and chase chains of NBMPS bulges ====
 468:              // *
 469:              for (INCOL = 3 * (1 - NBMPS) + KTOP - 1; (3 * NBMPS - 2 >= 0) ? (INCOL <= KBOT - 2) : (INCOL >= KBOT - 2); INCOL += 3 * NBMPS - 2)
 470:              {
 471:                  NDCOL = INCOL + KDU;
 472:                  if (ACCUM)
 473:                  {
 474:                      this._dlaset.Run("ALL", KDU, KDU, ZERO, ONE, ref U, offset_u
 475:                                       , LDU);
 476:                  }
 477:                  // *
 478:                  // *        ==== Near-the-diagonal bulge chase.  The following loop
 479:                  // *        .    performs the near-the-diagonal part of a small bulge
 480:                  // *        .    multi-shift QR sweep.  Each 6*NBMPS-2 column diagonal
 481:                  // *        .    chunk extends from column INCOL to column NDCOL
 482:                  // *        .    (including both column INCOL and column NDCOL). The
 483:                  // *        .    following loop chases a 3*NBMPS column long chain of
 484:                  // *        .    NBMPS bulges 3*NBMPS-2 columns to the right.  (INCOL
 485:                  // *        .    may be less than KTOP and and NDCOL may be greater than
 486:                  // *        .    KBOT indicating phantom columns from which to chase
 487:                  // *        .    bulges before they are actually introduced or to which
 488:                  // *        .    to chase bulges beyond column KBOT.)  ====
 489:                  // *
 490:                  for (KRCOL = INCOL; KRCOL <= Math.Min(INCOL + 3 * NBMPS - 3, KBOT - 2); KRCOL++)
 491:                  {
 492:                      // *
 493:                      // *           ==== Bulges number MTOP to MBOT are active double implicit
 494:                      // *           .    shift bulges.  There may or may not also be small
 495:                      // *           .    2-by-2 bulge, if there is room.  The inactive bulges
 496:                      // *           .    (if any) must wait until the active bulges have moved
 497:                      // *           .    down the diagonal to make room.  The phantom matrix
 498:                      // *           .    paradigm described above helps keep track.  ====
 499:                      // *
 500:                      MTOP = Math.Max(1, ((KTOP - 1) - KRCOL + 2) / 3 + 1);
 501:                      MBOT = Math.Min(NBMPS, (KBOT - KRCOL) / 3);
 502:                      M22 = MBOT + 1;
 503:                      BMP22 = (MBOT < NBMPS) && (KRCOL + 3 * (M22 - 1)) == (KBOT - 2);
 504:                      // *
 505:                      // *           ==== Generate reflections to chase the chain right
 506:                      // *           .    one column.  (The minimum value of K is KTOP-1.) ====
 507:                      // *
 508:                      for (M = MTOP; M <= MBOT; M++)
 509:                      {
 510:                          K = KRCOL + 3 * (M - 1);
 511:                          if (K == KTOP - 1)
 512:                          {
 513:                              this._dlaqr1.Run(3, H, KTOP+KTOP * LDH + o_h, LDH, SR[2 * M - 1 + o_sr], SI[2 * M - 1 + o_si], SR[2 * M + o_sr]
 514:                                               , SI[2 * M + o_si], ref V, 1+M * LDV + o_v);
 515:                              ALPHA = V[1+M * LDV + o_v];
 516:                              this._dlarfg.Run(3, ref ALPHA, ref V, 2+M * LDV + o_v, 1, ref V[1+M * LDV + o_v]);
 517:                          }
 518:                          else
 519:                          {
 520:                              BETA = H[K + 1+K * LDH + o_h];
 521:                              V[2+M * LDV + o_v] = H[K + 2+K * LDH + o_h];
 522:                              V[3+M * LDV + o_v] = H[K + 3+K * LDH + o_h];
 523:                              this._dlarfg.Run(3, ref BETA, ref V, 2+M * LDV + o_v, 1, ref V[1+M * LDV + o_v]);
 524:                              // *
 525:                              // *                 ==== A Bulge may collapse because of vigilant
 526:                              // *                 .    deflation or destructive underflow.  (The
 527:                              // *                 .    initial bulge is always collapsed.) Use
 528:                              // *                 .    the two-small-subdiagonals trick to try
 529:                              // *                 .    to get it started again. If V(2,M).NE.0 and
 530:                              // *                 .    V(3,M) = H(K+3,K+1) = H(K+3,K+2) = 0, then
 531:                              // *                 .    this bulge is collapsing into a zero
 532:                              // *                 .    subdiagonal.  It will be restarted next
 533:                              // *                 .    trip through the loop.)
 534:                              // *
 535:                              if (V[1+M * LDV + o_v] != ZERO && (V[3+M * LDV + o_v] != ZERO || (H[K + 3+(K + 1) * LDH + o_h] == ZERO && H[K + 3+(K + 2) * LDH + o_h] == ZERO)))
 536:                              {
 537:                                  // *
 538:                                  // *                    ==== Typical case: not collapsed (yet). ====
 539:                                  // *
 540:                                  H[K + 1+K * LDH + o_h] = BETA;
 541:                                  H[K + 2+K * LDH + o_h] = ZERO;
 542:                                  H[K + 3+K * LDH + o_h] = ZERO;
 543:                              }
 544:                              else
 545:                              {
 546:                                  // *
 547:                                  // *                    ==== Atypical case: collapsed.  Attempt to
 548:                                  // *                    .    reintroduce ignoring H(K+1,K).  If the
 549:                                  // *                    .    fill resulting from the new reflector
 550:                                  // *                    .    is too large, then abandon it.
 551:                                  // *                    .    Otherwise, use the new one. ====
 552:                                  // *
 553:                                  this._dlaqr1.Run(3, H, K + 1+(K + 1) * LDH + o_h, LDH, SR[2 * M - 1 + o_sr], SI[2 * M - 1 + o_si], SR[2 * M + o_sr]
 554:                                                   , SI[2 * M + o_si], ref VT, offset_vt);
 555:                                  SCL = Math.Abs(VT[1 + o_vt]) + Math.Abs(VT[2 + o_vt]) + Math.Abs(VT[3 + o_vt]);
 556:                                  if (SCL != ZERO)
 557:                                  {
 558:                                      VT[1 + o_vt] = VT[1 + o_vt] / SCL;
 559:                                      VT[2 + o_vt] = VT[2 + o_vt] / SCL;
 560:                                      VT[3 + o_vt] = VT[3 + o_vt] / SCL;
 561:                                  }
 562:                                  // *
 563:                                  // *                    ==== The following is the traditional and
 564:                                  // *                    .    conservative two-small-subdiagonals
 565:                                  // *                    .    test.  ====
 566:                                  // *                    .
 567:                                  if (Math.Abs(H[K + 1+K * LDH + o_h]) * (Math.Abs(VT[2 + o_vt]) + Math.Abs(VT[3 + o_vt])) > ULP * Math.Abs(VT[1 + o_vt]) * (Math.Abs(H[K+K * LDH + o_h]) + Math.Abs(H[K + 1+(K + 1) * LDH + o_h]) + Math.Abs(H[K + 2+(K + 2) * LDH + o_h])))
 568:                                  {
 569:                                      // *
 570:                                      // *                       ==== Starting a new bulge here would
 571:                                      // *                       .    create non-negligible fill.   If
 572:                                      // *                       .    the old reflector is diagonal (only
 573:                                      // *                       .    possible with underflows), then
 574:                                      // *                       .    change it to I.  Otherwise, use
 575:                                      // *                       .    it with trepidation. ====
 576:                                      // *
 577:                                      if (V[2+M * LDV + o_v] == ZERO && V[3+M * LDV + o_v] == ZERO)
 578:                                      {
 579:                                          V[1+M * LDV + o_v] = ZERO;
 580:                                      }
 581:                                      else
 582:                                      {
 583:                                          H[K + 1+K * LDH + o_h] = BETA;
 584:                                          H[K + 2+K * LDH + o_h] = ZERO;
 585:                                          H[K + 3+K * LDH + o_h] = ZERO;
 586:                                      }
 587:                                  }
 588:                                  else
 589:                                  {
 590:                                      // *
 591:                                      // *                       ==== Stating a new bulge here would
 592:                                      // *                       .    create only negligible fill.
 593:                                      // *                       .    Replace the old reflector with
 594:                                      // *                       .    the new one. ====
 595:                                      // *
 596:                                      ALPHA = VT[1 + o_vt];
 597:                                      this._dlarfg.Run(3, ref ALPHA, ref VT, 2 + o_vt, 1, ref VT[1 + o_vt]);
 598:                                      REFSUM = H[K + 1+K * LDH + o_h] + H[K + 2+K * LDH + o_h] * VT[2 + o_vt] + H[K + 3+K * LDH + o_h] * VT[3 + o_vt];
 599:                                      H[K + 1+K * LDH + o_h] = H[K + 1+K * LDH + o_h] - VT[1 + o_vt] * REFSUM;
 600:                                      H[K + 2+K * LDH + o_h] = ZERO;
 601:                                      H[K + 3+K * LDH + o_h] = ZERO;
 602:                                      V[1+M * LDV + o_v] = VT[1 + o_vt];
 603:                                      V[2+M * LDV + o_v] = VT[2 + o_vt];
 604:                                      V[3+M * LDV + o_v] = VT[3 + o_vt];
 605:                                  }
 606:                              }
 607:                          }
 608:                      }
 609:                      // *
 610:                      // *           ==== Generate a 2-by-2 reflection, if needed. ====
 611:                      // *
 612:                      K = KRCOL + 3 * (M22 - 1);
 613:                      if (BMP22)
 614:                      {
 615:                          if (K == KTOP - 1)
 616:                          {
 617:                              this._dlaqr1.Run(2, H, K + 1+(K + 1) * LDH + o_h, LDH, SR[2 * M22 - 1 + o_sr], SI[2 * M22 - 1 + o_si], SR[2 * M22 + o_sr]
 618:                                               , SI[2 * M22 + o_si], ref V, 1+M22 * LDV + o_v);
 619:                              BETA = V[1+M22 * LDV + o_v];
 620:                              this._dlarfg.Run(2, ref BETA, ref V, 2+M22 * LDV + o_v, 1, ref V[1+M22 * LDV + o_v]);
 621:                          }
 622:                          else
 623:                          {
 624:                              BETA = H[K + 1+K * LDH + o_h];
 625:                              V[2+M22 * LDV + o_v] = H[K + 2+K * LDH + o_h];
 626:                              this._dlarfg.Run(2, ref BETA, ref V, 2+M22 * LDV + o_v, 1, ref V[1+M22 * LDV + o_v]);
 627:                              H[K + 1+K * LDH + o_h] = BETA;
 628:                              H[K + 2+K * LDH + o_h] = ZERO;
 629:                          }
 630:                      }
 631:                      else
 632:                      {
 633:                          // *
 634:                          // *              ==== Initialize V(1,M22) here to avoid possible undefined
 635:                          // *              .    variable problems later. ====
 636:                          // *
 637:                          V[1+M22 * LDV + o_v] = ZERO;
 638:                      }
 639:                      // *
 640:                      // *           ==== Multiply H by reflections from the left ====
 641:                      // *
 642:                      if (ACCUM)
 643:                      {
 644:                          JBOT = Math.Min(NDCOL, KBOT);
 645:                      }
 646:                      else
 647:                      {
 648:                          if (WANTT)
 649:                          {
 650:                              JBOT = N;
 651:                          }
 652:                          else
 653:                          {
 654:                              JBOT = KBOT;
 655:                          }
 656:                      }
 657:                      for (J = Math.Max(KTOP, KRCOL); J <= JBOT; J++)
 658:                      {
 659:                          MEND = Math.Min(MBOT, (J - KRCOL + 2) / 3);
 660:                          for (M = MTOP; M <= MEND; M++)
 661:                          {
 662:                              K = KRCOL + 3 * (M - 1);
 663:                              REFSUM = V[1+M * LDV + o_v] * (H[K + 1+J * LDH + o_h] + V[2+M * LDV + o_v] * H[K + 2+J * LDH + o_h] + V[3+M * LDV + o_v] * H[K + 3+J * LDH + o_h]);
 664:                              H[K + 1+J * LDH + o_h] = H[K + 1+J * LDH + o_h] - REFSUM;
 665:                              H[K + 2+J * LDH + o_h] = H[K + 2+J * LDH + o_h] - REFSUM * V[2+M * LDV + o_v];
 666:                              H[K + 3+J * LDH + o_h] = H[K + 3+J * LDH + o_h] - REFSUM * V[3+M * LDV + o_v];
 667:                          }
 668:                      }
 669:                      if (BMP22)
 670:                      {
 671:                          K = KRCOL + 3 * (M22 - 1);
 672:                          for (J = Math.Max(K + 1, KTOP); J <= JBOT; J++)
 673:                          {
 674:                              REFSUM = V[1+M22 * LDV + o_v] * (H[K + 1+J * LDH + o_h] + V[2+M22 * LDV + o_v] * H[K + 2+J * LDH + o_h]);
 675:                              H[K + 1+J * LDH + o_h] = H[K + 1+J * LDH + o_h] - REFSUM;
 676:                              H[K + 2+J * LDH + o_h] = H[K + 2+J * LDH + o_h] - REFSUM * V[2+M22 * LDV + o_v];
 677:                          }
 678:                      }
 679:                      // *
 680:                      // *           ==== Multiply H by reflections from the right.
 681:                      // *           .    Delay filling in the last row until the
 682:                      // *           .    vigilant deflation check is complete. ====
 683:                      // *
 684:                      if (ACCUM)
 685:                      {
 686:                          JTOP = Math.Max(KTOP, INCOL);
 687:                      }
 688:                      else
 689:                      {
 690:                          if (WANTT)
 691:                          {
 692:                              JTOP = 1;
 693:                          }
 694:                          else
 695:                          {
 696:                              JTOP = KTOP;
 697:                          }
 698:                      }
 699:                      for (M = MTOP; M <= MBOT; M++)
 700:                      {
 701:                          if (V[1+M * LDV + o_v] != ZERO)
 702:                          {
 703:                              K = KRCOL + 3 * (M - 1);
 704:                              for (J = JTOP; J <= Math.Min(KBOT, K + 3); J++)
 705:                              {
 706:                                  REFSUM = V[1+M * LDV + o_v] * (H[J+(K + 1) * LDH + o_h] + V[2+M * LDV + o_v] * H[J+(K + 2) * LDH + o_h] + V[3+M * LDV + o_v] * H[J+(K + 3) * LDH + o_h]);
 707:                                  H[J+(K + 1) * LDH + o_h] = H[J+(K + 1) * LDH + o_h] - REFSUM;
 708:                                  H[J+(K + 2) * LDH + o_h] = H[J+(K + 2) * LDH + o_h] - REFSUM * V[2+M * LDV + o_v];
 709:                                  H[J+(K + 3) * LDH + o_h] = H[J+(K + 3) * LDH + o_h] - REFSUM * V[3+M * LDV + o_v];
 710:                              }
 711:                              // *
 712:                              if (ACCUM)
 713:                              {
 714:                                  // *
 715:                                  // *                    ==== Accumulate U. (If necessary, update Z later
 716:                                  // *                    .    with with an efficient matrix-matrix
 717:                                  // *                    .    multiply.) ====
 718:                                  // *
 719:                                  KMS = K - INCOL;
 720:                                  for (J = Math.Max(1, KTOP - INCOL); J <= KDU; J++)
 721:                                  {
 722:                                      REFSUM = V[1+M * LDV + o_v] * (U[J+(KMS + 1) * LDU + o_u] + V[2+M * LDV + o_v] * U[J+(KMS + 2) * LDU + o_u] + V[3+M * LDV + o_v] * U[J+(KMS + 3) * LDU + o_u]);
 723:                                      U[J+(KMS + 1) * LDU + o_u] = U[J+(KMS + 1) * LDU + o_u] - REFSUM;
 724:                                      U[J+(KMS + 2) * LDU + o_u] = U[J+(KMS + 2) * LDU + o_u] - REFSUM * V[2+M * LDV + o_v];
 725:                                      U[J+(KMS + 3) * LDU + o_u] = U[J+(KMS + 3) * LDU + o_u] - REFSUM * V[3+M * LDV + o_v];
 726:                                  }
 727:                              }
 728:                              else
 729:                              {
 730:                                  if (WANTZ)
 731:                                  {
 732:                                      // *
 733:                                      // *                    ==== U is not accumulated, so update Z
 734:                                      // *                    .    now by multiplying by reflections
 735:                                      // *                    .    from the right. ====
 736:                                      // *
 737:                                      for (J = ILOZ; J <= IHIZ; J++)
 738:                                      {
 739:                                          REFSUM = V[1+M * LDV + o_v] * (Z[J+(K + 1) * LDZ + o_z] + V[2+M * LDV + o_v] * Z[J+(K + 2) * LDZ + o_z] + V[3+M * LDV + o_v] * Z[J+(K + 3) * LDZ + o_z]);
 740:                                          Z[J+(K + 1) * LDZ + o_z] = Z[J+(K + 1) * LDZ + o_z] - REFSUM;
 741:                                          Z[J+(K + 2) * LDZ + o_z] = Z[J+(K + 2) * LDZ + o_z] - REFSUM * V[2+M * LDV + o_v];
 742:                                          Z[J+(K + 3) * LDZ + o_z] = Z[J+(K + 3) * LDZ + o_z] - REFSUM * V[3+M * LDV + o_v];
 743:                                      }
 744:                                  }
 745:                              }
 746:                          }
 747:                      }
 748:                      // *
 749:                      // *           ==== Special case: 2-by-2 reflection (if needed) ====
 750:                      // *
 751:                      K = KRCOL + 3 * (M22 - 1);
 752:                      if (BMP22 && (V[1+M22 * LDV + o_v] != ZERO))
 753:                      {
 754:                          for (J = JTOP; J <= Math.Min(KBOT, K + 3); J++)
 755:                          {
 756:                              REFSUM = V[1+M22 * LDV + o_v] * (H[J+(K + 1) * LDH + o_h] + V[2+M22 * LDV + o_v] * H[J+(K + 2) * LDH + o_h]);
 757:                              H[J+(K + 1) * LDH + o_h] = H[J+(K + 1) * LDH + o_h] - REFSUM;
 758:                              H[J+(K + 2) * LDH + o_h] = H[J+(K + 2) * LDH + o_h] - REFSUM * V[2+M22 * LDV + o_v];
 759:                          }
 760:                          // *
 761:                          if (ACCUM)
 762:                          {
 763:                              KMS = K - INCOL;
 764:                              for (J = Math.Max(1, KTOP - INCOL); J <= KDU; J++)
 765:                              {
 766:                                  REFSUM = V[1+M22 * LDV + o_v] * (U[J+(KMS + 1) * LDU + o_u] + V[2+M22 * LDV + o_v] * U[J+(KMS + 2) * LDU + o_u]);
 767:                                  U[J+(KMS + 1) * LDU + o_u] = U[J+(KMS + 1) * LDU + o_u] - REFSUM;
 768:                                  U[J+(KMS + 2) * LDU + o_u] = U[J+(KMS + 2) * LDU + o_u] - REFSUM * V[2+M22 * LDV + o_v];
 769:                              }
 770:                          }
 771:                          else
 772:                          {
 773:                              if (WANTZ)
 774:                              {
 775:                                  for (J = ILOZ; J <= IHIZ; J++)
 776:                                  {
 777:                                      REFSUM = V[1+M22 * LDV + o_v] * (Z[J+(K + 1) * LDZ + o_z] + V[2+M22 * LDV + o_v] * Z[J+(K + 2) * LDZ + o_z]);
 778:                                      Z[J+(K + 1) * LDZ + o_z] = Z[J+(K + 1) * LDZ + o_z] - REFSUM;
 779:                                      Z[J+(K + 2) * LDZ + o_z] = Z[J+(K + 2) * LDZ + o_z] - REFSUM * V[2+M22 * LDV + o_v];
 780:                                  }
 781:                              }
 782:                          }
 783:                      }
 784:                      // *
 785:                      // *           ==== Vigilant deflation check ====
 786:                      // *
 787:                      MSTART = MTOP;
 788:                      if (KRCOL + 3 * (MSTART - 1) < KTOP) MSTART = MSTART + 1;
 789:                      MEND = MBOT;
 790:                      if (BMP22) MEND = MEND + 1;
 791:                      if (KRCOL == KBOT - 2) MEND = MEND + 1;
 792:                      for (M = MSTART; M <= MEND; M++)
 793:                      {
 794:                          K = Math.Min(KBOT - 1, KRCOL + 3 * (M - 1));
 795:                          // *
 796:                          // *              ==== The following convergence test requires that
 797:                          // *              .    the tradition small-compared-to-nearby-diagonals
 798:                          // *              .    criterion and the Ahues & Tisseur (LAWN 122, 1997)
 799:                          // *              .    criteria both be satisfied.  The latter improves
 800:                          // *              .    accuracy in some examples. Falling back on an
 801:                          // *              .    alternate convergence criterion when TST1 or TST2
 802:                          // *              .    is zero (as done here) is traditional but probably
 803:                          // *              .    unnecessary. ====
 804:                          // *
 805:                          if (H[K + 1+K * LDH + o_h] != ZERO)
 806:                          {
 807:                              TST1 = Math.Abs(H[K+K * LDH + o_h]) + Math.Abs(H[K + 1+(K + 1) * LDH + o_h]);
 808:                              if (TST1 == ZERO)
 809:                              {
 810:                                  if (K >= KTOP + 1) TST1 = TST1 + Math.Abs(H[K+(K - 1) * LDH + o_h]);
 811:                                  if (K >= KTOP + 2) TST1 = TST1 + Math.Abs(H[K+(K - 2) * LDH + o_h]);
 812:                                  if (K >= KTOP + 3) TST1 = TST1 + Math.Abs(H[K+(K - 3) * LDH + o_h]);
 813:                                  if (K <= KBOT - 2) TST1 = TST1 + Math.Abs(H[K + 2+(K + 1) * LDH + o_h]);
 814:                                  if (K <= KBOT - 3) TST1 = TST1 + Math.Abs(H[K + 3+(K + 1) * LDH + o_h]);
 815:                                  if (K <= KBOT - 4) TST1 = TST1 + Math.Abs(H[K + 4+(K + 1) * LDH + o_h]);
 816:                              }
 817:                              if (Math.Abs(H[K + 1+K * LDH + o_h]) <= Math.Max(SMLNUM, ULP * TST1))
 818:                              {
 819:                                  H12 = Math.Max(Math.Abs(H[K + 1+K * LDH + o_h]), Math.Abs(H[K+(K + 1) * LDH + o_h]));
 820:                                  H21 = Math.Min(Math.Abs(H[K + 1+K * LDH + o_h]), Math.Abs(H[K+(K + 1) * LDH + o_h]));
 821:                                  H11 = Math.Max(Math.Abs(H[K + 1+(K + 1) * LDH + o_h]), Math.Abs(H[K+K * LDH + o_h] - H[K + 1+(K + 1) * LDH + o_h]));
 822:                                  H22 = Math.Min(Math.Abs(H[K + 1+(K + 1) * LDH + o_h]), Math.Abs(H[K+K * LDH + o_h] - H[K + 1+(K + 1) * LDH + o_h]));
 823:                                  SCL = H11 + H12;
 824:                                  TST2 = H22 * (H11 / SCL);
 825:                                  // *
 826:                                  if (TST2 == ZERO || H21 * (H12 / SCL) <= Math.Max(SMLNUM, ULP * TST2)) H[K + 1+K * LDH + o_h] = ZERO;
 827:                              }
 828:                          }
 829:                      }
 830:                      // *
 831:                      // *           ==== Fill in the last row of each bulge. ====
 832:                      // *
 833:                      MEND = Math.Min(NBMPS, (KBOT - KRCOL - 1) / 3);
 834:                      for (M = MTOP; M <= MEND; M++)
 835:                      {
 836:                          K = KRCOL + 3 * (M - 1);
 837:                          REFSUM = V[1+M * LDV + o_v] * V[3+M * LDV + o_v] * H[K + 4+(K + 3) * LDH + o_h];
 838:                          H[K + 4+(K + 1) * LDH + o_h] =  - REFSUM;
 839:                          H[K + 4+(K + 2) * LDH + o_h] =  - REFSUM * V[2+M * LDV + o_v];
 840:                          H[K + 4+(K + 3) * LDH + o_h] = H[K + 4+(K + 3) * LDH + o_h] - REFSUM * V[3+M * LDV + o_v];
 841:                      }
 842:                      // *
 843:                      // *           ==== End of near-the-diagonal bulge chase. ====
 844:                      // *
 845:                  }
 846:                  // *
 847:                  // *        ==== Use U (if accumulated) to update far-from-diagonal
 848:                  // *        .    entries in H.  If required, use U to update Z as
 849:                  // *        .    well. ====
 850:                  // *
 851:                  if (ACCUM)
 852:                  {
 853:                      if (WANTT)
 854:                      {
 855:                          JTOP = 1;
 856:                          JBOT = N;
 857:                      }
 858:                      else
 859:                      {
 860:                          JTOP = KTOP;
 861:                          JBOT = KBOT;
 862:                      }
 863:                      if ((!BLK22) || (INCOL < KTOP) || (NDCOL > KBOT) || (NS <= 2))
 864:                      {
 865:                          // *
 866:                          // *              ==== Updates not exploiting the 2-by-2 block
 867:                          // *              .    structure of U.  K1 and NU keep track of
 868:                          // *              .    the location and size of U in the special
 869:                          // *              .    cases of introducing bulges and chasing
 870:                          // *              .    bulges off the bottom.  In these special
 871:                          // *              .    cases and in case the number of shifts
 872:                          // *              .    is NS = 2, there is no 2-by-2 block
 873:                          // *              .    structure to exploit.  ====
 874:                          // *
 875:                          K1 = Math.Max(1, KTOP - INCOL);
 876:                          NU = (KDU - Math.Max(0, NDCOL - KBOT)) - K1 + 1;
 877:                          // *
 878:                          // *              ==== Horizontal Multiply ====
 879:                          // *
 880:                          for (JCOL = Math.Min(NDCOL, KBOT) + 1; (NH >= 0) ? (JCOL <= JBOT) : (JCOL >= JBOT); JCOL += NH)
 881:                          {
 882:                              JLEN = Math.Min(NH, JBOT - JCOL + 1);
 883:                              this._dgemm.Run("C", "N", NU, JLEN, NU, ONE
 884:                                              , U, K1+K1 * LDU + o_u, LDU, H, INCOL + K1+JCOL * LDH + o_h, LDH, ZERO, ref WH, offset_wh
 885:                                              , LDWH);
 886:                              this._dlacpy.Run("ALL", NU, JLEN, WH, offset_wh, LDWH, ref H, INCOL + K1+JCOL * LDH + o_h
 887:                                               , LDH);
 888:                          }
 889:                          // *
 890:                          // *              ==== Vertical multiply ====
 891:                          // *
 892:                          for (JROW = JTOP; (NV >= 0) ? (JROW <= Math.Max(KTOP, INCOL) - 1) : (JROW >= Math.Max(KTOP, INCOL) - 1); JROW += NV)
 893:                          {
 894:                              JLEN = Math.Min(NV, Math.Max(KTOP, INCOL) - JROW);
 895:                              this._dgemm.Run("N", "N", JLEN, NU, NU, ONE
 896:                                              , H, JROW+(INCOL + K1) * LDH + o_h, LDH, U, K1+K1 * LDU + o_u, LDU, ZERO, ref WV, offset_wv
 897:                                              , LDWV);
 898:                              this._dlacpy.Run("ALL", JLEN, NU, WV, offset_wv, LDWV, ref H, JROW+(INCOL + K1) * LDH + o_h
 899:                                               , LDH);
 900:                          }
 901:                          // *
 902:                          // *              ==== Z multiply (also vertical) ====
 903:                          // *
 904:                          if (WANTZ)
 905:                          {
 906:                              for (JROW = ILOZ; (NV >= 0) ? (JROW <= IHIZ) : (JROW >= IHIZ); JROW += NV)
 907:                              {
 908:                                  JLEN = Math.Min(NV, IHIZ - JROW + 1);
 909:                                  this._dgemm.Run("N", "N", JLEN, NU, NU, ONE
 910:                                                  , Z, JROW+(INCOL + K1) * LDZ + o_z, LDZ, U, K1+K1 * LDU + o_u, LDU, ZERO, ref WV, offset_wv
 911:                                                  , LDWV);
 912:                                  this._dlacpy.Run("ALL", JLEN, NU, WV, offset_wv, LDWV, ref Z, JROW+(INCOL + K1) * LDZ + o_z
 913:                                                   , LDZ);
 914:                              }
 915:                          }
 916:                      }
 917:                      else
 918:                      {
 919:                          // *
 920:                          // *              ==== Updates exploiting U's 2-by-2 block structure.
 921:                          // *              .    (I2, I4, J2, J4 are the last rows and columns
 922:                          // *              .    of the blocks.) ====
 923:                          // *
 924:                          I2 = (KDU + 1) / 2;
 925:                          I4 = KDU;
 926:                          J2 = I4 - I2;
 927:                          J4 = KDU;
 928:                          // *
 929:                          // *              ==== KZS and KNZ deal with the band of zeros
 930:                          // *              .    along the diagonal of one of the triangular
 931:                          // *              .    blocks. ====
 932:                          // *
 933:                          KZS = (J4 - J2) - (NS + 1);
 934:                          KNZ = NS + 1;
 935:                          // *
 936:                          // *              ==== Horizontal multiply ====
 937:                          // *
 938:                          for (JCOL = Math.Min(NDCOL, KBOT) + 1; (NH >= 0) ? (JCOL <= JBOT) : (JCOL >= JBOT); JCOL += NH)
 939:                          {
 940:                              JLEN = Math.Min(NH, JBOT - JCOL + 1);
 941:                              // *
 942:                              // *                 ==== Copy bottom of H to top+KZS of scratch ====
 943:                              // *                  (The first KZS rows get multiplied by zero.) ====
 944:                              // *
 945:                              this._dlacpy.Run("ALL", KNZ, JLEN, H, INCOL + 1 + J2+JCOL * LDH + o_h, LDH, ref WH, KZS + 1+1 * LDWH + o_wh
 946:                                               , LDWH);
 947:                              // *
 948:                              // *                 ==== Multiply by U21' ====
 949:                              // *
 950:                              this._dlaset.Run("ALL", KZS, JLEN, ZERO, ZERO, ref WH, offset_wh
 951:                                               , LDWH);
 952:                              this._dtrmm.Run("L", "U", "C", "N", KNZ, JLEN
 953:                                              , ONE, U, J2 + 1+(1 + KZS) * LDU + o_u, LDU, ref WH, KZS + 1+1 * LDWH + o_wh, LDWH);
 954:                              // *
 955:                              // *                 ==== Multiply top of H by U11' ====
 956:                              // *
 957:                              this._dgemm.Run("C", "N", I2, JLEN, J2, ONE
 958:                                              , U, offset_u, LDU, H, INCOL + 1+JCOL * LDH + o_h, LDH, ONE, ref WH, offset_wh
 959:                                              , LDWH);
 960:                              // *
 961:                              // *                 ==== Copy top of H bottom of WH ====
 962:                              // *
 963:                              this._dlacpy.Run("ALL", J2, JLEN, H, INCOL + 1+JCOL * LDH + o_h, LDH, ref WH, I2 + 1+1 * LDWH + o_wh
 964:                                               , LDWH);
 965:                              // *
 966:                              // *                 ==== Multiply by U21' ====
 967:                              // *
 968:                              this._dtrmm.Run("L", "L", "C", "N", J2, JLEN
 969:                                              , ONE, U, 1+(I2 + 1) * LDU + o_u, LDU, ref WH, I2 + 1+1 * LDWH + o_wh, LDWH);
 970:                              // *
 971:                              // *                 ==== Multiply by U22 ====
 972:                              // *
 973:                              this._dgemm.Run("C", "N", I4 - I2, JLEN, J4 - J2, ONE
 974:                                              , U, J2 + 1+(I2 + 1) * LDU + o_u, LDU, H, INCOL + 1 + J2+JCOL * LDH + o_h, LDH, ONE, ref WH, I2 + 1+1 * LDWH + o_wh
 975:                                              , LDWH);
 976:                              // *
 977:                              // *                 ==== Copy it back ====
 978:                              // *
 979:                              this._dlacpy.Run("ALL", KDU, JLEN, WH, offset_wh, LDWH, ref H, INCOL + 1+JCOL * LDH + o_h
 980:                                               , LDH);
 981:                          }
 982:                          // *
 983:                          // *              ==== Vertical multiply ====
 984:                          // *
 985:                          for (JROW = JTOP; (NV >= 0) ? (JROW <= Math.Max(INCOL, KTOP) - 1) : (JROW >= Math.Max(INCOL, KTOP) - 1); JROW += NV)
 986:                          {
 987:                              JLEN = Math.Min(NV, Math.Max(INCOL, KTOP) - JROW);
 988:                              // *
 989:                              // *                 ==== Copy right of H to scratch (the first KZS
 990:                              // *                 .    columns get multiplied by zero) ====
 991:                              // *
 992:                              this._dlacpy.Run("ALL", JLEN, KNZ, H, JROW+(INCOL + 1 + J2) * LDH + o_h, LDH, ref WV, 1+(1 + KZS) * LDWV + o_wv
 993:                                               , LDWV);
 994:                              // *
 995:                              // *                 ==== Multiply by U21 ====
 996:                              // *
 997:                              this._dlaset.Run("ALL", JLEN, KZS, ZERO, ZERO, ref WV, offset_wv
 998:                                               , LDWV);
 999:                              this._dtrmm.Run("R", "U", "N", "N", JLEN, KNZ
1000:                                              , ONE, U, J2 + 1+(1 + KZS) * LDU + o_u, LDU, ref WV, 1+(1 + KZS) * LDWV + o_wv, LDWV);
1001:                              // *
1002:                              // *                 ==== Multiply by U11 ====
1003:                              // *
1004:                              this._dgemm.Run("N", "N", JLEN, I2, J2, ONE
1005:                                              , H, JROW+(INCOL + 1) * LDH + o_h, LDH, U, offset_u, LDU, ONE, ref WV, offset_wv
1006:                                              , LDWV);
1007:                              // *
1008:                              // *                 ==== Copy left of H to right of scratch ====
1009:                              // *
1010:                              this._dlacpy.Run("ALL", JLEN, J2, H, JROW+(INCOL + 1) * LDH + o_h, LDH, ref WV, 1+(1 + I2) * LDWV + o_wv
1011:                                               , LDWV);
1012:                              // *
1013:                              // *                 ==== Multiply by U21 ====
1014:                              // *
1015:                              this._dtrmm.Run("R", "L", "N", "N", JLEN, I4 - I2
1016:                                              , ONE, U, 1+(I2 + 1) * LDU + o_u, LDU, ref WV, 1+(1 + I2) * LDWV + o_wv, LDWV);
1017:                              // *
1018:                              // *                 ==== Multiply by U22 ====
1019:                              // *
1020:                              this._dgemm.Run("N", "N", JLEN, I4 - I2, J4 - J2, ONE
1021:                                              , H, JROW+(INCOL + 1 + J2) * LDH + o_h, LDH, U, J2 + 1+(I2 + 1) * LDU + o_u, LDU, ONE, ref WV, 1+(1 + I2) * LDWV + o_wv
1022:                                              , LDWV);
1023:                              // *
1024:                              // *                 ==== Copy it back ====
1025:                              // *
1026:                              this._dlacpy.Run("ALL", JLEN, KDU, WV, offset_wv, LDWV, ref H, JROW+(INCOL + 1) * LDH + o_h
1027:                                               , LDH);
1028:                          }
1029:                          // *
1030:                          // *              ==== Multiply Z (also vertical) ====
1031:                          // *
1032:                          if (WANTZ)
1033:                          {
1034:                              for (JROW = ILOZ; (NV >= 0) ? (JROW <= IHIZ) : (JROW >= IHIZ); JROW += NV)
1035:                              {
1036:                                  JLEN = Math.Min(NV, IHIZ - JROW + 1);
1037:                                  // *
1038:                                  // *                    ==== Copy right of Z to left of scratch (first
1039:                                  // *                    .     KZS columns get multiplied by zero) ====
1040:                                  // *
1041:                                  this._dlacpy.Run("ALL", JLEN, KNZ, Z, JROW+(INCOL + 1 + J2) * LDZ + o_z, LDZ, ref WV, 1+(1 + KZS) * LDWV + o_wv
1042:                                                   , LDWV);
1043:                                  // *
1044:                                  // *                    ==== Multiply by U12 ====
1045:                                  // *
1046:                                  this._dlaset.Run("ALL", JLEN, KZS, ZERO, ZERO, ref WV, offset_wv
1047:                                                   , LDWV);
1048:                                  this._dtrmm.Run("R", "U", "N", "N", JLEN, KNZ
1049:                                                  , ONE, U, J2 + 1+(1 + KZS) * LDU + o_u, LDU, ref WV, 1+(1 + KZS) * LDWV + o_wv, LDWV);
1050:                                  // *
1051:                                  // *                    ==== Multiply by U11 ====
1052:                                  // *
1053:                                  this._dgemm.Run("N", "N", JLEN, I2, J2, ONE
1054:                                                  , Z, JROW+(INCOL + 1) * LDZ + o_z, LDZ, U, offset_u, LDU, ONE, ref WV, offset_wv
1055:                                                  , LDWV);
1056:                                  // *
1057:                                  // *                    ==== Copy left of Z to right of scratch ====
1058:                                  // *
1059:                                  this._dlacpy.Run("ALL", JLEN, J2, Z, JROW+(INCOL + 1) * LDZ + o_z, LDZ, ref WV, 1+(1 + I2) * LDWV + o_wv
1060:                                                   , LDWV);
1061:                                  // *
1062:                                  // *                    ==== Multiply by U21 ====
1063:                                  // *
1064:                                  this._dtrmm.Run("R", "L", "N", "N", JLEN, I4 - I2
1065:                                                  , ONE, U, 1+(I2 + 1) * LDU + o_u, LDU, ref WV, 1+(1 + I2) * LDWV + o_wv, LDWV);
1066:                                  // *
1067:                                  // *                    ==== Multiply by U22 ====
1068:                                  // *
1069:                                  this._dgemm.Run("N", "N", JLEN, I4 - I2, J4 - J2, ONE
1070:                                                  , Z, JROW+(INCOL + 1 + J2) * LDZ + o_z, LDZ, U, J2 + 1+(I2 + 1) * LDU + o_u, LDU, ONE, ref WV, 1+(1 + I2) * LDWV + o_wv
1071:                                                  , LDWV);
1072:                                  // *
1073:                                  // *                    ==== Copy the result back to Z ====
1074:                                  // *
1075:                                  this._dlacpy.Run("ALL", JLEN, KDU, WV, offset_wv, LDWV, ref Z, JROW+(INCOL + 1) * LDZ + o_z
1076:                                                   , LDZ);
1077:                              }
1078:                          }
1079:                      }
1080:                  }
1081:              }
1082:              // *
1083:              // *     ==== End of DLAQR5 ====
1084:              // *
1085:   
1086:              #endregion
1087:   
1088:          }
1089:      }
1090:  }