docs: private void orthes (double[][] V, double[][] H) {
  		double[] ort=new double[n];
  		
        //  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.
     
        int low = 0;
        int 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 = scale + Math.abs(H[i][m-1]);
           }
           if (scale != 0.0) {
     
              // Compute Householder transformation.
     
              double h = 0.0;
              for (int i = high; i >= m; i--) {
                 ort[i] = H[i][m-1]/scale;
                 h += ort[i] * ort[i];
              }
              double g = Math.sqrt(h);
              if (ort[m] > 0) {
                 g = -g;
              }
              h = h - ort[m] * g;
              ort[m] = ort[m] - 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 = f/h;
                 for (int i = m; i <= high; i++) {
                    H[i][j] -= f*ort[i];
                 }
             }
     
             for (int i = 0; i <= high; i++) {
                 double f = 0.0;
                 for (int j = high; j >= m; j--) {
                    f += ort[j]*H[i][j];
                 }
                 f = f/h;
                 for (int j = m; j <= high; j++) {
                    H[i][j] -= f*ort[j];
                 }
              }
              ort[m] = scale*ort[m];
              H[m][m-1] = scale*g;
           }
        }
     
        // Accumulate transformations (Algol's ortran).
  
        for (int i = 0; i < n; i++) {
           for (int j = 0; j < n; j++) {
              V[i][j] = (i == j ? 1.0 : 0.0);
           }
        }
  
        for (int m = high-1; m >= low+1; m--) {
           if (H[m][m-1] != 0.0) {
              for (int i = m+1; i <= high; i++) {
                 ort[i] = H[i][m-1];
              }
              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];
                 }
              }
           }
        }
     }