Rank-k operation on real symmetric matrix.
972074
def cublasSsyrk(handle, uplo, trans, n, k, alpha, A, lda, beta, C, ldc):
    
    
    status = _libcublas.cublasSsyrk_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       n, k, ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda, 
                                       ctypes.byref(ctypes.c_float(beta)), 
                                       int(C), ldc)
    cublasCheckStatus(status)
96453
def _eig_complex_symmetric(M: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    
    if not np.allclose(M, M.transpose()):
        raise np.linalg.LinAlgError('Not a symmetric matrix')

    # The matrix of eigenvectors should be orthogonal.
    # But the standard 'eig' method will fail to return an orthogonal
    # eigenvector matrix when the eigenvalues are degenerate. However,
    # both the real and
    # imaginary part of M must be symmetric with the same orthogonal
    # matrix of eigenvectors. But either the real or imaginary part could
    # vanish. So we use a randomized algorithm where we diagonalize a
    # random linear combination of real and imaginary parts to find the
    # eigenvectors, taking advantage of the 'eigh' subroutine for
    # diagonalizing symmetric matrices.
    # This can fail if we're very unlucky with our random coefficient, so we
    # give the algorithm a few chances to succeed.

    # Empirically, never seems to fail on randomly sampled complex
    # symmetric 4x4 matrices.
    # If failure rate is less than 1 in a million, then 16 rounds
    # will have overall failure rate less than 1 in a googol.
    # However, cannot (yet) guarantee that there aren't special cases
    # which have much higher failure rates.

    # GEC 2018

    max_attempts = 16
    for _ in range(max_attempts):
        c = np.random.uniform(0, 1)
        matrix = c * M.real + (1-c) * M.imag
        _, eigvecs = np.linalg.eigh(matrix)
        eigvecs = np.array(eigvecs, dtype=complex)
        eigvals = np.diag(eigvecs.transpose() @ M @ eigvecs)

        # Finish if we got a correct answer.
        reconstructed = eigvecs @ np.diag(eigvals) @ eigvecs.transpose()
        if np.allclose(M, reconstructed):
            return eigvals, eigvecs

    # Should never happen. Hopefully.
    raise np.linalg.LinAlgError(
        'Cannot diagonalize complex symmetric matrix.')
1263626
private void orthes(Matrix h, Matrix v, Vector ort) {
        // 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 n = ort.length();
        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.get(i, m - 1));
            }

            if (scale != 0.0) {

                // Compute Householder transformation.

                double hh = 0.0;
                for (int i = high; i >= m; i--) {
                    ort.set(i, h.get(i, m - 1) / scale);
                    hh += ort.get(i) * ort.get(i);
                }

                double g = Math.sqrt(hh);

                if (ort.get(m) > Matrices.EPS) {
                    g = -g;
                }

                hh = hh - ort.get(m) * g;
                ort.updateAt(m, Vectors.asMinusFunction(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.get(i) * h.get(i, j);
                    }
                    f = f / hh;
                    for (int i = m; i <= high; i++) {
                        h.updateAt(i, j, Matrices.asMinusFunction(f * ort.get(i)));
                    }
                }

                for (int i = 0; i <= high; i++) {
                    double f = 0.0;
                    for (int j = high; j >= m; j--) {
                        f += ort.get(j) * h.get(i, j);
                    }
                    f = f / hh;
                    for (int j = m; j <= high; j++) {
                        h.updateAt(i, j, Matrices.asMinusFunction(f * ort.get(j)));
                    }
                }
                ort.set(m, scale * ort.get(m));
                h.set(m,  m - 1, scale * g);
            }
        }

        // Accumulate transformations (Algol's ortran).

        for (int m = high - 1; m >= low + 1; m--) {
            if (Math.abs(h.get(m, m - 1)) > Matrices.EPS) {
                for (int i = m + 1; i <= high; i++) {
                    ort.set(i, h.get(i, m - 1));
                }
                for (int j = m; j <= high; j++) {
                    double g = 0.0;
                    for (int i = m; i <= high; i++) {
                        g += ort.get(i) * v.get(i, j);
                    }
                    // Double division avoids possible underflow
                    g = (g / ort.get(m)) / h.get(m, m - 1);
                    for (int i = m; i <= high; i++) {
                        v.updateAt(i, j, Matrices.asPlusFunction(g * ort.get(i)));
                    }
                }
            }
        }
    }
465968
private static final int bestPivot(int rank, int m1, int m2, int m3, int m4, int m5) {
    if(rank < m1) {
      return m1;
    }
    if(rank > m5) {
      return m5;
    }
    if(rank < m2) {
      return m2;
    }
    if(rank > m4) {
      return m4;
    }
    return m3;
  }
163316
def Ak(self, k, m_pred, P_pred):  # returns state iteration matrix
        

        return self.A[:, :, self.index[self.A_time_var_index, k]]
116888
public static void symmRankTranA_sub(DMatrixRMaj a , DMatrixRMaj c ,
                                         int startIndexC )
    {
        // TODO update so that it doesn't modify/read parts that it shouldn't
        final double dataA[] = a.data;
        final double dataC[] = c.data;

//        for( int i = 0; i < a.numCols; i++ ) {
//            for( int k = 0; k < a.numRows; k++ ) {
//                double valA = dataA[k*a.numCols+i];
//
//                for( int j = i; j < a.numCols; j++ ) {
//                    dataC[startIndexC+i*c.numCols+j] -= valA * dataA[k*a.numCols+j];
//                }
//            }
//        }

        final int strideC = c.numCols + 1;
        for( int i = 0; i < a.numCols; i++ ) {
            int indexA = i;
            int endR = a.numCols;

            for( int k = 0; k < a.numRows; k++ , indexA += a.numCols , endR += a.numCols) {
                int indexC = startIndexC;
                final double valA = dataA[indexA];
                int indexR = indexA;

                while( indexR < endR ) {
                    dataC[indexC++] -= valA * dataA[indexR++];
                }
            }
            startIndexC += strideC;
        }

    }
67264
function(k, U) {
            	var U = U || pca$1(this).U;
            	var Ureduce= U.slice(1, U.rows(), 1, k);
            	return {Z: this.x(Ureduce), U: U};
                }
50673
def schur_complement(mat, row, col):
    
    a = mat[:row, :col]
    b = mat[:row, col:]
    c = mat[row:, :col]
    d = mat[row:, col:]
    return a - b.dot(d.I).dot(c)
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);
	}
906378
def K(self,X,X2,target):
        """"""
        AX = np.dot(X,self.transform)
        if X2 is None:
            X2 = X
            AX2 = AX
        else:
            AX2 = np.dot(X2, self.transform)
        self.k.K(X,X2,target)
        self.k.K(AX,X2,target)
        self.k.K(X,AX2,target)
        self.k.K(AX,AX2,target)
258913
def K(self, parm):
        
        return RQ_K_matrix(self.X, parm) + np.identity(self.X.shape[0])*(10**-10)
1237858
def cublasStbmv(handle, uplo, trans, diag, n, k, A, lda, x, incx):
    
    
    status = _libcublas.cublasStbmv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       n, k, int(A), lda, int(x), incx)
    cublasCheckStatus(status)
682373
def _build_matrix_non_uniform(p, q, coords, k):
    """"""
    A = [[1] * (p+q+1)]
    for i in range(1, p + q + 1):
        line = [(coords[k+j] - coords[k])**i for j in range(-p, q+1)]
        A.append(line)
    return np.array(A)
395013
def skew_matrix(w):
    ''''''
    return np.array([[0, -w[2], w[1]],
                     [w[2], 0, -w[0]],
                     [-w[1], w[0], 0]])
68520
private Matrix calculateEdgeSimMatrix(
            final List<Edge> edgeList, final SparseMatrix sm) {

        final int numEdges = edgeList.size();
        final Matrix edgeSimMatrix = 
            new SparseSymmetricMatrix(
                new SparseHashMatrix(numEdges, numEdges));

        Object key = workQueue.registerTaskGroup(numEdges);
        for (int i = 0; i < numEdges; ++i) {
            final int row = i;
            workQueue.add(key, new Runnable() {
                    public void run() {
                        for (int j = row; j < numEdges; ++j) {
                            Edge e1 = edgeList.get(row);
                            Edge e2 = edgeList.get(j);
                            
                            double sim = getEdgeSimilarity(sm, e1, e2);
                
                            if (sim > 0) {
                                // The symmetric matrix handles the (j,i) case
                                edgeSimMatrix.set(row, j, sim);
                            }
                        }
                    }
                });            
        }
        workQueue.await(key);
        return edgeSimMatrix;
    }
929959
def cublasZsyrk(handle, uplo, trans, n, k, alpha, A, lda, beta, C, ldc):
    
    
    status = _libcublas.cublasZsyrk_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       n, k, ctypes.byref(cuda.cuDoubleComplex(alpha.real,    
                                                                               alpha.imag)),
                                       int(A), lda,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real,
                                                                         beta.imag)),
                                       int(C), ldc)
    cublasCheckStatus(status)
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);
            }
        }
    }