Matrix-vector product for complex general matrix.
552428
def cublasCgemv(handle, trans, m, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasCgemv_v2(handle,
                                       _CUBLAS_OP[trans], m, n,
                                       ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(y), incy)
    cublasCheckStatus(status)
1361461
def eighCovarianceMatrix(cfile):
    
    # precompute eigenvalue decomposition
    K = np.loadtxt(cfile+'.cov')
    K+= 1e-4*sp.eye(K.shape[0])
    S,U = la.eigh(K); S=S[::-1]; U=U[:,::-1]
    np.savetxt(cfile+'.cov.eval',S,fmt='%.6f')
    np.savetxt(cfile+'.cov.evec',U,fmt='%.6f')
964528
def cublasZsymv(handle, uplo, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasZsymv_v2(handle, 
                                        _CUBLAS_FILL_MODE[uplo], n, 
                                        ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                            alpha.imag)), 
                                        int(A), lda, int(x), incx, 
                                        ctypes.byref(cuda.cuDoubleComplex(beta.real,
                                                                            beta.imag)), 
                                        int(y), incy)
    cublasCheckStatus(status)
989122
def G(w, Xs):
    ''''''
    n = len(Xs)
    P = projection_matrix(w)
    Ys = [np.dot(P, X) for X in Xs]
    A = calc_A(Ys)
    A_hat = calc_A_hat(A, skew_matrix(w))

    
    u = sum(np.dot(Y, Y) for Y in Ys) / n
    v = np.dot(A_hat, sum(np.dot(Y, Y) * Y for Y in Ys)) / np.trace(np.dot(A_hat, A))

    return sum((np.dot(Y, Y) - u - 2 * np.dot(Y, v)) ** 2 for Y in Ys)
1127369
def vector_projection(v1, v2):
    '''
    '''
    return scalar_projection(v1, v2) * v2 / np.linalg.norm(v2)
1195500
static void compute_B(DMatrix2x2 B , DMatrixRMaj R_v , double v1 , double v2 ) {
		B.a11 = R_v.data[0]+R_v.data[6]*-v1;
		B.a12 = R_v.data[1]+R_v.data[7]*-v1;
		B.a21 = R_v.data[3]+R_v.data[6]*-v2;
		B.a22 = R_v.data[4]+R_v.data[7]*-v2;
	}
1064236
def _compute_C_matrix(self):
        
        self.C = np.repeat(np.diag(self.Lp)[:, np.newaxis],
                           self.Lp.shape[0], axis=1)
        self.C += np.repeat(np.diag(self.Lp)[np.newaxis, :],
                            self.Lp.shape[0], axis=0)
        self.C -= 2*self.Lp
        # the following is much slower
        # self.C = np.zeros(self.Lp.shape)
        # for i in range(self.Lp.shape[0]):
        #     for j in range(self.Lp.shape[1]):
        #         self.C[i, j] = self.Lp[i, i] + self.Lp[j, j] - 2*self.Lp[i, j]
        volG = np.sum(self.z)
        self.C *= volG
        settings.mt(0, 'computed commute distance matrix')
        self.distances_dpt = self.C
103306
def cublasStrmm(handle, side, uplo, trans, diag, m, n, alpha, A, lda, B, ldb, C, ldc):
    
    
    status = _libcublas.cublasStrmm_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, int(C), ldc)
    cublasCheckStatus(status)
1169693
public Complex contraction() {
        int i, j;
        Complex result = new Complex(0d, 0d);
        for (i = 0; i < rows; i++)
            for (j = 0; j < columns; j++) {
                result.real += realmatrix[i][j];
                result.imag += imagmatrix[i][j];
            }
        return result;
    }
489044
public IMatrix mul(Complex a) {
        IMatrix result = new IMatrix(rows, columns);
        mul(a, result);
        return result;
    }
1266891
public static ComplexDoubleMatrix gerc(ComplexDouble alpha, ComplexDoubleMatrix x,
			ComplexDoubleMatrix y, ComplexDoubleMatrix a) {
		NativeBlas.zgerc(a.rows, a.columns, alpha, x.data, 0, 1, y.data, 0, 1, a.data,
				0, a.rows);
		return a;
	}
304697
def matvec(a, b, compression=False):
    """"""
    acrs = _vector.vector.to_list(a.tt)
    bcrs = _vector.vector.to_list(b)
    ccrs = []
    d = b.d

    def get_core(i):
        acr = _np.reshape(
            acrs[i],
            (a.tt.r[i],
             a.n[i],
             a.m[i],
             a.tt.r[
                 i + 1]),
            order='F')
        acr = acr.transpose([3, 0, 1, 2])  # a(R_{i+1}, R_i, n_i, m_i)
        bcr = bcrs[i].transpose([1, 0, 2])  # b(m_i, r_i, r_{i+1})
        # c(R_{i+1}, R_i, n_i, r_i, r_{i+1})
        ccr = _np.tensordot(acr, bcr, axes=(3, 0))
        ccr = ccr.transpose([1, 3, 2, 0, 4]).reshape(
            (a.tt.r[i] * b.r[i], a.n[i], a.tt.r[i + 1] * b.r[i + 1]), order='F')
        return ccr

    if compression:  # the compression is laaaaazy and one-directioned
        # calculate norm of resulting _vector first
        nrm = _np.array([[1.0]])  # 1 x 1
        v = _np.array([[1.0]])
        for i in xrange(d):
            ccr = get_core(i)
            # print(str(ccr.shape) + " -> "),
            # minimal loss compression
            ccr = _np.tensordot(v, ccr, (1, 0))
            rl, n, rr = ccr.shape
            if i < d - 1:
                u, s, v = _np.linalg.svd(
                    ccr.reshape(
                        (rl * n, rr), order='F'), full_matrices=False)
                newr = min(rl * n, rr)
                ccr = u[:, :newr].reshape((rl, n, newr), order='F')
                v = _np.dot(_np.diag(s[:newr]), v[:newr, :])
            # print(ccr.shape)
            # r x r . r x n x R -> r x n x R
            nrm = _np.tensordot(nrm, ccr, (0, 0))
            # r x n x R . r x n x R -> n x R x n x R
            nrm = _np.tensordot(nrm, _np.conj(ccr), (0, 0))
            nrm = nrm.diagonal(axis1=0, axis2=2)  # n x R x n x R -> R x R x n
            nrm = nrm.sum(axis=2)  # R x R x n -> R x R
        if nrm.size > 1:
            raise Exception('too many numbers in norm')
        # print("Norm calculated:", nrm)
        nrm = _np.sqrt(_np.linalg.norm(nrm))
        # print("Norm predicted:", nrm)
        compression = compression * nrm / _np.sqrt(d - 1)
        v = _np.array([[1.0]])

    for i in xrange(d):
        ccr = get_core(i)
        rl, n, rr = ccr.shape
        if compression:
            ccr = _np.tensordot(v, ccr, (1, 0))  # c(s_i, n_i, r_i, r_{i+1})
            if i < d - 1:
                rl = v.shape[0]
                u, s, v = _np.linalg.svd(
                    ccr.reshape(
                        (rl * n, rr), order='F'), full_matrices=False)
                ss = _np.cumsum(s[::-1])[::-1]
                newr = max(min([r for r in range(ss.size) if ss[
                    r] <= compression] + [min(rl * n, rr)]), 1)
                # print("Rank % 4d replaced by % 4d" % (rr, newr))
                ccr = u[:, :newr].reshape((rl, n, newr), order='F')
                v = _np.dot(_np.diag(s[:newr]), v[:newr, :])
        ccrs.append(ccr)
    result = _vector.vector.from_list(ccrs)
    if compression:
        # print(result)
        print("Norm actual:", result.norm(), " mean rank:", result.rmean())
        # print("Norm very actual:", matvec(a,b).norm())
    return result
719273
def cublasCgbmv(handle, trans, m, n, kl, ku, alpha, A, lda,
                x, incx, beta, y, incy):
    

    status = _libcublas.cublasCgbmv_v2(handle,
                                       trans, m, n, kl, ku,
                                       ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(y), incy)
    cublasCheckStatus(status)
395013
def skew_matrix(w):
    ''''''
    return np.array([[0, -w[2], w[1]],
                     [w[2], 0, -w[0]],
                     [-w[1], w[0], 0]])
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];
}
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)
1001103
def _complex_even(f, fx, x, h):
        
        n = len(x)
        ee = np.diag(h)
        hess = 2. * np.outer(h, h)
        for i in range(n):
            for j in range(i, n):
                hess[i, j] = (f(x + 1j * ee[i] + ee[j])
                              - f(x + 1j * ee[i] - ee[j])).imag / hess[j, i]
                hess[j, i] = hess[i, j]
        return hess