Python linalg.get_blas_funcs函数代码示例

    def test_spr_hpr(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES+COMPLEX_DTYPES):
            n = 3
            A = rand(n, n).astype(dtype)
            if ind > 1:
                A += rand(n, n)*1j
            A = A.astype(dtype)
            A = A + A.T if ind < 4 else A + A.conj().T
            c, r = tril_indices(n)
            Ap = A[r, c]
            x = rand(n).astype(dtype)
            alpha = (DTYPES+COMPLEX_DTYPES)[mod(ind, 4)](2.5)

            if ind > 3:
                func, = get_blas_funcs(('hpr',), dtype=dtype)
                y2 = alpha * x[:, None].dot(x[None, :].conj()) + A
            else:
                func, = get_blas_funcs(('spr',), dtype=dtype)
                y2 = alpha * x[:, None].dot(x[None, :]) + A

            y1 = func(n=n, alpha=alpha, ap=Ap, x=x)
            y1f = zeros((3, 3), dtype=dtype)
            y1f[r, c] = y1
            y1f[c, r] = y1.conj() if ind > 3 else y1
            assert_array_almost_equal(y1f, y2)

    def test_spmv_hpmv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES+COMPLEX_DTYPES):
            n = 3
            A = rand(n, n).astype(dtype)
            if ind > 1:
                A += rand(n, n)*1j
            A = A.astype(dtype)
            A = A + A.T if ind < 4 else A + A.conj().T
            c, r = tril_indices(n)
            Ap = A[r, c]
            x = rand(n).astype(dtype)
            y = rand(n).astype(dtype)
            xlong = arange(2*n).astype(dtype)
            ylong = ones(2*n).astype(dtype)
            alpha, beta = dtype(1.25), dtype(2)

            if ind > 3:
                func, = get_blas_funcs(('hpmv',), dtype=dtype)
            else:
                func, = get_blas_funcs(('spmv',), dtype=dtype)
            y1 = func(n=n, alpha=alpha, ap=Ap, x=x, y=y, beta=beta)
            y2 = alpha * A.dot(x) + beta * y
            assert_array_almost_equal(y1, y2)

            # Test inc and offsets
            y1 = func(n=n-1, alpha=alpha, beta=beta, x=xlong, y=ylong, ap=Ap,
                      incx=2, incy=2, offx=n, offy=n)
            y2 = (alpha * A[:-1, :-1]).dot(xlong[3::2]) + beta * ylong[3::2]
            assert_array_almost_equal(y1[3::2], y2)
            assert_almost_equal(y1[4], ylong[4])

def _have_blas_gemm():
    try:
        linalg.get_blas_funcs(["gemm"])
        return True
    except (AttributeError, ValueError):
        warnings.warn("Could not import BLAS, falling back to np.dot")
        return False

def test_get_blas_funcs():
    # check that it returns Fortran code for arrays that are
    # fortran-ordered
    f1, f2, f3 = get_blas_funcs(
        ('axpy', 'axpy', 'axpy'),
        (np.empty((2,2), dtype=np.complex64, order='F'),
         np.empty((2,2), dtype=np.complex128, order='C'))
        )

    # get_blas_funcs will choose libraries depending on most generic
    # array
    assert_equal(f1.typecode, 'z')
    assert_equal(f1.module_name, 'cblas')
    assert_equal(f2.typecode, 'z')
    assert_equal(f2.module_name, 'cblas')

    # check defaults.
    f1 = get_blas_funcs('rotg')
    assert_equal(f1.typecode, 'd')

    # check also dtype interface
    f1 = get_blas_funcs('gemm', dtype=np.complex64)
    assert_equal(f1.typecode, 'c')
    f1 = get_blas_funcs('gemm', dtype='F')
    assert_equal(f1.typecode, 'c')

    def test_sbmv_hbmv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 6
            k = 2
            A = zeros((n, n), dtype=dtype)
            Ab = zeros((k+1, n), dtype=dtype)

            # Form the array and its packed banded storage
            A[arange(n), arange(n)] = rand(n)
            for ind2 in range(1, k+1):
                temp = rand(n-ind2)
                A[arange(n-ind2), arange(ind2, n)] = temp
                Ab[-1-ind2, ind2:] = temp
            A = A.astype(dtype)
            A = A + A.T if ind < 2 else A + A.conj().T
            Ab[-1, :] = diag(A)
            x = rand(n).astype(dtype)
            y = rand(n).astype(dtype)
            alpha, beta = dtype(1.25), dtype(3)

            if ind > 1:
                func, = get_blas_funcs(('hbmv',), dtype=dtype)
            else:
                func, = get_blas_funcs(('sbmv',), dtype=dtype)
            y1 = func(k=k, alpha=alpha, a=Ab, x=x, y=y, beta=beta)
            y2 = alpha * A.dot(x) + beta * y
            assert_array_almost_equal(y1, y2)

    def test_spr2_hpr2(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 3
            A = rand(n, n).astype(dtype)
            if ind > 1:
                A += rand(n, n)*1j
            A = A.astype(dtype)
            A = A + A.T if ind < 2 else A + A.conj().T
            c, r = tril_indices(n)
            Ap = A[r, c]
            x = rand(n).astype(dtype)
            y = rand(n).astype(dtype)
            alpha = dtype(2)

            if ind > 1:
                func, = get_blas_funcs(('hpr2',), dtype=dtype)
            else:
                func, = get_blas_funcs(('spr2',), dtype=dtype)

            u = alpha.conj() * x[:, None].dot(y[None, :].conj())
            y2 = A + u + u.conj().T
            y1 = func(n=n, alpha=alpha, x=x, y=y, ap=Ap)
            y1f = zeros((3, 3), dtype=dtype)
            y1f[r, c] = y1
            y1f[[1, 2, 2], [0, 0, 1]] = y1[[1, 3, 4]].conj()
            assert_array_almost_equal(y1f, y2)

def test_get_blas_funcs_alias():
    # check alias for get_blas_funcs
    f, g = get_blas_funcs(('nrm2', 'dot'), dtype=np.complex64)
    assert f.typecode == 'c'
    assert g.typecode == 'c'

    f, g, h = get_blas_funcs(('dot', 'dotc', 'dotu'), dtype=np.float64)
    assert f is g
    assert f is h

def test_fast_dot():
    """Check fast dot blas wrapper function"""
    rng = np.random.RandomState(42)
    A = rng.random_sample([2, 10])
    B = rng.random_sample([2, 10])

    try:
        linalg.get_blas_funcs('gemm')
        has_blas = True
    except AttributeError, ValueError:
        has_blas = False

def py_rect_maxvol(A, tol = 1., maxK = None, min_add_K = None, minK = None, start_maxvol_iters = 10, identity_submatrix = True):
    """Python implementation of rectangular 2-volume maximization. For information see :py:func:`rect_maxvol` function"""
    # tol2 - square of parameter tol
    tol2 = tol**2
    # N - number of rows, r - number of columns of matrix A
    N, r = A.shape
    if N <= r:
        return np.arange(N, dtype = np.int32), np.eye(N, dtype = A.dtype)
    if maxK is None or maxK > N:
        maxK = N
    if maxK < r:
        maxK = r
    if minK is None or minK < r:
        minK = r
    if minK > N:
        minK = N
    if min_add_K is not None:
        minK = max(minK, r + min_add_K) 
    if minK > maxK:
        minK = maxK
        #raise ValueError('minK value cannot be greater than maxK value')
    index = np.zeros(N, dtype = np.int32)
    chosen = np.ones(N)
    tmp_index, C = py_maxvol(A, tol = 1, max_iters = start_maxvol_iters)
    index[:r] = tmp_index
    chosen[tmp_index] = 0
    C = np.asfortranarray(C)
    # compute square 2-norms of each row in matrix C
    row_norm_sqr = np.array([chosen[i]*np.linalg.norm(C[i], 2)**2 for i in xrange(N)])
    # find maximum value in row_norm_sqr
    i = np.argmax(row_norm_sqr)
    K = r
    # set cgeru or zgeru for complex numbers and dger or sger for float numbers
    try:
        ger = get_blas_funcs('geru', [C])
    except:
        ger = get_blas_funcs('ger', [C])
    while (row_norm_sqr[i] > tol2 and K < maxK) or K < minK:
        # add i to index and recompute C and square norms of each row by SVM-formula
        index[K] = i
        chosen[i] = 0
        c = C[i].copy()
        v = C.dot(c.conj())
        l = 1.0/(1+v[i])
        ger(-l,v,c,a=C,overwrite_a=1)
        C = np.hstack([C, l*v.reshape(-1,1)])
        row_norm_sqr -= (l*v*v.conj()).real
        row_norm_sqr *= chosen
        # find maximum value in row_norm_sqr
        i = row_norm_sqr.argmax()
        K += 1
    if identity_submatrix:
        C[index[:K]] = np.eye(K, dtype = C.dtype)
    return index[:K].copy(), C

def test_inplace_swap_column():
    X = np.array([[0, 3, 0],
                  [2, 4, 0],
                  [0, 0, 0],
                  [9, 8, 7],
                  [4, 0, 5]], dtype=np.float64)
    X_csr = sp.csr_matrix(X)
    X_csc = sp.csc_matrix(X)

    swap = linalg.get_blas_funcs(('swap',), (X,))
    swap = swap[0]
    X[:, 0], X[:, -1] = swap(X[:, 0], X[:, -1])
    inplace_swap_column(X_csr, 0, -1)
    inplace_swap_column(X_csc, 0, -1)
    assert_array_equal(X_csr.toarray(), X_csc.toarray())
    assert_array_equal(X, X_csc.toarray())
    assert_array_equal(X, X_csr.toarray())

    X[:, 0], X[:, 1] = swap(X[:, 0], X[:, 1])
    inplace_swap_column(X_csr, 0, 1)
    inplace_swap_column(X_csc, 0, 1)
    assert_array_equal(X_csr.toarray(), X_csc.toarray())
    assert_array_equal(X, X_csc.toarray())
    assert_array_equal(X, X_csr.toarray())
    assert_raises(TypeError, inplace_swap_column, X_csr.tolil())

    X = np.array([[0, 3, 0],
                  [2, 4, 0],
                  [0, 0, 0],
                  [9, 8, 7],
                  [4, 0, 5]], dtype=np.float32)
    X_csr = sp.csr_matrix(X)
    X_csc = sp.csc_matrix(X)
    swap = linalg.get_blas_funcs(('swap',), (X,))
    swap = swap[0]
    X[:, 0], X[:, -1] = swap(X[:, 0], X[:, -1])
    inplace_swap_column(X_csr, 0, -1)
    inplace_swap_column(X_csc, 0, -1)
    assert_array_equal(X_csr.toarray(), X_csc.toarray())
    assert_array_equal(X, X_csc.toarray())
    assert_array_equal(X, X_csr.toarray())
    X[:, 0], X[:, 1] = swap(X[:, 0], X[:, 1])
    inplace_swap_column(X_csr, 0, 1)
    inplace_swap_column(X_csc, 0, 1)
    assert_array_equal(X_csr.toarray(), X_csc.toarray())
    assert_array_equal(X, X_csc.toarray())
    assert_array_equal(X, X_csr.toarray())
    assert_raises(TypeError, inplace_swap_column, X_csr.tolil())

def norm2(q):
    """
    Compute the euclidean norm of an array ``q`` by calling the BLAS routine
    """
    q = np.asarray(q)
    nrm2 = get_blas_funcs('nrm2', dtype=q.dtype)
    return nrm2(q)

    def test_trsv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 15
            A = (rand(n, n)+eye(n)).astype(dtype)
            x = rand(n).astype(dtype)
            func, = get_blas_funcs(('trsv',), dtype=dtype)

            y1 = func(a=A, x=x)
            y2 = solve(triu(A), x)
            assert_array_almost_equal(y1, y2)

            y1 = func(a=A, x=x, lower=1)
            y2 = solve(tril(A), x)
            assert_array_almost_equal(y1, y2)

            y1 = func(a=A, x=x, diag=1)
            A[arange(n), arange(n)] = dtype(1)
            y2 = solve(triu(A), x)
            assert_array_almost_equal(y1, y2)

            y1 = func(a=A, x=x, diag=1, trans=1)
            y2 = solve(triu(A).T, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(a=A, x=x, diag=1, trans=2)
            y2 = solve(triu(A).conj().T, x)
            assert_array_almost_equal(y1, y2)

    def test_tpsv(self):
        seed(1234)
        for ind, dtype in enumerate(DTYPES):
            n = 10
            x = rand(n).astype(dtype)
            # Upper triangular array
            A = triu(rand(n, n)) if ind < 2 else triu(rand(n, n)+rand(n, n)*1j)
            A += eye(n)
            # Form the packed storage
            c, r = tril_indices(n)
            Ap = A[r, c]
            func, = get_blas_funcs(('tpsv',), dtype=dtype)

            y1 = func(n=n, ap=Ap, x=x)
            y2 = solve(A, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1)
            A[arange(n), arange(n)] = dtype(1)
            y2 = solve(A, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1, trans=1)
            y2 = solve(A.T, x)
            assert_array_almost_equal(y1, y2)

            y1 = func(n=n, ap=Ap, x=x, diag=1, trans=2)
            y2 = solve(A.conj().T, x)
            assert_array_almost_equal(y1, y2)

def _update_dict_slow(X, A, B, G, Q, Q_idx, idx, fit_intercept,
                      components_range, norm, impute=True):
    Q_idx = Q[:, idx]

    if impute:
        old_sub_G = Q_idx.dot(Q_idx.T)

    ger, = linalg.get_blas_funcs(('ger',), (A, Q_idx))
    R = B[:, idx] - np.dot(Q_idx.T, A).T

    # norm = np.sqrt(np.sum(Q_idx ** 2, axis=1))
    norm = np.sqrt(np.sum(Q ** 2, axis=1))
    # Intercept on first column
    for j in components_range:
        ger(1.0, A[j], Q_idx[j], a=R, overwrite_a=True)
        Q_idx[j] = R[j] / A[j, j]
        # new_norm = np.sqrt(np.sum(Q_idx[j] ** 2))
        # if new_norm > norm[j]:
        #     Q_idx[j] /= new_norm / norm[j]
        Q[j, idx] = Q_idx[j]
        new_norm = np.sqrt(np.sum(Q[j] ** 2))
        if new_norm > 1:
            Q_idx[j] /= new_norm
            Q[j] /= new_norm

        ger(-1.0, A[j], Q_idx[j], a=R, overwrite_a=True)

    Q[:, idx] = Q_idx

    if impute:
        G += Q_idx.dot(Q_idx.T) - old_sub_G

def cool_syrk(fact, X): 
    syrk = get_blas_funcs("syrk", [X])
    R = syrk(fact, X)
    d = np.diag(R).copy()
    size = mat_to_upper_F(R)
    R.resize([size,])
    return R,d