docs: private boolean _decompose() {
          double h[] = QH.data;
  
          for( int k = 0; k < N-2; k++ ) {
              // find the largest value in this column
              // this is used to normalize the column and mitigate overflow/underflow
              double max = 0;
  
              for( int i = k+1; i < N; i++ ) {
                  // copy the householder vector to vector outside of the matrix to reduce caching issues
                  // big improvement on larger matrices and a relatively small performance hit on small matrices.
                  double val = u[i] = h[i*N+k];
                  val = Math.abs(val);
                  if( val > max )
                      max = val;
              }
  
              if( max > 0 ) {
                  // -------- set up the reflector Q_k
  
                  double tau = 0;
                  // normalize to reduce overflow/underflow
                  // and compute tau for the reflector
                  for( int i = k+1; i < N; i++ ) {
                      double val = u[i] /= max;
                      tau += val*val;
                  }
  
                  tau = Math.sqrt(tau);
  
                  if( u[k+1] < 0 )
                      tau = -tau;
  
                  // write the reflector into the lower left column of the matrix
                  double nu = u[k+1] + tau;
                  u[k+1] = 1.0;
  
                  for( int i = k+2; i < N; i++ ) {
                      h[i*N+k] = u[i] /= nu;
                  }
  
                  double gamma = nu/tau;
                  gammas[k] = gamma;
  
                  // ---------- multiply on the left by Q_k
                  QrHelperFunctions_DDRM.rank1UpdateMultR(QH, u, gamma, k + 1, k + 1, N, b);
  
                  // ---------- multiply on the right by Q_k
                  QrHelperFunctions_DDRM.rank1UpdateMultL(QH, u, gamma, 0, k + 1, N);
  
                  // since the first element in the householder vector is known to be 1
                  // store the full upper hessenberg
                  h[(k+1)*N+k] = -tau*max;
  
              } else {
                  gammas[k] = 0;
              }
  
          }
  
          return true;
      }