Accumulate the stabilized elementary similarity transformations used
in the reduction of a real nonsymmetric matrix to upper Hessenberg form by elmhes.
545209
private static void eltran(DenseMatrix A, DenseMatrix V, int[] perm) {
        int n = A.nrows();
        for (int mp = n - 2; mp > 0; mp--) {
            for (int k = mp + 1; k < n; k++) {
                V.set(k, mp, A.get(k, mp - 1));
            }
            int i = perm[mp];
            if (i != mp) {
                for (int j = mp; j < n; j++) {
                    V.set(mp, j, V.get(i, j));
                    V.set(i, j, 0.0);
                }
                V.set(i, mp, 1.0);
            }
        }
    }
541691
def cublasStrsm(handle, side, uplo, trans, diag, m, n, alpha, A, lda, B, ldb):
    
    
    status = _libcublas.cublasStrsm_v2(handle, 
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       m, n, ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda, int(B), ldb)
    cublasCheckStatus(status)
351348
function lu_custom(A) {
    // Outer product LU decomposition with partial pivoting
    // Ref: Algo 3.4.1 Golub and Van Loan
    if (A.shape.length != 2) {
        throw new Error('Input is not a Matrix (2D)');
    }
    if (A.shape[0] !== A.shape[1]) {
        throw new Error('Input is not a Square Matrix');
    }
    var n = A.shape[0];
    var perm = new common_1.NDArray({ shape: [n] });
    for (var i = 0; i < n; i++) {
        perm.set(i, i);
    }
    function recordPermutation(ri, rj) {
        var tmp = perm.get(ri);
        perm.set(ri, perm.get(rj));
        perm.set(rj, tmp);
    }
    for (var k = 0; k < n - 1; k++) {
        // Find the maximum absolute entry in k'th column
        var ipivot = 0;
        var pivot = -Infinity;
        for (var i = k; i < n; i++) {
            var val = Math.abs(A.get(i, k)); // TODO:Complex
            if (val > pivot) {
                pivot = val;
                ipivot = i;
            }
        }
        // Swap rows k and ipivot
        A.swaprows(k, ipivot);
        recordPermutation(k, ipivot);
        if (common_1.iszero(pivot)) {
            throw new Error('Can\'t perform LU decomp. 0 on diagonal');
        }
        for (var i = k + 1; i < n; i++) {
            A.set(i, k, A.get(i, k) / pivot); // TODO:Complex
        }
        for (var i = k + 1; i < n; i++) {
            for (var j = n - 1; j > k; j--) {
                // TODO:Complex
                A.set(i, j, A.get(i, j) - A.get(i, k) * A.get(k, j));
            }
        }
    }
    return perm;
}
507755
def advance(self):
        """"""
        if self.iter >= self.maxiter:
            raise ArgumentError('Maximum number of iterations reached.')
        if self.invariant:
            raise ArgumentError('Krylov subspace was found to be invariant '
                                'in the previous iteration.')

        N = self.V.shape[0]
        k = self.iter

        # the matrix-vector multiplication
        Av = self.A * self.V[:, [k]]

        if self.ortho == 'house':
            # Householder
            for j in range(k+1):
                Av[j:] = self.houses[j].apply(Av[j:])
                Av[j] *= numpy.conj(self.houses[j].alpha)
            if k+1 < N:
                house = House(Av[k+1:])
                self.houses.append(house)
                Av[k+1:] = house.apply(Av[k+1:]) * numpy.conj(house.alpha)
                self.H[:k+2, [k]] = Av[:k+2]
            else:
                self.H[:k+1, [k]] = Av[:k+1]
            # next line is safe due to the multiplications with alpha
            self.H[k+1, k] = numpy.abs(self.H[k+1, k])
            if self.H[k+1, k] / numpy.linalg.norm(self.H[:k+2, :k+1], 2)\
                    <= 1e-14:
                self.invariant = True
            else:
                vnew = numpy.zeros((N, 1), dtype=self.dtype)
                vnew[k+1] = 1
                for j in range(k+1, -1, -1):
                    vnew[j:] = self.houses[j].apply(vnew[j:])
                self.V[:, [k+1]] = vnew * self.houses[-1].alpha
        else:
            # determine vectors for orthogonalization
            start = 0

            # Lanczos?
            if self.ortho == 'lanczos':
                start = k
                if k > 0:
                    self.H[k-1, k] = self.H[k, k-1]
                    if self.M is not None \
                            and not isinstance(self.M, IdentityLinearOperator):
                        Av -= self.H[k, k-1] * self.P[:, [k-1]]
                    else:
                        Av -= self.H[k, k-1] * self.V[:, [k-1]]

            # (double) modified Gram-Schmidt
            for reortho in range(self.reorthos+1):
                # orthogonalize
                for j in range(start, k+1):
                    alpha = inner(self.V[:, [j]], Av, ip_B=self.ip_B)[0, 0]
                    if self.ortho == 'lanczos':
                        # check if alpha is real
                        if abs(alpha.imag) > 1e-10:
                            warnings.warn(
                                'Iter {0}: abs(alpha.imag) = {1} > 1e-10. '
                                'Is your operator self-adjoint in the '
                                'provided inner product?'
                                .format(self.iter, abs(alpha.imag)))
                        alpha = alpha.real
                    self.H[j, k] += alpha
                    if self.M is not None:
                        Av -= alpha * self.P[:, [j]]
                    else:
                        Av -= alpha * self.V[:, [j]]
            if self.M is not None:
                MAv = self.M * Av
                self.H[k+1, k] = norm(Av, MAv, ip_B=self.ip_B)
            else:
                self.H[k+1, k] = norm(Av, ip_B=self.ip_B)
            if self.H[k+1, k] / numpy.linalg.norm(self.H[:k+2, :k+1], 2)\
                    <= 1e-14:
                self.invariant = True
            else:
                if self.M is not None:
                    self.P[:, [k+1]] = Av / self.H[k+1, k]
                    self.V[:, [k+1]] = MAv / self.H[k+1, k]
                else:
                    self.V[:, [k+1]] = Av / self.H[k+1, k]

        # increase iteration counter
        self.iter += 1
860019
public void decompose( DMatrixRMaj A ) {

        this.QR = A.copy();

        int N = Math.min(A.numCols,A.numRows);

        gammas = new double[ A.numCols ];

        DMatrixRMaj A_small = new DMatrixRMaj(A.numRows,A.numCols);
        DMatrixRMaj A_mod = new DMatrixRMaj(A.numRows,A.numCols);
        DMatrixRMaj v = new DMatrixRMaj(A.numRows,1);
        DMatrixRMaj Q_k = new DMatrixRMaj(A.numRows,A.numRows);

        for( int i = 0; i < N; i++ ) {
            // reshape temporary variables
            A_small.reshape(QR.numRows-i,QR.numCols-i);
            v.reshape(A_small.numRows,1);
            Q_k.reshape(v.getNumElements(),v.getNumElements());

            // use extract matrix to get the column that is to be zeroed
            CommonOps_DDRM.extract(QR,i,QR.numRows,i,i+1,v,0,0);

            double max = CommonOps_DDRM.elementMaxAbs(v);

            if( max > 0 && v.getNumElements() > 1 ) {
                // normalize to reduce overflow issues
                CommonOps_DDRM.divide(v,max);

                // compute the magnitude of the vector
                double tau = NormOps_DDRM.normF(v);

                if( v.get(0) < 0 )
                    tau *= -1.0;

                double u_0 = v.get(0) + tau;
                double gamma = u_0/tau;

                CommonOps_DDRM.divide(v,u_0);
                v.set(0,1.0);

                // extract the submatrix of A which is being operated on
                CommonOps_DDRM.extract(QR,i,QR.numRows,i,QR.numCols,A_small,0,0);

                // A = (I - &gamma;*u*u<sup>T</sup>)A
                CommonOps_DDRM.setIdentity(Q_k);
                CommonOps_DDRM.multAddTransB(-gamma,v,v,Q_k);
                CommonOps_DDRM.mult(Q_k,A_small,A_mod);

                // save the results
                CommonOps_DDRM.insert(A_mod, QR, i,i);
                CommonOps_DDRM.insert(v, QR, i,i);
                QR.unsafe_set(i,i,-tau*max);

                // save gamma for recomputing Q later on
                gammas[i] = gamma;
            }
        }
    }
509502
public static ComplexFloatMatrix gerc(ComplexFloat alpha, ComplexFloatMatrix x,
			ComplexFloatMatrix y, ComplexFloatMatrix a) {
		NativeBlas.cgerc(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
				0, a.rows);
		return a;
	}
619428
private void hqr2 (double[][] V, double[][] H) {
   
      //  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 = Math.pow(2.0,-52.0);
      double exshift = 0.0;
      double p=0,q=0,r=0,s=0,z=0,t,w,x,y;
   
      // Store roots isolated by balanc and compute matrix norm
   
      double norm = 0.0;
      for (int i = 0; i < nn; i++) {
         if (i < low | i > high) {
            d[i] = H[i][i];
            e[i] = 0.0;
         }
         for (int j = Math.max(i-1,0); j < nn; j++) {
            norm = norm + Math.abs(H[i][j]);
         }
      }
   
      // 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 = Math.abs(H[l-1][l-1]) + Math.abs(H[l][l]);
            if (s == 0.0) {
               s = norm;
            }
            if (Math.abs(H[l][l-1]) < eps * s) {
               break;
            }
            l--;
         }
       
         // Check for convergence
         // One root found
   
         if (l == n) {
            H[n][n] = H[n][n] + exshift;
            d[n] = H[n][n];
            e[n] = 0.0;
            n--;
            iter = 0;
   
         // Two roots found
   
         } else if (l == n-1) {
            w = H[n][n-1] * H[n-1][n];
            p = (H[n-1][n-1] - H[n][n]) / 2.0;
            q = p * p + w;
            z = Math.sqrt(Math.abs(q));
            H[n][n] = H[n][n] + exshift;
            H[n-1][n-1] = H[n-1][n-1] + exshift;
            x = H[n][n];
   
            // Real pair
   
            if (q >= 0) {
               if (p >= 0) {
                  z = p + z;
               } else {
                  z = p - z;
               }
               d[n-1] = x + z;
               d[n] = d[n-1];
               if (z != 0.0) {
                  d[n] = x - w / z;
               }
               e[n-1] = 0.0;
               e[n] = 0.0;
               x = H[n][n-1];
               s = Math.abs(x) + Math.abs(z);
               p = x / s;
               q = z / s;
               r = Math.sqrt(p * p+q * q);
               p = p / r;
               q = q / r;
   
               // Row modification
   
               for (int j = n-1; j < nn; j++) {
                  z = H[n-1][j];
                  H[n-1][j] = q * z + p * H[n][j];
                  H[n][j] = q * H[n][j] - p * z;
               }
   
               // Column modification
   
               for (int i = 0; i <= n; i++) {
                  z = H[i][n-1];
                  H[i][n-1] = q * z + p * H[i][n];
                  H[i][n] = q * H[i][n] - p * z;
               }
   
               // Accumulate transformations
   
               for (int i = low; i <= high; i++) {
                  z = V[i][n-1];
                  V[i][n-1] = q * z + p * V[i][n];
                  V[i][n] = q * V[i][n] - p * z;
               }
   
            // Complex pair
   
            } else {
               d[n-1] = x + p;
               d[n] = x + p;
               e[n-1] = z;
               e[n] = -z;
            }
            n = n - 2;
            iter = 0;
   
         // No convergence yet
   
         } else {
   
            // Form shift
   
            x = H[n][n];
            y = 0.0;
            w = 0.0;
            if (l < n) {
               y = H[n-1][n-1];
               w = H[n][n-1] * H[n-1][n];
            }
   
            // Wilkinson's original ad hoc shift
   
            if (iter == 10) {
               exshift += x;
               for (int i = low; i <= n; i++) {
                  H[i][i] -= x;
               }
               s = Math.abs(H[n][n-1]) + Math.abs(H[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 = Math.sqrt(s);
                    if (y < x) {
                       s = -s;
                    }
                    s = x - w / ((y - x) / 2.0 + s);
                    for (int i = low; i <= n; i++) {
                       H[i][i] -= 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[m][m];
               r = x - z;
               s = y - z;
               p = (r * s - w) / H[m+1][m] + H[m][m+1];
               q = H[m+1][m+1] - z - r - s;
               r = H[m+2][m+1];
               s = Math.abs(p) + Math.abs(q) + Math.abs(r);
               p = p / s;
               q = q / s;
               r = r / s;
               if (m == l) {
                  break;
               }
               if (Math.abs(H[m][m-1]) * (Math.abs(q) + Math.abs(r)) <
                  eps * (Math.abs(p) * (Math.abs(H[m-1][m-1]) + Math.abs(z) +
                  Math.abs(H[m+1][m+1])))) {
                     break;
               }
               m--;
            }
   
            for (int i = m+2; i <= n; i++) {
               H[i][i-2] = 0.0;
               if (i > m+2) {
                  H[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[k][k-1];
                  q = H[k+1][k-1];
                  r = (notlast ? H[k+2][k-1] : 0.0);
                  x = Math.abs(p) + Math.abs(q) + Math.abs(r);
                  if (x != 0.0) {
                     p = p / x;
                     q = q / x;
                     r = r / x;
                  }
               }
               if (x == 0.0) {
                  break;
               }
               s = Math.sqrt(p * p + q * q + r * r);
               if (p < 0) {
                  s = -s;
               }
               if (s != 0) {
                  if (k != m) {
                     H[k][k-1] = -s * x;
                  } else if (l != m) {
                     H[k][k-1] = -H[k][k-1];
                  }
                  p = p + s;
                  x = p / s;
                  y = q / s;
                  z = r / s;
                  q = q / p;
                  r = r / p;
   
                  // Row modification
   
                  for (int j = k; j < nn; j++) {
                     p = H[k][j] + q * H[k+1][j];
                     if (notlast) {
                        p = p + r * H[k+2][j];
                        H[k+2][j] = H[k+2][j] - p * z;
                     }
                     H[k][j] = H[k][j] - p * x;
                     H[k+1][j] = H[k+1][j] - p * y;
                  }
   
                  // Column modification
   
                  for (int i = 0; i <= Math.min(n,k+3); i++) {
                     p = x * H[i][k] + y * H[i][k+1];
                     if (notlast) {
                        p = p + z * H[i][k+2];
                        H[i][k+2] = H[i][k+2] - p * r;
                     }
                     H[i][k] = H[i][k] - p;
                     H[i][k+1] = H[i][k+1] - p * q;
                  }
   
                  // Accumulate transformations
   
                  for (int i = low; i <= high; i++) {
                     p = x * V[i][k] + y * V[i][k+1];
                     if (notlast) {
                        p = p + z * V[i][k+2];
                        V[i][k+2] = V[i][k+2] - p * r;
                     }
                     V[i][k] = V[i][k] - p;
                     V[i][k+1] = V[i][k+1] - p * q;
                  }
               }  // (s != 0)
            }  // k loop
         }  // check convergence
      }  // while (n >= low)
      
      // Backsubstitute to find vectors of upper triangular form

      if (norm == 0.0) {
         return;
      }
   
      for (n = nn-1; n >= 0; n--) {
         p = d[n];
         q = e[n];
   
         // Real vector
   
         if (q == 0) {
            int l = n;
            H[n][n] = 1.0;
            for (int i = n-1; i >= 0; i--) {
               w = H[i][i] - p;
               r = 0.0;
               for (int j = l; j <= n; j++) {
                  r = r + H[i][j] * H[j][n];
               }
               if (e[i] < 0.0) {
                  z = w;
                  s = r;
               } else {
                  l = i;
                  if (e[i] == 0.0) {
                     if (w != 0.0) {
                        H[i][n] = -r / w;
                     } else {
                        H[i][n] = -r / (eps * norm);
                     }
   
                  // Solve real equations
   
                  } else {
                     x = H[i][i+1];
                     y = H[i+1][i];
                     q = (d[i] - p) * (d[i] - p) + e[i] * e[i];
                     t = (x * s - z * r) / q;
                     H[i][n] = t;
                     if (Math.abs(x) > Math.abs(z)) {
                        H[i+1][n] = (-r - w * t) / x;
                     } else {
                        H[i+1][n] = (-s - y * t) / z;
                     }
                  }
   
                  // Overflow control
   
                  t = Math.abs(H[i][n]);
                  if ((eps * t) * t > 1) {
                     for (int j = i; j <= n; j++) {
                        H[j][n] = H[j][n] / t;
                     }
                  }
               }
            }
   
         // Complex vector
   
         } else if (q < 0) {
            int l = n-1;

            // Last vector component imaginary so matrix is triangular
   
            if (Math.abs(H[n][n-1]) > Math.abs(H[n-1][n])) {
               H[n-1][n-1] = q / H[n][n-1];
               H[n-1][n] = -(H[n][n] - p) / H[n][n-1];
            } else {
               cdiv(0.0,-H[n-1][n],H[n-1][n-1]-p,q);
               H[n-1][n-1] = cdivr;
               H[n-1][n] = cdivi;
            }
            H[n][n-1] = 0.0;
            H[n][n] = 1.0;
            for (int i = n-2; i >= 0; i--) {
               double ra,sa,vr,vi;
               ra = 0.0;
               sa = 0.0;
               for (int j = l; j <= n; j++) {
                  ra = ra + H[i][j] * H[j][n-1];
                  sa = sa + H[i][j] * H[j][n];
               }
               w = H[i][i] - p;
   
               if (e[i] < 0.0) {
                  z = w;
                  r = ra;
                  s = sa;
               } else {
                  l = i;
                  if (e[i] == 0) {
                     cdiv(-ra,-sa,w,q);
                     H[i][n-1] = cdivr;
                     H[i][n] = cdivi;
                  } else {
   
                     // Solve complex equations
   
                     x = H[i][i+1];
                     y = H[i+1][i];
                     vr = (d[i] - p) * (d[i] - p) + e[i] * e[i] - q * q;
                     vi = (d[i] - p) * 2.0 * q;
                     if (vr == 0.0 & vi == 0.0) {
                        vr = eps * norm * (Math.abs(w) + Math.abs(q) +
                        Math.abs(x) + Math.abs(y) + Math.abs(z));
                     }
                     cdiv(x*r-z*ra+q*sa,x*s-z*sa-q*ra,vr,vi);
                     H[i][n-1] = cdivr;
                     H[i][n] = cdivi;
                     if (Math.abs(x) > (Math.abs(z) + Math.abs(q))) {
                        H[i+1][n-1] = (-ra - w * H[i][n-1] + q * H[i][n]) / x;
                        H[i+1][n] = (-sa - w * H[i][n] - q * H[i][n-1]) / x;
                     } else {
                        cdiv(-r-y*H[i][n-1],-s-y*H[i][n],z,q);
                        H[i+1][n-1] = cdivr;
                        H[i+1][n] = cdivi;
                     }
                  }
   
                  // Overflow control

                  t = Math.max(Math.abs(H[i][n-1]),Math.abs(H[i][n]));
                  if ((eps * t) * t > 1) {
                     for (int j = i; j <= n; j++) {
                        H[j][n-1] = H[j][n-1] / t;
                        H[j][n] = H[j][n] / t;
                     }
                  }
               }
            }
         }
      }
   
      // Vectors of isolated roots
   
      for (int i = 0; i < nn; i++) {
         if (i < low | i > high) {
            for (int j = i; j < nn; j++) {
               V[i][j] = H[i][j];
            }
         }
      }
   
      // Back transformation to get eigenvectors of original matrix
   
      for (int j = nn-1; j >= low; j--) {
         for (int i = low; i <= high; i++) {
            z = 0.0;
            for (int k = low; k <= Math.min(j,high); k++) {
               z = z + V[i][k] * H[k][j];
            }
            V[i][j] = z;
         }
      }
   }
314543
public static void update(DMatrixRMaj H ,
							  DMatrixRMaj s ,
							  DMatrixRMaj y ,
							  DMatrixRMaj tempV0, DMatrixRMaj tempV1) {
		double p = VectorVectorMult_DDRM.innerProd(y,s);
		if( p == 0 )
			return;

		p = 1.0/p;

		double sBs = VectorVectorMult_DDRM.innerProdA(s,H,s);
		if( sBs == 0 )
			return;

		CommonOps_DDRM.mult(H,s,tempV0);
		CommonOps_DDRM.multTransA(s,H,tempV1);

		VectorVectorMult_DDRM.rank1Update(-p, H , tempV0, y);
		VectorVectorMult_DDRM.rank1Update(-p, H , y, tempV1);
		VectorVectorMult_DDRM.rank1Update(p*(p*sBs+1), H , y, y);
	}
141015
function eig(A) {
    if (A.shape.length !== 2) {
        throw new Error('Input is not matrix');
    }
    if (A.shape[0] !== A.shape[1]) {
        throw new Error('Input is not square matrix');
    }
    var n = A.shape[0];
    var copyA = A.clone();
    copyA.swapOrder();
    var _a = lapack.geev(copyA.data, n, true, true), WR = _a[0], WI = _a[1], VL = _a[2], VR = _a[3];
    var eigval = new common_1.NDArray({ shape: [n] });
    var eigvecL = new common_1.NDArray({ shape: [n, n] });
    var eigvecR = new common_1.NDArray({ shape: [n, n] });
    for (var i = 0; i < n; i++) {
        if (common_1.iszero(WI[i])) {
            eigval.set(i, WR[i]);
        }
        else {
            eigval.set(i, new common_1.Complex(WR[i], WI[i]));
        }
    }
    for (var i = 0; i < n; i++) {
        for (var j = 0; j < n;) {
            if (common_1.iszero(WI[j])) {
                // We want to extract eigen-vectors in i'th column of VL and VR
                eigvecL.set(i, j, VL[j * n + i]);
                eigvecR.set(i, j, VR[j * n + i]);
                j += 1;
            }
            else {
                // i-th eigen vector will be complex and
                // i+1-th eigen vector will be its conjugate
                // There are n eigen-vectors for i'th eigen-value
                eigvecL.set(i, j, new common_1.Complex(VL[j * n + i], VL[(j + 1) * n + i]));
                eigvecL.set(i, j + 1, new common_1.Complex(VL[j * n + i], -VL[(j + 1) * n + i]));
                eigvecR.set(i, j, new common_1.Complex(VR[j * n + i], VR[(j + 1) * n + i]));
                eigvecR.set(i, j + 1, new common_1.Complex(VR[j * n + i], -VR[(j + 1) * n + i]));
                j += 2;
            }
        }
    }
    return [eigval, eigvecL, eigvecR];
}
234474
private function tred2(): void
    {
        //  This is derived from the Algol procedures tred2 by
        //  Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
        //  Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
        //  Fortran subroutine in EISPACK.
        $this->d = $this->V[$this->n - 1];
        // Householder reduction to tridiagonal form.
        for ($i = $this->n - 1; $i > 0; --$i) {
            $i_ = $i - 1;
            // Scale to avoid under/overflow.
            $h = $scale = 0.0;
            $scale += array_sum(array_map('abs', $this->d));
            if ($scale == 0.0) {
                $this->e[$i] = $this->d[$i_];
                $this->d = array_slice($this->V[$i_], 0, $this->n - 1);
                for ($j = 0; $j < $i; ++$j) {
                    $this->V[$j][$i] = $this->V[$i][$j] = 0.0;
                }
            } else {
                // Generate Householder vector.
                for ($k = 0; $k < $i; ++$k) {
                    $this->d[$k] /= $scale;
                    $h += $this->d[$k] ** 2;
                }

                $f = $this->d[$i_];
                $g = $h ** .5;
                if ($f > 0) {
                    $g = -$g;
                }

                $this->e[$i] = $scale * $g;
                $h -= $f * $g;
                $this->d[$i_] = $f - $g;

                for ($j = 0; $j < $i; ++$j) {
                    $this->e[$j] = 0.0;
                }

                // Apply similarity transformation to remaining columns.
                for ($j = 0; $j < $i; ++$j) {
                    $f = $this->d[$j];
                    $this->V[$j][$i] = $f;
                    $g = $this->e[$j] + $this->V[$j][$j] * $f;

                    for ($k = $j + 1; $k <= $i_; ++$k) {
                        $g += $this->V[$k][$j] * $this->d[$k];
                        $this->e[$k] += $this->V[$k][$j] * $f;
                    }

                    $this->e[$j] = $g;
                }

                $f = 0.0;

                if ($h == 0.0) {
                    $h = 1e-32;
                }

                for ($j = 0; $j < $i; ++$j) {
                    $this->e[$j] /= $h;
                    $f += $this->e[$j] * $this->d[$j];
                }

                $hh = $f / (2 * $h);
                for ($j = 0; $j < $i; ++$j) {
                    $this->e[$j] -= $hh * $this->d[$j];
                }

                for ($j = 0; $j < $i; ++$j) {
                    $f = $this->d[$j];
                    $g = $this->e[$j];
                    for ($k = $j; $k <= $i_; ++$k) {
                        $this->V[$k][$j] -= ($f * $this->e[$k] + $g * $this->d[$k]);
                    }

                    $this->d[$j] = $this->V[$i - 1][$j];
                    $this->V[$i][$j] = 0.0;
                }
            }

            $this->d[$i] = $h;
        }

        // Accumulate transformations.
        for ($i = 0; $i < $this->n - 1; ++$i) {
            $this->V[$this->n - 1][$i] = $this->V[$i][$i];
            $this->V[$i][$i] = 1.0;
            $h = $this->d[$i + 1];
            if ($h != 0.0) {
                for ($k = 0; $k <= $i; ++$k) {
                    $this->d[$k] = $this->V[$k][$i + 1] / $h;
                }

                for ($j = 0; $j <= $i; ++$j) {
                    $g = 0.0;
                    for ($k = 0; $k <= $i; ++$k) {
                        $g += $this->V[$k][$i + 1] * $this->V[$k][$j];
                    }

                    for ($k = 0; $k <= $i; ++$k) {
                        $this->V[$k][$j] -= $g * $this->d[$k];
                    }
                }
            }

            for ($k = 0; $k <= $i; ++$k) {
                $this->V[$k][$i + 1] = 0.0;
            }
        }

        $this->d = $this->V[$this->n - 1];
        $this->V[$this->n - 1] = array_fill(0, $this->n, 0.0);
        $this->V[$this->n - 1][$this->n - 1] = 1.0;
        $this->e[0] = 0.0;
    }
1224290
def _augmenting_row_reduction(self):
        
        unassigned = np.where(self._x == -1)[0]
        for i in unassigned:
            for _ in range(self.c.size):
                # Time in this loop can be proportional to 1/epsilon
                # This step is not strictly necessary, so cutoff early
                # to avoid near-infinite loops

                # find smallest 2 values and indices
                temp = self.c[i] - self._v
                j1 = np.argmin(temp)
                u1 = temp[j1]
                temp[j1] = np.inf
                j2 = np.argmin(temp)
                u2 = temp[j2]

                if u1 < u2:
                    self._v[j1] -= u2 - u1
                elif self._y[j1] != -1:
                    j1 = j2
                k = self._y[j1]
                if k != -1:
                    self._x[k] = -1
                    self._x[i] = j1
                    self._y[j1] = i
                    i = k
                if k == -1 or abs(u1 - u2) < self.epsilon:
                    break
66810
def apply_M(self, ax, ay):
        
        jac = numpy.array(
            [[self.dx.dot(ax), self.dy.dot(ax)], [self.dx.dot(ay), self.dy.dot(ay)]]
        )

        # jacs and J are of shape (2, 2, k). M must be of the same shape and
        # contain the result of the k 2x2 dot products. Perhaps there's a
        # dot() for this.
        M = numpy.einsum("ijl,jkl->ikl", jac, self.J)
        # M = numpy.array([
        #     [
        #         jac[0][0]*self.J[0][0] + jac[0][1]*self.J[1][0],
        #         jac[0][0]*self.J[0][1] + jac[0][1]*self.J[1][1],
        #     ],
        #     [
        #         jac[1][0]*self.J[0][0] + jac[1][1]*self.J[1][0],
        #         jac[1][0]*self.J[0][1] + jac[1][1]*self.J[1][1],
        #     ],
        #     ])

        # One could use
        #
        #     M = numpy.moveaxis(M, -1, 0)
        #     _, sigma, _ = numpy.linalg.svd(M)
        #
        # but computing the singular values explicitly via
        # <https://scicomp.stackexchange.com/a/14103/3980> is faster and more
        # explicit.
        a = (M[0, 0] + M[1, 1]) / 2
        b = (M[0, 0] - M[1, 1]) / 2
        c = (M[1, 0] + M[0, 1]) / 2
        d = (M[1, 0] - M[0, 1]) / 2

        return a, b, c, d
841948
public Se3_F64 decompose( DMatrixRMaj H )
	{
		// step through each calibration grid and compute its parameters
		DMatrixRMaj h[] = SpecializedOps_DDRM.splitIntoVectors(H, true);

		// lambda = 1/norm(inv(K)*h1) or 1/norm(inv(K)*h2)
		// use the average to attempt to reduce error
		CommonOps_DDRM.mult(K_inv,h[0],temp);
		double lambda = NormOps_DDRM.normF(temp);
		CommonOps_DDRM.mult(K_inv,h[1],temp);
		lambda += NormOps_DDRM.normF(temp);
		lambda = 2.0/lambda;

		// compute the column in the rotation matrix
		CommonOps_DDRM.mult(lambda,K_inv,h[0],r1);
		CommonOps_DDRM.mult(lambda,K_inv,h[1],r2);
		CommonOps_DDRM.mult(lambda,K_inv,h[2],t);

		Vector3D_F64 v1 = UtilVector3D_F64.convert(r1);
		Vector3D_F64 v2 = UtilVector3D_F64.convert(r2);
		Vector3D_F64 v3 = v1.cross(v2);

		UtilVector3D_F64.createMatrix(R, v1, v2, v3);

		Se3_F64 ret = new Se3_F64();
		// the R matrix is probably not a real rotation matrix.  So find
		// the closest real rotation matrix
		ConvertRotation3D_F64.approximateRotationMatrix(R,ret.getR());
		ret.getT().set(t.data[0],t.data[1],t.data[2]);

		return ret;
	}
423540
function hm(n, v, incv, offv, tau, b, offb) {
    // Apply householder transformation b = b - tau v v' b
    // LAPACK: DLARF

    if (tau === 0.0) {
      return;
    }

    // compute v' b
    var d = b[offb + 0];
    for (var i = 1; i < n; ++i) {
      d += b[offb + i] * v[offv + incv*i];
    }

    // compute b = b - tau v v' b
    b[offb + 0] -= tau * d;
    for (var i = 1; i < n; ++i) {
      b[offb + i] -= tau * v[offv + incv*i] * d;
    }
  }
421721
public function gaussianElimination() : array
    {
        $minDim = min($this->shape());

        $b = $this->a;

        $swaps = 0;

        for ($k = 0; $k < $minDim; $k++) {
            $index = $k;

            for ($i = $k; $i < $this->m; $i++) {
                if (abs($b[$i][$k]) > abs($b[$index][$k])) {
                    $index = $i;
                }
            }

            if ($b[$index][$k] == 0) {
                throw new RuntimeException('Cannot compute row echelon form'
                    . ' of a singular matrix.');
            }

            if ($k !== $index) {
                $temp = $b[$k];

                $b[$k] = $b[$index];
                $b[$index] = $temp;

                $swaps++;
            }

            $diag = $b[$k][$k];

            for ($i = $k + 1; $i < $this->m; $i++) {
                $scale = $diag != 0 ? $b[$i][$k] / $diag : 1;
                
                for ($j = $k + 1; $j < $this->n; $j++) {
                    $b[$i][$j] -= $scale * $b[$k][$j];
                }
                
                $b[$i][$k] = 0;
            }
        }

        $b = self::quick($b);

        return [$b, $swaps];
    }
1304808
private SimpleMatrix computeAffineH( List<AssociatedPair> observations ,
										 DMatrixRMaj H , DMatrixRMaj Hzero ) {

		SimpleMatrix A = new SimpleMatrix(observations.size(),3);
		SimpleMatrix b = new SimpleMatrix(A.numRows(),1);

		Point2D_F64 c = new Point2D_F64();
		Point2D_F64 k = new Point2D_F64();

		for( int i = 0; i < observations.size(); i++ ) {
			AssociatedPair a = observations.get(i);

			GeometryMath_F64.mult(Hzero, a.p1, k);
			GeometryMath_F64.mult(H,a.p2,c);

			A.setRow(i,0,k.x,k.y,1);
			b.set(i,0,c.x);
		}

		SimpleMatrix x = A.solve(b);

		SimpleMatrix Ha = SimpleMatrix.identity(3);
		Ha.setRow(0,0,x.getDDRM().data);

		return Ha;
	}
854401
def _graph_reduction(adj, x, g, f):
    
    as_list = set()
    as_nodes = {v for v in adj if len(adj[v]) <= f and is_almost_simplicial(adj, v)}
    while as_nodes:
        as_list.union(as_nodes)
        for n in as_nodes:

            # update g and f
            dv = len(adj[n])
            if dv > g:
                g = dv
            if g > f:
                f = g

            # eliminate v
            x.append(n)
            _elim_adj(adj, n)

        # see if we have any more simplicial nodes
        as_nodes = {v for v in adj if len(adj[v]) <= f and is_almost_simplicial(adj, v)}

    return g, f, as_list