docs: private void hqr2()
      {
  
          //  This is derived from the Algol procedure hqr2,
          //  by Martin and Wilkinson, Handbook for Auto. Comp.,
          //  Vol.ii-Linear Algebra, and the corresponding
          //  Fortran subroutine in EISPACK.
  
          // Initialize
  
          int nn = this.n;
          int n = nn - 1;
          int low = 0;
          int high = nn - 1;
          double eps = pow(2.0, -52.0);
          double exshift = 0.0;
          double p = 0, q = 0, r = 0, s = 0, z = 0, t, w, x, y;
          /**
           * Output from complex division
           */
          final double[] cr = new double[2];
          double norm = hqr2GetNormStart(nn, low, high);
  
          // Outer loop over eigenvalue index
  
          int iter = 0;
          while (n >= low)
          {
  
              // Look for single small sub-diagonal element
  
              int l = n;
              while (l > low)
              {
                  s = abs(H.get(l-1, l-1)) + abs(H.get(l, l));
                  if (s == 0.0)
                  {
                      s = norm;
                  }
                  if (abs(H.get(l, l-1)) < eps * s)
                  {
                      break;
                  }
                  l--;
              }
  
              // Check for convergence
              // One root found
  
              if (l == n)
              {
                  H.increment(n, n, exshift);
                  d[n] = H.get(n, n);
                  e[n] = 0.0;
                  n--;
                  iter = 0;
              }
              else if (l == n - 1) // Two roots found
              {
                  hqr2FoundTwoRoots(exshift, n, nn, low, high);
                  n = n - 2;
                  iter = 0;
  
                  // No convergence yet
              }
              else
              {
  
                  // Form shift
                  
                  x = H.get(n, n);
                  y = 0.0;
                  w = 0.0;
                  if (l < n)
                  {
                      y = H.get(n-1, n-1);
                      w = pow(H.get(n, n-1), 2);
                  }
  
                  // Wilkinson's original ad hoc shift
  
                  if (iter == 10)
                  {
                      exshift += x;
                      RowColumnOps.addDiag(H, low, n+1, -x);
                      s = abs(H.get(n, n-1)) + abs(H.get(n-1, n-2));
                      x = y = 0.75 * s;
                      w = -0.4375 * s * s;
                  }
  
                  // MATLAB's new ad hoc shift
  
                  if (iter == 30)
                  {
                      s = (y - x) / 2.0;
                      s = s * s + w;
                      if (s > 0)
                      {
                          s = sqrt(s);
                          if (y < x)
                          {
                              s = -s;
                          }
                          s = x - w / ((y - x) / 2.0 + s);
                          RowColumnOps.addDiag(H, low, n+1, -s);
                          exshift += s;
                          x = y = w = 0.964;
                      }
                  }
  
                  iter = iter + 1;   // (Could check iteration count here.)
  
                  // Look for two consecutive small sub-diagonal elements
  
                  int m = n - 2;
                  while (m >= l)
                  {
                      z = H.get(m, m);
                      r = x - z;
                      s = y - z;
                      p = (r * s - w) / H.get(m+1, m) + H.get(m, m+1);
                      q = H.get(m+1, m+1) - z - r - s;
                      r = H.get(m+2, m+1);
                      s = abs(p) + abs(q) + abs(r);
                      p = p / s;
                      q = q / s;
                      r = r / s;
                      if (m == l)
                      {
                          break;
                      }
                      if (abs(H.get(m, m-1)) * (abs(q) + abs(r))
                              < eps * (abs(p) * (abs(H.get(m-1, m-1)) + abs(z)
                              + abs(H.get(m+1, m+1)))))
                      {
                          break;
                      }
                      m--;
                  }
  
                  for (int i = m + 2; i <= n; i++)
                  {
                      H.set(i, i-2, 0.0);
                      if (i > m + 2)
                      {
                          H.set(i, i-3, 0.0);
                      }
                  }
  
                  // Double QR step involving rows l:n and columns m:n
  
                  for (int k = m; k <= n - 1; k++)
                  {
                      boolean notlast = (k != n - 1);
                      if (k != m)
                      {
                          p = H.get(k, k-1);
                          q = H.get(k+1, k-1);
                          r = (notlast ? H.get(k+2, k-1) : 0.0);
                          x = abs(p) + abs(q) + abs(r);
                          if (x != 0.0)
                          {
                              p = p / x;
                              q = q / x;
                              r = r / x;
                          }
                      }
                      
                      if (x == 0.0)
                          break;
                      
                      s = sqrt(p * p + q * q + r * r);
                      if (p < 0)
                      {
                          s = -s;
                      }
                      if (s != 0)
                      {
                          if (k != m)
                          {
                              H.set(k, k-1, -s*x);
                          }
                          else if (l != m)
                          {
                              H.set(k, k-1, -H.get(k, k-1));
                          }
                          p = p + s;
                          x = p / s;
                          y = q / s;
                          z = r / s;
                          q = q / p;
                          r = r / p;
  
                          // Row modification
                          rowOpTransform2(H, k, nn - 1, x, k, y, notlast, z, r, q);
  
                          // Column modification
                          columnOpTransform2(H, 0, min(n, k + 3), x, k, y, notlast, z, r, q);
  
                          // Accumulate transformations
                          columnOpTransform2(V, low, high, x, k, y, notlast, z, r, q);
                          
                      }  // (s != 0)
                  }  // k loop
              }  // check convergence
          }  // while (n >= low)
          
          // Backsubstitute to find vectors of upper triangular form
  
          if (norm == 0.0)
              return;
          
          backsubtituteFindVectors(nn, z, s, eps, norm, cr);
          
          // Vectors of isolated roots
  
          for (int i = 0; i < nn; i++)
              if (i < low | i > high)
              {
                  for(int j = i; j < nn-1; j++)
                      H.set(i, j, V.get(i, j));
              }
          backtransform(nn, low, high);
      }