MKL разреженный матричный продукт не имеет смысла - PullRequest
Новая
       41

MKL разреженный матричный продукт не имеет смысла

1 голос
/ 05 ноября 2019

Я составил следующий код из фрагментов, которые я нашел в Интернете.

Пока код выполняется, но результаты не имеют смысла.

У меня три вопроса:

a) Как мне составить выходную матрицу C, еслитолько данные матрицы возвращаются по ссылке?

b) Почему данные не похожи на результаты scipy (которые являются правильными)

c) Почему код MKL намного медленнее? 

import numpy as np
from time import time
from scipy.sparse import csc_matrix, random, diags
from ctypes import *

mkl = cdll.LoadLibrary("libmkl_rt.so")


def sparse_csc_mm(A, B):
    # from: https://groups.google.com/a/continuum.io/forum/#!topic/anaconda/TpuMdSVSf4E
    # A must be in CSC format, could enforce here with a check

    m = A.shape[0]
    n = A.shape[1]
    Bn = B.shape[1]

    Ax = A.data.ctypes.data_as(POINTER(c_double))
    Ai = A.indices.ctypes.data_as(POINTER(c_int))
    Ap =      A.indptr.ctypes.data_as(POINTER(c_int))
    eindptr = A.indptr.ctypes.data_as(POINTER(c_int))
    eindptr_ = cast(pointer(eindptr), POINTER(c_void_p))
    eindptr_.contents.value += sizeof(eindptr._type_)

    C = np.zeros(m, dtype=np.float64)

    Bx = B.data.ctypes.data_as(POINTER(c_double))
    Cx = C.ctypes.data_as(POINTER(c_double))

    alpha_beta = byref(c_double(1.))
    ldb_ldc = byref(c_int(0))
    no_transpose = byref(c_char(b'N'))
    general_matrix_c_order = c_char_p(b'G**C')

    #  mkl_dcscmm(transa, m, n, k,
    #             alpha, matdescra, val, indx, pntrb, pntre,
    #             b, ldb, beta, c, ldc)
    # docs: https://software.intel.com/en-us/mkl-developer-reference-fortran-mkl-cscmm
    ret = mkl.mkl_dcscmm(no_transpose, byref(c_int(m)), byref(c_int(Bn)), byref(c_int(n)),
                         alpha_beta, general_matrix_c_order, Ax, Ai, Ap, eindptr,
                         Bx, ldb_ldc, alpha_beta, Cx, ldb_ldc)

    return C


def test_csc_mm():
    np.random.seed(0)
    k = 10
    m, n = k, k

    A = csc_matrix(random(m, n, density=0.01)) + diags(np.ones(n))
    B = csc_matrix(random(m, n, density=0.01)) + diags(np.ones(n))

    t1 = time()
    C1 = A * B
    t_scipy = time() - t1
    print('time in scipy', t_scipy)

    t1 = time()
    C2 = sparse_csc_mm(A, B)
    t_mkl = time() - t1
    print('time in mkl',t_mkl)

    print('data in scipy', C1.data)
    print('data in mkl', C2)


if __name__ == '__main__':

    test_csc_mm()

Вывод на консоль выглядит следующим образом:

время в scipy 8.702278137207031e-05

время в mkl 0.021214723587036133

данные в scipy [1. 0.96896177 1. 0.09609841 1. 1. 1. 1. 1. 1. 1. 1.]

данные в мкл [10.09609841 10.09609841 9.78273334 10.09609841 10.09609841 10.09609841 10.09609841 10.09609841 10.09609841 10.09609841]

1024 *
...