Matrix-vector product for Hermitian-packed matrix.
6969
def cublasZhpmv(handle, uplo, n, alpha, AP, x, incx, beta, y, incy):
    
    
    status = _libcublas.cublasZhpmv_v2(handle,
                                       _CUBLAS_FILL_MODE[uplo], 
                                       n, ctypes.byref(cuda.cuDoubleComplex(alpha.real,
                                                                            alpha.imag)),
                                       int(AP), int(x), incx,
                                       ctypes.byref(cuda.cuDoubleComplex(beta.real, 
                                                                         beta.imag)),
                                       int(y), incy)
    cublasCheckStatus(status)
801886
def _M2_sparse_sym(Xvar, mask_X, Yvar, mask_Y, weights=None, column_selection=None):
    
    assert len(mask_X) == len(mask_Y), 'X and Y need to have equal sizes for symmetrization'

    if column_selection is None:
        mask_Xk = mask_X
        mask_Yk = mask_Y
        Xvark = Xvar
        Yvark = Yvar
    else:
        mask_Xk = mask_X[column_selection]
        mask_Yk = mask_Y[column_selection]
        Xvark = Xvar[:, _filter_variable_indices(mask_X, column_selection)]
        Yvark = Yvar[:, _filter_variable_indices(mask_Y, column_selection)]

    Cxxyy = np.zeros((len(mask_X), len(mask_Yk)))
    Cxxyy[np.ix_(mask_X, mask_Xk)] = _M2_dense(Xvar, Xvark, weights=weights)
    Cxxyy[np.ix_(mask_Y, mask_Yk)] += _M2_dense(Yvar, Yvark, weights=weights)

    Cxyyx = np.zeros((len(mask_X), len(mask_Yk)))
    Cxy = _M2_dense(Xvar, Yvark, weights=weights)
    Cyx = _M2_dense(Yvar, Xvark, weights=weights)
    Cxyyx[np.ix_(mask_X, mask_Yk)] = Cxy
    Cxyyx[np.ix_(mask_Y, mask_Xk)] += Cyx

    return Cxxyy, Cxyyx
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)
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)
1131398
def cross_product_matrix(vec):
    """"""
    return np.array([[0, -vec[2], vec[1]],
                     [vec[2], 0, -vec[0]],
                     [-vec[1], vec[0], 0]])
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);
        }
    }
1278575
def cublasZtrmm(handle, side, uplo, trans, diag, m, n, alpha, A, lda, B, ldb, C, ldc):
    
    
    status = _libcublas.cublasZtrmm_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, int(C), ldc)
    cublasCheckStatus(status)
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
500618
function svdJs() {
                var A = this;
                var V = Matrix$b.I(A.rows());
                var S = A.transpose();
                var U = Matrix$b.I(A.cols());
                var err = Number.MAX_VALUE;
                var i = 0;
                var maxLoop = 100;

                while(err > 2.2737e-13 && i < maxLoop) {
                    var qr = S.transpose().qrJs();
                    S = qr.R;
                    V = V.x(qr.Q);
                    qr = S.transpose().qrJs();
                    U = U.x(qr.Q);
                    S = qr.R;

                    var e = S.triu(1).unroll().norm();
                    var f = S.diagonal().norm();

                    if(f == 0)
                        f = 1;

                    err = e / f;

                    i++;
                }

                var ss = S.diagonal();
                var s = [];

                for(var i = 1; i <= ss.cols(); i++) {
                    var ssn = ss.e(i);
                    s.push(Math.abs(ssn));

                    if(ssn < 0) {
                        for(var j = 0; j < U.rows(); j++) {
                            V.elements[j][i - 1] = -(V.elements[j][i - 1]);
                        }
                    }
                }

                return {U: U, S: $V(s).toDiagonalMatrix(), V: V};
            }
408799
def tf_min_eig_vec(self):
    """"""
    # Full eigen decomposition requires the explicit psd matrix M
    _, matrix_m = self.dual_object.get_full_psd_matrix()
    [eig_vals, eig_vectors] = tf.self_adjoint_eig(matrix_m)
    index = tf.argmin(eig_vals)
    return tf.reshape(
        eig_vectors[:, index], shape=[eig_vectors.shape[0].value, 1])
669277
def _kvmatrix2d(km,vm):
    '''
        
        
    '''
    d = {}
    kmwfs = get_kmwfs(km)
    vmwfs = elel.get_wfs(vm)
    lngth = vmwfs.__len__()
    for i in range(0,lngth):
        value = elel.getitem_via_pathlist(vm,vmwfs[i])
        cond = elel.is_leaf(value)
        if(cond):
            _setitem_via_pathlist(d,kmwfs[i],value)
        else:
            _setdefault_via_pathlist(d,kmwfs[i])
    return(d)
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.)
576046
def cublasChpmv(handle, uplo, n, alpha, AP, x, incx, beta, y, incy):
    
    
    status = _libcublas.cublasChpmv_v2(handle, 
                                       _CUBLAS_FILL_MODE[uplo], 
                                       n, ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                           alpha.imag)),
                                       int(AP), int(x), incx,
                                       ctypes.byref(cuda.cuFloatComplex(beta.real,
                                                                        beta.imag)),
                                       int(y), incy)
    cublasCheckStatus(status)
449262
def _gill_king(mat, beta, delta):
    """"""
    size = mat.shape[0]

    # initialize d_vec and lowtri
    if scipy.sparse.issparse(mat):
        lowtri = scipy.sparse.eye(*mat.shape)
    else:
        lowtri = numpy.eye(size)

    d_vec = numpy.zeros(size, dtype=float)

    # there are no inner for loops, everything implemented with
    # vector operations for a reasonable level of efficiency
    for idx in range(size):
        if idx == 0:
            idz = []     # column index: all columns to left of diagonal
                        # d_vec(idz) doesn't work in case idz is empty
        else:
            idz = numpy.s_[:idx]

        djtemp = mat[idx, idx] - numpy.dot(
            lowtri[idx, idz], d_vec[idz]*lowtri[idx, idz].T)
        # C(idx, idx) in book

        if idx < size - 1:
            idy = numpy.s_[idx+1:size]
            # row index: all rows below diagonal
            ccol = mat[idy, idx] - numpy.dot(
                lowtri[idy, idz], d_vec[idz]*lowtri[idx, idz].T)
            # C(idy, idx) in book
            theta = abs(ccol).max()
            # guarantees d_vec(idx) not too small and lowtri(idy, idx) not too
            # big in sufficiently positive definite case, d_vec(idx) = djtemp
            d_vec[idx] = max(abs(djtemp), (theta/beta)**2, delta)
            lowtri[idy, idx] = ccol/d_vec[idx]

        else:
            d_vec[idx] = max(abs(djtemp), delta)

    # convert to usual output format: replace lowtri by lowtri*sqrt(D) and
    # transpose
    for idx in range(size):
        lowtri[:, idx] = lowtri[:, idx]*numpy.sqrt(d_vec[idx])
        # lowtri = lowtri*diag(sqrt(d_vec)) bad in sparse case

    return lowtri
124553
def cublasCsyr(handle, uplo, n, alpha, x, incx, A, lda):
    

    status = _libcublas.cublasCsyr_v2(handle,
                                        _CUBLAS_FILL_MODE[uplo], n, 
                                        ctypes.byref(cuda.cuFloatComplex(alpha.real,
                                                                        alpha.imag)),
                                        int(x), incx, int(A), lda)
    cublasCheckStatus(status)
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)
1028992
def numbaGaussian2d(psf, sy, sx):
    '''
    
    '''
    ps0, ps1 = psf.shape
    c0,c1 = ps0//2, ps1//2
    ssx = 2*sx**2
    ssy = 2*sy**2
    for i in range(ps0):
        for j in range(ps1):
            psf[i,j]=exp( -( (i-c0)**2/ssy
                            +(j-c1)**2/ssx) )
    psf/=psf.sum()