Symmetric tridiagonal QL algorithm.
766162
private void tql2() {
    // This is derived from the Algol procedures tql2, by
    // Bowdler, Martin, Reinsch, and Wilkinson, Handbook for
    // Auto. Comp., Vol.ii-Linear Algebra, and the corresponding
    // Fortran subroutine in EISPACK.

    System.arraycopy(e, 1, e, 0, n - 1);
    e[n - 1] = 0.;

    double f = 0., tst1 = 0.;
    for(int l = 0; l < n; l++) {
      // Find small subdiagonal element
      tst1 = max(tst1, abs(d[l]) + abs(e[l]));
      int m = l;
      for(; m < n; ++m) {
        if(abs(e[m]) <= 0x1P-52 * tst1) {
          break;
        }
      }

      // If m == l, d[l] is an eigenvalue,
      // otherwise, iterate.
      if(m > l) {
        do {
          // Compute implicit shift
          f += tql2ComputeImplicitShift(l);
          // Implicit QL transformation.
          tql2ImplicitQL(l, m, d[l + 1]);
          // Check for convergence.
          // TODO: Iteration limit?
        }
        while(abs(e[l]) > 0x1P-52 * tst1);
      }
      d[l] += f;
      e[l] = 0.;
    }
  }
7858
private void computeSolutions(DMatrix4x4 Q) {
		DMatrixRMaj w_i = new DMatrixRMaj(3,3);

		for (int i = 0; i < cameras.size; i++) {
			computeW(cameras.get(i),Q,w_i);
			Intrinsic calib = solveForCalibration(w_i);
			if( sanityCheck(calib)) {
				solutions.add(calib);
			}
		}
	}
1116694
public DMatrixRMaj getQ(DMatrixRMaj Q ) {
        Q = UtilDecompositons_DDRM.checkIdentity(Q,N,N);

        for( int i = 0; i < N; i++ ) w[i] = 0;

        for( int j = N-2; j >= 0; j-- ) {
            w[j+1] = 1;
            for( int i = j+2; i < N; i++ ) {
                w[i] = QT.get(j,i);
            }
            QrHelperFunctions_DDRM.rank1UpdateMultR(Q, w, gammas[j + 1], j + 1, j + 1, N, b);
//            Q.print();
        }

        return Q;
    }
490009
def precompute_sharp_round(nxk, nyk, xc, yc):
    

    # Create arrays for the two- and four-fold symmetry computations:
    s4m = np.ones((nyk,nxk),dtype=np.int16)
    s4m[yc, xc] = 0

    s2m = np.ones((nyk,nxk),dtype=np.int16)
    s2m[yc, xc] = 0
    s2m[yc:nyk, 0:xc]     = -1;
    s2m[0:yc+1, xc+1:nxk] = -1;

    return s2m, s4m
1068849
public void recomputeQ( DMatrix4x4 Q ) {
		CommonOps_DDF3.multTransB(k,k,w);

		Q.a11 = w.a11;Q.a12 = w.a12;Q.a13 = w.a13;
		Q.a21 = w.a21;Q.a22 = w.a22;Q.a23 = w.a23;
		Q.a31 = w.a31;Q.a32 = w.a32;Q.a33 = w.a33;

		CommonOps_DDF3.mult(w,p,t);
		CommonOps_DDF3.scale(-1,t);

		Q.a14 = t.a1;Q.a24 = t.a2;Q.a34 = t.a3;
		Q.a41 = t.a1;Q.a42 = t.a2;Q.a43 = t.a3;

		Q.a44 = -CommonOps_DDF3.dot(t,p);
	}
1262849
public void householderSymmetric( int row , double gamma )
    {
        int startU = row*N;

        // compute v = -gamma*A*u
        for( int i = row+1; i < N; i++ ) {
            double totalReal = 0;
            double totalImag = 0;

            // the lower triangle is not written to so it needs to traverse upwards
            // to get the information.  Reduces the number of matrix writes need
            // improving large matrix performance
            for( int j = row+1; j < i; j++ ) {
                double realA = QT.data[(j*N+i)*2];
                double imagA = -QT.data[(j*N+i)*2+1];

                double realU = QT.data[(startU+j)*2];
                double imagU = QT.data[(startU+j)*2+1];

                totalReal += realA*realU - imagA*imagU;
                totalImag += realA*imagU + imagA*realU;
            }
            for( int j = i; j < N; j++ ) {
                double realA = QT.data[(i*N+j)*2];
                double imagA = QT.data[(i*N+j)*2+1];

                double realU = QT.data[(startU+j)*2];
                double imagU = QT.data[(startU+j)*2+1];

                totalReal += realA*realU - imagA*imagU;
                totalImag += realA*imagU + imagA*realU;
            }
            w[i*2]   = -gamma*totalReal;
            w[i*2+1] = -gamma*totalImag;
        }

        // alpha = -0.5*gamma*u^T*v
        double realAplha = 0;
        double imageAlpha = 0;

        for( int i = row+1; i < N; i++ ) {
            double realU = QT.data[(startU+i)*2];
            double imagU = -QT.data[(startU+i)*2+1];

            double realV = w[i*2];
            double imagV = w[i*2+1];

            realAplha += realU*realV - imagU*imagV;
            imageAlpha += realU*imagV + imagU*realV;
        }
        realAplha *= -0.5*gamma;
        imageAlpha *= -0.5*gamma;

        // w = v + alpha*u
        for( int i = row+1; i < N; i++ ) {
            double realU = QT.data[(startU+i)*2];
            double imagU = QT.data[(startU+i)*2+1];

            w[i*2]   += realAplha*realU - imageAlpha*imagU;
            w[i*2+1] += realAplha*imagU + imageAlpha*realU;
        }

        // A = A + w*u^T + u*w^T
        for( int i = row+1; i < N; i++ ) {

            double realWW = w[i*2];
            double imagWW = w[i*2+1];

            double realUU = QT.data[(startU+i)*2];
            double imagUU = QT.data[(startU+i)*2+1];

            int indA = (i*N+i)*2;
            for( int j = i; j < N; j++ ) {
                // only write to the upper portion of the matrix
                // this reduces the number of cache misses

                double realU = QT.data[(startU+j)*2];
                double imagU = -QT.data[(startU+j)*2+1];

                double realW = w[j*2];
                double imagW = -w[j*2+1];

                QT.data[indA++] += realWW*realU - imagWW*imagU + realW*realUU - imagW*imagUU;
                QT.data[indA++] += realWW*imagU + imagWW*realU + realW*imagUU + imagW*realUU;
            }
        }
    }
403717
function d3_geom_hullCCW(i1, i2, i3, v) {
  var t, a, b, c, d, e, f;
  t = v[i1]; a = t[0]; b = t[1];
  t = v[i2]; c = t[0]; d = t[1];
  t = v[i3]; e = t[0]; f = t[1];
  return ((f-b)*(c-a) - (d-b)*(e-a)) > 0;
}
25483
def dBinaryRochedr (r, D, q, F):
    

    r2 = (r*r).sum()
    r1 = np.sqrt(r2)

    return -1./r2 - q*(r1-r[0]/r1*D)/((r[0]-D)*(r[0]-D)+r[1]*r[1]+r[2]*r[2])**1.5 - q*r[0]/r1/D/D + F*F*(1+q)*(1-r[2]*r[2]/r2)*r1
1240381
def q12d(grid, spacing=None, E=1e5, nu=0.3, dirichlet_boundary=True,
         format=None):
    
    X, Y = tuple(grid)

    if X < 1 or Y < 1:
        raise ValueError('invalid grid shape')

    if dirichlet_boundary:
        X += 1
        Y += 1

    pts = np.mgrid[0:X+1, 0:Y+1]
    pts = np.hstack((pts[0].T.reshape(-1, 1) - X / 2.0,
                     pts[1].T.reshape(-1, 1) - Y / 2.0))

    if spacing is None:
        DX, DY = 1, 1
    else:
        DX, DY = tuple(spacing)
        pts *= [DX, DY]

    # compute local stiffness matrix
    lame = E * nu / ((1 + nu) * (1 - 2*nu))  # Lame's first parameter
    mu = E / (2 + 2*nu)                   # shear modulus

    vertices = np.array([[0, 0], [DX, 0], [DX, DY], [0, DY]])
    K = q12d_local(vertices, lame, mu)

    nodes = np.arange((X+1)*(Y+1)).reshape(X+1, Y+1)
    LL = nodes[:-1, :-1]
    Id = (2*LL).repeat(K.size).reshape(-1, 8, 8)
    J = Id.copy()
    Id += np.tile([0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8, 1))
    J += np.tile([0, 1, 2, 3, 2*X + 4, 2*X + 5, 2*X + 2, 2*X + 3], (8, 1)).T
    V = np.tile(K, (X*Y, 1))

    Id = np.ravel(Id)
    J = np.ravel(J)
    V = np.ravel(V)

    # sum duplicates
    A = coo_matrix((V, (Id, J)), shape=(pts.size, pts.size)).tocsr()
    A = A.tobsr(blocksize=(2, 2))

    del Id, J, V, LL, nodes

    B = np.zeros((2 * (X+1)*(Y+1), 3))
    B[0::2, 0] = 1
    B[1::2, 1] = 1
    B[0::2, 2] = -pts[:, 1]
    B[1::2, 2] = pts[:, 0]

    if dirichlet_boundary:
        mask = np.zeros((X+1, Y+1), dtype='bool')
        mask[1:-1, 1:-1] = True
        mask = np.ravel(mask)
        data = np.zeros(((X-1)*(Y-1), 2, 2))
        data[:, 0, 0] = 1
        data[:, 1, 1] = 1
        indices = np.arange((X-1)*(Y-1))
        indptr = np.concatenate((np.array([0]), np.cumsum(mask)))
        P = bsr_matrix((data, indices, indptr),
                       shape=(2*(X+1)*(Y+1), 2*(X-1)*(Y-1)))
        Pt = P.T
        A = P.T * A * P

        B = Pt * B

    return A.asformat(format), B
280319
def quatInv(q):
    """"""
    qConj = -q
    qConj[0] = -qConj[0]
    normSqr = multiplyQuats(q, qConj)[0]
    return qConj/normSqr
89504
def makeQ(r1, r2, r3, r4=0):
    
    Q = np.asarray([
        [r4, -r3, r2, r1],
        [r3, r4, -r1, r2],
        [-r2, r1, r4, r3],
        [-r1, -r2, -r3, r4]])
    return Q
987488
public static int cusolverSpScsrlsvqr(
        cusolverSpHandle handle, 
        int m, 
        int nnz, 
        cusparseMatDescr descrA, 
        Pointer csrVal, 
        Pointer csrRowPtr, 
        Pointer csrColInd, 
        Pointer b, 
        float tol, 
        int reorder, 
        Pointer x, 
        int[] singularity)
    {
        return checkResult(cusolverSpScsrlsvqrNative(handle, m, nnz, descrA, csrVal, csrRowPtr, csrColInd, b, tol, reorder, x, singularity));
    }
154021
@Override
    public QR qr() {
        // Initialize.
        int m = nrows();
        int n = ncols();
        double[] rDiagonal = new double[n];

        // Main loop.
        for (int k = 0; k < n; k++) {
            // Compute 2-norm of k-th column without under/overflow.
            double nrm = 0.0;
            for (int i = k; i < m; i++) {
                nrm = Math.hypot(nrm, get(i, k));
            }

            if (nrm != 0.0) {
                // Form k-th Householder vector.
                if (get(k, k) < 0) {
                    nrm = -nrm;
                }
                for (int i = k; i < m; i++) {
                    div(i, k, nrm);
                }
                add(k, k, 1.0);

                // Apply transformation to remaining columns.
                for (int j = k + 1; j < n; j++) {
                    double s = 0.0;
                    for (int i = k; i < m; i++) {
                        s += get(i, k) * get(i, j);
                    }
                    s = -s / get(k, k);
                    for (int i = k; i < m; i++) {
                        add(i, j, s * get(i, k));
                    }
                }
            }
            rDiagonal[k] = -nrm;
        }

        boolean singular = false;
        for (int j = 0; j < rDiagonal.length; j++) {
            if (rDiagonal[j] == 0) {
                singular = true;
                break;
            }
        }

        return new QR(this, rDiagonal, singular);
    }
298315
def rsdl(self):
        """"""

        return np.linalg.norm((self.X - self.Yprv).ravel())
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;
    }
  }
1175084
def _map_tril_1d_on_2d(indices, dims):
    

    N = (dims * dims - dims) / 2

    m = np.ceil(np.sqrt(2 * N))
    c = m - np.round(np.sqrt(2 * (N - indices))) - 1
    r = np.mod(indices + (c + 1) * (c + 2) / 2 - 1, m) + 1

    return np.array([r, c], dtype=np.int64)
699930
def _multicomplex2(f, fx, x, h):
        """"""
        n = len(x)
        ee = np.diag(h)
        hess = np.outer(h, h)
        cmplx_wrap = Bicomplex.__array_wrap__
        for i in range(n):
            for j in range(i, n):
                zph = Bicomplex(x + 1j * ee[i, :], ee[j, :])
                hess[i, j] = cmplx_wrap(f(zph)).imag12 / hess[j, i]
                hess[j, i] = hess[i, j]
        return hess