Matrix-vector product for real symmetric-packed matrix.
650
def cublasSspmv(handle, uplo, n, alpha, AP, x, incx, beta, y, incy):
    

    status = _libcublas.cublasSspmv_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       n,
                                       ctypes.byref(ctypes.c_float(alpha)),
                                       ctypes.byref(ctypes.c_float(AP)),
                                       int(x),
                                       incx,
                                       ctypes.byref(ctypes.c_float(beta)),
                                       int(y),
                                       incy)
    cublasCheckStatus(status)
1083012
def cublasSsbmv(handle, uplo, n, k, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasSsbmv_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], n, k,
                                       ctypes.byref(ctypes.c_float(alpha)),
                                       int(A), lda, int(x), incx,
                                       ctypes.byref(ctypes.c_float(beta)),
                                       int(y), incy)
    cublasCheckStatus(status)
1218098
def cublasDtpsv(handle, uplo, trans, diag, n, AP, x, incx):
    

    status = _libcublas.cublasDtpsv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       n, int(AP), int(x), incx)
    cublasCheckStatus(status)
86373
def cublasCtpmv(handle, uplo, trans, diag, n, AP, x, incx):
    
    
    status = _libcublas.cublasCtpmv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       n, int(AP), int(x), incx)
    cublasCheckStatus(status)
396675
def cublasCsymv(handle, uplo, n, alpha, A, lda, x, incx, beta, y, incy):
    

    status = _libcublas.cublasCsymv_v2(handle, 
                                        _CUBLAS_FILL_MODE[uplo], 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)
988274
def cublasCher2k(handle, uplo, trans, n, k, alpha, A, lda, B, ldb, beta, C, ldc):
    
    
    status = _libcublas.cublasCher2k_v2(handle, 
                                        _CUBLAS_FILL_MODE[uplo], 
                                        _CUBLAS_OP[trans], 
                                        n, k, ctypes.byref(cuda.cuFloatComplex(alpha.real,                 
                                                                               alpha.imag)),
                                        int(A), lda, int(B), ldb, 
                                        ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                         beta.imag)),
                                        int(C), ldc)
    cublasCheckStatus(status)
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)
193322
def cublasChpr2(handle, uplo, n, alpha, x, inx, y, incy, AP):
    

    status = _libcublas.cublasChpr2_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       n, ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                           alpha.imag)),
                                       int(x), incx, int(y), incy, int(AP))
    cublasCheckStatus(status)
942549
def d2Sbr_dV2(Cbr, Ybr, V, lam):
    
    nb = len(V)

    diaglam = spdiag(lam)
    diagV = spdiag(V)

    A = Ybr.H * diaglam * Cbr
    B = conj(diagV) * A * diagV
    D = spdiag(mul((A*V), conj(V)))
    E = spdiag(mul((A.T * conj(V)), V))
    F = B + B.T
    G = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))

    Haa = F - D - E
    Hva = 1j * G * (B - B.T - D + E)
    Hav = Hva.T
    Hvv = G * F * G

    return Haa, Hav, Hva, Hvv
221722
private void svdImpute(double[][] raw, double[][] data) {
        SVD svd = Matrix.newInstance(data).svd();

        int d = data[0].length;

        for (int i = 0; i < raw.length; i++) {
            int missing = 0;
            for (int j = 0; j < d; j++) {
                if (Double.isNaN(raw[i][j])) {
                    missing++;
                } else {
                    data[i][j] = raw[i][j];
                }
            }

            if (missing == 0) {
                continue;
            }

            DenseMatrix A = Matrix.zeros(d - missing, k);
            double[] b = new double[d - missing];

            for (int j = 0, m = 0; j < d; j++) {
                if (!Double.isNaN(raw[i][j])) {
                    for (int l = 0; l < k; l++) {
                        A.set(m, l, svd.getV().get(j, l));
                    }
                    b[m++] = raw[i][j];
                }
            }

            double[] s = new double[k];
            QR qr = A.qr();
            qr.solve(b, s);

            for (int j = 0; j < d; j++) {
                if (Double.isNaN(raw[i][j])) {
                    data[i][j] = 0;
                    for (int l = 0; l < k; l++) {
                        data[i][j] += s[l] * svd.getV().get(j, l);
                    }
                }
            }
        }
    }
625225
function matrix(rows, columns, values) {
    // console.log(values);
    const res = [];
    let row = [];
    let real;
    values.forEach((v, i) => {
        if (i % 2 == 0) {
            real = v;
        } else {
            row.push({
                real: real,
                imag: v
            });
            if ((i + 1) / 2 % columns == 0) {
                res.push(row);
                row = [];
            }
        }
    });
    return res;
}
216205
def rotvec2mat(u, phi):
    """"""

    phi = np.squeeze(phi)
    norm_u = np.linalg.norm(u)

    if norm_u < 1e-12:
        raise Exception("the rotation vector is equal to zero")

    u = u / norm_u
    # http://en.wikipedia.org/wiki/Rotation_matrix
    s = np.sin(phi)
    c = np.cos(phi)
    t = 1 - c

    ux = u[0]
    uy = u[1]
    uz = u[2]
    res = np.array([[t * ux * ux + c, t * ux * uy - s * uz, t * ux * uz + s * uy],
                    [t * ux * uy + s * uz, t * uy * uy + c, t * uy * uz - s * ux],
                    [t * ux * uz - s * uy, t * uy * uz + s * ux, t * uz * uz + c]])

    return res
241475
def d2Ibr_dV2(Ybr, V, lam):
    
    nb = len(V)
    diaginvVm = spdiag(div(matrix(1.0, (nb, 1)), abs(V)))

    Haa = spdiag(mul(-(Ybr.T * lam), V))
    Hva = -1j * Haa * diaginvVm
    Hav = Hva
    Hvv = spmatrix([], [], [], (nb, nb))

    return Haa, Hav, Hva, Hvv
581199
def _jmomentsurfaceMCIntegrand(vz,vR,vT,R,z,df,sigmaR1,gamma,sigmaz1,mvT,n,m,o):
    """"""
    return df(R,vR*sigmaR1,vT*sigmaR1*gamma,z,vz*sigmaz1,use_physical=False,
              func=(lambda x,y,z: x**n*y**m*z**o))\
              *numpy.exp(vR**2./2.+(vT-mvT)**2./2.+vz**2./2.)
1145532
def multidimensional_discretization(rho, sigma, N=3, method='rouwenhorst', m=2):
    

    # rho is assumed to be a scalar
    # sigma is a positive symmetric matrix
    # N number of points in each non-degenerate dimension
    # m : standard deviations to approximate

    import scipy.linalg
    from itertools import product

    d = sigma.shape[1]

    sigma = sigma.copy()

    zero_columns = np.where(sigma.sum(axis=0)==0)[0]
    for i in zero_columns:
        sigma[i,i] = 1

    L = scipy.linalg.cholesky(sigma)

    N = int(N)

    if method=='tauchen':
        [nodes_1d, probas_1d] = tauchen(N, 0, rho, 1, m=m)
    elif method=='rouwenhorst':
        [nodes_1d, probas_1d] = rouwenhorst(rho, 1, N)

    markov_nodes = np.array( list( product( *([nodes_1d]*d))) ).T
    markov_indices = np.array( list( product( *([range(N)]*d)  )  ) ).T

    markov_nodes = np.dot(L, markov_nodes)


    transition_matrix = 1
    for i in range(d):
        transition_matrix = np.kron(transition_matrix, probas_1d)

    markov_nodes = np.ascontiguousarray(markov_nodes.T)

    for i in zero_columns:
        markov_nodes[:,i] = 0

    return [markov_nodes, transition_matrix]
212810
def cublasZtrsm(handle, side, uplo, transa, diag, m, n, alpha, A, lda, B, ldb):
    
    
    status = _libcublas.cublasZtrsm_v2(handle, 
                                       _CUBLAS_SIDE_MODE[side], 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       _CUBLAS_OP[trans], 
                                       _CUBLAS_DIAG[diag], 
                                       m, n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,                    
                                                                               alpha.imag)),
                                       int(A), lda, int(B), ldb)
    cublasCheckStatus(status)
1357713
def cublasCgeru(handle, m, n, alpha, x, incx, y, incy, A, lda):
    

    status = _libcublas.cublasCgeru_v2(handle,
                                       m, n, ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                              alpha.imag)),
                                       int(x), incx, int(y), incy, int(A), lda)
    cublasCheckStatus(status)