Rank-k operation on Hermitian matrix.
353641
def cublasZherk(handle, uplo, trans, n, k, alpha, A, lda, beta, C, ldc):
    
    
    status = _libcublas.cublasZherk_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans],
                                       n, k, ctypes.byref(ctypes.c_double(alpha)),
                                       int(A), lda, 
                                       ctypes.byref(ctypes.c_double(beta)),
                                       int(C), ldc)
    cublasCheckStatus(status)
1340920
private static Matrix permuteMatrix(Matrix mat, int cpRows, int cpCols) {
		//fix up cp to be positive
		if(cpRows == 0 && cpCols == 0) {
			return mat.copy();
		}
		if(cpRows < 0) {
			cpRows = mat.getRowDimension()+cpRows;
		}
		if(cpRows < 0 || cpRows >= mat.getRowDimension()) {
			throw new ArrayIndexOutOfBoundsException( String.format(
					"Can't permute rows by %d: only %d rows.",
					cpRows, mat.getRowDimension() )
			);
		}

		if(cpCols < 0) {
			cpCols = mat.getColumnDimension()+cpCols;
		}
		if(cpCols < 0 || cpCols >= mat.getColumnDimension()) {
			throw new ArrayIndexOutOfBoundsException( String.format(
					"Can't permute cols by %d: only %d rows.",
					cpCols, mat.getColumnDimension() )
			);
		}

		int[] rows = new int[mat.getRowDimension()];
		for(int i=0;i<rows.length;i++) {
			rows[i] = (i+cpRows)%rows.length;
		}
		int[] cols = new int[mat.getColumnDimension()];
		for(int i=0;i<cols.length;i++) {
			cols[i] = (i+cpCols)%cols.length;
		}

		Matrix newMat = mat.getMatrix(rows, cols);
		assert(newMat.getRowDimension() == mat.getRowDimension());
		assert(newMat.getColumnDimension() == mat.getColumnDimension());
		assert(newMat.get(0, 0) ==
			mat.get(cpRows%mat.getRowDimension(), cpCols%mat.getColumnDimension()));


		return newMat;
	}
1111590
private void gather(double[] z, SparseVector v, double taui, int d) {
        // Number of entries in the lower and upper part of the original matrix
        int nl = 0, nu = 0;
        for (VectorEntry e : v) {
            if (e.index() < d)
                nl++;
            else if (e.index() > d)
                nu++;
        }
        v.zero();

        // Entries in the L part of the vector
        lower.clear();
        for (int i = 0; i < d; ++i)
            if (Math.abs(z[i]) > taui)
                lower.add(new IntDoubleEntry(i, z[i]));

        // Entries in the U part of the vector
        upper.clear();
        for (int i = d + 1; i < z.length; ++i)
            if (Math.abs(z[i]) > taui)
                upper.add(new IntDoubleEntry(i, z[i]));

        // Sort in descending order
        Collections.sort(lower);
        Collections.sort(upper);

        // Always keep the diagonal
        v.set(d, z[d]);

        // Keep at most nl+p lower entries
        for (int i = 0; i < Math.min(nl + p, lower.size()); ++i) {
            IntDoubleEntry e = lower.get(i);
            v.set(e.index, e.value);
        }

        // Keep at most nu+p upper entries
        for (int i = 0; i < Math.min(nu + p, upper.size()); ++i) {
            IntDoubleEntry e = upper.get(i);
            v.set(e.index, e.value);
        }
    }
291673
private void tred2()
    {
        for(int j = 0; j < n; j++)
            d[j] = V.get(n-1, j);

        // Householder reduction to tridiagonal form.

        for (int i = n - 1; i > 0; i--)
        {

            // Scale to avoid under/overflow.

            double scale = 0.0;
            double h = 0.0;
            for (int k = 0; k < i; k++)
            {
                scale = scale + abs(d[k]);
            }
            if (scale == 0.0)
            {
                e[i] = d[i - 1];
                for (int j = 0; j < i; j++)
                {
                    d[j] = V.get(i-1, j);
                    V.set(i, j, 0.0);
                    V.set(j, i, 0.0);
                }
            }
            else
            {

                // Generate Householder vector.

                for (int k = 0; k < i; k++)
                {
                    d[k] /= scale;
                    h += d[k] * d[k];
                }
                double f = d[i - 1];
                double g = sqrt(h);
                if (f > 0)
                    g = -g;
                e[i] = scale * g;
                h -= f * g;
                d[i - 1] = f - g;
                Arrays.fill(e, 0, i, 0.0);

                // Apply similarity transformation to remaining columns.

                for (int j = 0; j < i; j++)
                {
                    f = d[j];
                    V.set(j, i, f);
                    g = e[j] + V.get(j, j) * f;
                    for (int k = j + 1; k <= i - 1; k++)
                    {
                        g += V.get(k,j) * d[k];
                        e[k] += V.get(k, j) * f;
                    }
                    e[j] = g;
                }
                f = 0.0;
                for (int j = 0; j < i; j++)
                {
                    e[j] /= h;
                    f += e[j] * d[j];
                }
                double hh = f / (h + h);
                for (int j = 0; j < i; j++)
                {
                    e[j] -= hh * d[j];
                }
                for (int j = 0; j < i; j++)
                {
                    f = d[j];
                    g = e[j];
                    
                    for (int k = j; k <= i - 1; k++)
                    {
                        V.increment(k, j, -(f * e[k] + g * d[k]));
                    }
                    d[j] = V.get(i-1, j);
                    V.set(i, j, 0.0);
                }
            }
            d[i] = h;
        }


        // Accumulate transformations.

        for (int i = 0; i < n - 1; i++)
        {
            V.set(n-1, i, V.get(i, i));
            V.set(i, i, 1.0);
            double h = d[i + 1];
            if (h != 0.0)
            {
                for (int k = 0; k <= i; k++)
                {
                    d[k] = V.get(k, i+1) / h;
                }
                for (int j = 0; j <= i; j++)
                {
                    double g = 0.0;
                    for (int k = 0; k <= i; k++)
                    {
                        g += V.get(k, i+1) * V.get(k, j);
                    }
                    
                    RowColumnOps.addMultCol(V, j, 0, i+1, -g, d);
                }
            }

            RowColumnOps.fillCol(V, i+1, 0, i+1, 0.0);
            
        }
        for (int j = 0; j < n; j++)
        {
            d[j] = V.get(n-1, j);
            V.set(n-1, j, 0.0);
        }
        V.set(n-1, n-1, 1.0);
        e[0] = 0.0;
    }
1303325
public static void symmRankNMinus_U(int blockLength ,
                                        DSubmatrixD1 A , DSubmatrixD1 B )
    {

        int heightB = B.row1-B.row0;
        if( heightB > blockLength )
            throw new IllegalArgumentException("Height of B cannot be greater than the block length");

        int N = B.col1-B.col0;

        if( A.col1-A.col0 != N )
            throw new IllegalArgumentException("A does not have the expected number of columns based on B's width");
        if( A.row1-A.row0 != N )
            throw new IllegalArgumentException("A does not have the expected number of rows based on B's width");


        for( int i = B.col0; i < B.col1; i += blockLength ) {

            int indexB_i = B.row0*B.original.numCols + i*heightB;
            int widthB_i = Math.min(blockLength,B.col1-i);

            int rowA = i-B.col0+A.row0;
            int heightA = Math.min( blockLength , A.row1 - rowA);

            for( int j = i; j < B.col1; j += blockLength ) {

                int widthB_j = Math.min(blockLength,B.col1-j);

                int indexA = rowA * A.original.numCols + (j-B.col0+A.col0)*heightA;
                int indexB_j = B.row0*B.original.numCols + j*heightB;

                if( i == j ) {
                    // only the upper portion of this block needs to be modified since it is along a diagonal
                    multTransABlockMinus_U( B.original.data,A.original.data,
                            indexB_i,indexB_j,indexA,heightB,widthB_i,widthB_j);
                } else {
                    multTransABlockMinus( B.original.data,A.original.data,
                            indexB_i,indexB_j,indexA,heightB,widthB_i,widthB_j);
                }
            }
        }
    }
650218
def solve_factorized_schur(z, Fval, LU, G, A):
    M, N=G.shape
    P, N=A.shape

    """"""
    m = M    

    """Primal variable"""
    x = z[0:N]

    """Multiplier for equality constraints"""
    nu = z[N:N+P]

    """Multiplier for inequality constraints"""
    l = z[N+P:N+P+M]

    """Slacks"""
    s = z[N+P+M:]

    """Dual infeasibility"""
    rd = Fval[0:N]
    
    """Primal infeasibility"""
    rp1 = Fval[N:N+P]
    rp2 = Fval[N+P:N+P+M]

    """Centrality"""
    rc = Fval[N+P+M:]

    """Sigma matrix"""
    SIG = diags(l/s, 0)

    """Assemble right hand side of augmented system"""
    r1 = rd + mydot(G.T, mydot(SIG, rp2)) - mydot(G.T, rc/s)
    r2 = rp1

    """Unpack LU-factors"""
    LU_S, LU_H = LU

    """Assemble right hand side for normal equation"""
    b = r2 - mydot(A, mysolve(LU_H, r1))  

    """Solve for dnu"""
    dnu = mysolve(LU_S, b)
       
    """Solve for dx"""
    dx = mysolve(LU_H, -(r1 + mydot(A.T, dnu)))    
    
    """Obtain search directions for l and s"""
    ds = -rp2 - mydot(G, dx)
    dl = -mydot(SIG, ds) - rc/s

    dz = np.hstack((dx, dnu, dl, ds))
    return dz
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)
417841
private static void balbak(DenseMatrix V, double[] scale) {
        int n = V.nrows();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                V.mul(i, j, scale[i]);
            }
        }
    }
916154
def choi(self) -> bk.BKTensor:
        """"""
        # Put superop axes in [ok, ib, ob, ik] and reshape to matrix
        N = self.qubit_nb
        return bk.reshape(self.sharp.tensor, [2**(N*2)] * 2)
1035019
private static int[] elmhes(DenseMatrix A) {
        int n = A.nrows();
        int[] perm = new int[n];

        for (int m = 1; m < n - 1; m++) {
            double x = 0.0;
            int i = m;
            for (int j = m; j < n; j++) {
                if (Math.abs(A.get(j, m - 1)) > Math.abs(x)) {
                    x = A.get(j, m - 1);
                    i = j;
                }
            }
            perm[m] = i;
            if (i != m) {
                for (int j = m - 1; j < n; j++) {
                    double swap = A.get(i, j);
                    A.set(i, j, A.get(m, j));
                    A.set(m, j, swap);
                }
                for (int j = 0; j < n; j++) {
                    double swap = A.get(j, i);
                    A.set(j, i, A.get(j, m));
                    A.set(j, m, swap);
                }
            }
            if (x != 0.0) {
                for (i = m + 1; i < n; i++) {
                    double y = A.get(i, m - 1);
                    if (y != 0.0) {
                        y /= x;
                        A.set(i, m - 1, y);
                        for (int j = m; j < n; j++) {
                            A.sub(i, j, y * A.get(m, j));
                        }
                        for (int j = 0; j < n; j++) {
                            A.add(j, m, y * A.get(j, i));
                        }
                    }
                }
            }
        }

        return perm;
    }
1271544
def K(self, X, X2, target):
        """"""
        if (X2 is None) or (X2 is X):
            target[np.diag_indices_from(target)] += self._Kdiag(X)
108359
def cublasDsyr2k(handle, uplo, trans, n, k, alpha, A, lda, B, ldb, beta, C, ldc):
    

    status = _libcublas.cublasDsyr2k_v2(handle, 
                                        _CUBLAS_FILL_MODE[uplo], 
                                        _CUBLAS_OP[trans], 
                                        n, k, ctypes.byref(ctypes.c_double(alpha)),
                                        int(A), lda, int(B), ldb, 
                                        ctypes.byref(ctypes.c_double(beta)), 
                                        int(C), ldc)
    cublasCheckStatus(status)
204942
def _svd_sym_koopman(K, C00_train, Ctt_train):
    
    from pyemma._ext.variational.solvers.direct import spd_inv_sqrt
    # reweight operator to empirical distribution
    C0t_re = mdot(C00_train, K)
    # symmetrized operator and SVD
    K_sym = mdot(spd_inv_sqrt(C00_train), C0t_re, spd_inv_sqrt(Ctt_train))
    U, S, Vt = np.linalg.svd(K_sym, compute_uv=True, full_matrices=False)
    # projects back to singular functions of K
    U = mdot(spd_inv_sqrt(C00_train), U)
    Vt = mdot(Vt,spd_inv_sqrt(Ctt_train))
    return U, S, Vt.T
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;
        }

    }
1167745
def is_hermitian(operator):
    
    if isinstance(operator, qt.Qobj):
        return (operator.dag() - operator).norm(FROBENIUS) / operator.norm(FROBENIUS) < EPS
    if isinstance(operator, np.ndarray):
        return np.linalg.norm(operator.T.conj() - operator) / np.linalg.norm(operator) < EPS
    return spnorm(operator.H - operator) / spnorm(operator) < EPS
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);
	}
1210915
def _K(m):
    
    M = m*(m - 1)//2
    K = np.zeros((M, m**2), dtype=np.int64)
    row = 0
    for j in range(1, m):
        col = (j - 1)*m + j
        s = m - j
        K[row:(row+s), col:(col+s)] = np.eye(s)
        row += s
    return K