docs: private void orthes() {
      // FIXME: does this fail on NaN/inf values?
  
      // This is derived from the Algol procedures orthes and ortran,
      // by Martin and Wilkinson, Handbook for Auto. Comp.,
      // Vol.ii-Linear Algebra, and the corresponding
      // Fortran subroutines in EISPACK.
  
      final int low = 0, high = n - 1;
      for(int m = low + 1; m <= high - 1; m++) {
        // Scale column.
        double scale = 0.0;
        for(int i = m; i <= high; i++) {
          scale += abs(H[i][m - 1]);
        }
        if(scale > 0.0 || scale < 0.0) {
          // Compute Householder transformation.
          double h = 0.0;
          for(int i = high; i >= m; i--) {
            double oi = ort[i] = H[i][m - 1] / scale;
            h += oi * oi;
          }
          double g = FastMath.sqrt(h), om = ort[m];
          g = (om > 0) ? -g : g;
          h -= om * g;
          ort[m] = om - g;
  
          // Apply Householder similarity transformation
          // H = (I-u*u'/h)*H*(I-u*u')/h)
  
          for(int j = m; j < n; j++) {
            double f = 0.0;
            for(int i = high; i >= m; i--) {
              f += ort[i] * H[i][j];
            }
            f /= h;
            for(int i = m; i <= high; i++) {
              H[i][j] -= f * ort[i];
            }
          }
  
          for(int i = 0; i <= high; i++) {
            final double[] Hi = H[i];
            double f = 0.0;
            for(int j = high; j >= m; j--) {
              f += ort[j] * Hi[j];
            }
            f /= h;
            for(int j = m; j <= high; j++) {
              Hi[j] -= f * ort[j];
            }
          }
          ort[m] *= scale;
          H[m][m - 1] = scale * g;
        }
      }
  
      // Accumulate transformations (Algol's ortran).
      for(int i = 0; i < n; i++) {
        Arrays.fill(V[i], 0.);
        V[i][i] = 1.;
      }
  
      for(int m = high - 1; m >= low + 1; m--) {
        final double[] H_m = H[m];
        if(H_m[m - 1] != 0.0) {
          for(int i = m + 1; i <= high; i++) {
            ort[i] = H[i][m - 1];
          }
          final double ort_m = ort[m];
          for(int j = m; j <= high; j++) {
            double g = 0.0;
            for(int i = m; i <= high; i++) {
              g += ort[i] * V[i][j];
            }
            // Double division avoids possible underflow
            g = (g / ort_m) / H_m[m - 1];
            for(int i = m; i <= high; i++) {
              V[i][j] += g * ort[i];
            }
          }
        }
      }
    }