Lapack的DSYSV

问题描述

我在Cython中编写了一种算法，该算法经常执行Lapack和Blas例程。由于对这些例程的单独调用彼此独立，因此我尝试通过使用OpenMP并行化计算来加快处理速度。不幸的是，与单线程版本相比，使用多线程会大大降低运行速度，而且我不清楚发生了什么。

最小的工作示例

我为您创建了一个最小的工作示例，以重现该问题。

在文件main.py中创建一些伪数据，并启动在Cython函数run_test中实现的算法：

import numpy as np
import cProfile
from test import run_test

if __name__ == '__main__':
    num_ordinates = 500
    np.random.seed(1)
    m = np.ascontiguousarray(np.where(np.random.randint(2,size=(num_ordinates,)) > 0,1,-1),dtype=float)
    ordinates = np.zeros(num_ordinates,dtype=float,order='C')
    num_coefficients = (num_ordinates * (num_ordinates + 1) // 2)
    coefficients = np.zeros(num_coefficients,order='C')
    k = 0

    for i in range(num_ordinates):
        v = m[i]
        ordinates[i] += v / (num_ordinates + 1)

        for j in range(i):
            coefficients[k] -= -m[j] * v / pow(num_ordinates + 1,2)
            k += 1

        coefficients[k] += abs(v) * num_ordinates / pow(num_ordinates + 1,2)
        k += 1

    num_threads = 8
    cProfile.run('run_test(coefficients,ordinates,num_threads)')

该算法在名为test.pyx的文件中实现。它对给定的数据执行相同计算的多次迭代。在每次迭代中，都会调用Lapack的DSYSV例程以及Blas的DDOT和DSPMV例程。外循环使用prange命令能够在指定数量的线程之间分布迭代。

from libc.stdlib cimport malloc,free
from scipy.linalg.cython_lapack cimport dsysv
from scipy.linalg.cython_blas cimport ddot,dspmv
from cython.parallel cimport prange

cpdef double run_test(double[::1] coefficients,double[::1] ordinates,int num_threads):
    cdef double *coefficients_ptr = &coefficients[0]
    cdef double *ordinates_ptr = &ordinates[0]
    cdef int n = ordinates.shape[0]
    cdef double result = 0
    cdef double *out
    cdef double *out2
    cdef double tmp
    cdef int i

    for i in prange(1000,nogil=True,num_threads=num_threads):
        out = run_dsysv2(coefficients_ptr,ordinates_ptr,n)
        tmp = run_ddot(out,n)
        out2 = run_dspmv(coefficients_ptr,out,n)
        result += run_ddot(out,out2,n)
        free(out)
        free(out2)

    return result

cdef double* run_dsysv(double *coefficients,double *ordinates,int n) nogil except +:
    cdef double *a = <double*>malloc(n * n * sizeof(double))
    cdef double *b = <double*>malloc(n * sizeof(double))
    cdef int i,j,offset
    cdef int k = 0

    for i in range(n):
        b[i] = ordinates[i]
        offset = i * n

        for j in range(i + 1):
            a[offset + j] = coefficients[k]
            k += 1

    cdef char* uplo = 'U'
    cdef int nrhs = 1
    cdef int info
    cdef double worksize
    cdef int lwork = -1
    dsysv(uplo,&n,&nrhs,a,<int*>0,b,&worksize,&lwork,&info)
    lwork = <int>worksize
    cdef double *work = <double*>malloc(lwork * sizeof(double))
    cdef int *ipiv = <int*>malloc(n * sizeof(int))
    dsysv(uplo,ipiv,work,&info)

    free(ipiv)
    free(work)
    free(a)

    if info != 0:
        with gil:
            raise ArithmeticError('Non-zero info code: ' + str(info))

    return b

cdef double run_ddot(double *x,double *y,int n) nogil:
    cdef int inc = 1
    return ddot(&n,x,&inc,y,&inc)

cdef double* run_dspmv(double* m,double* a,int n) nogil:
    cdef char* uplo = 'U'
    cdef double alpha = 1
    cdef int inc = 1
    cdef double beta = 0
    cdef double *out = <double*>malloc(n * sizeof(double))
    dspmv(uplo,&alpha,m,&beta,&inc)
    return out

为完整起见，这也是setup.py文件：

import numpy
from Cython.Build import cythonize
from setuptools import setup,Extension

extensions = [
    Extension(name='*',sources=['*.pyx'],language='c++',extra_compile_args=['-fopenmp'],extra_link_args=['-fopenmp'],define_macros=[("NPY_NO_DEPRECATED_API","NPY_1_7_API_VERSION")])
]

compiler_directives = {
    'boundscheck': False,'wraparound': False,'cdivision': True,'initializedcheck': False
}

setup(name='test',packages=['test'],install_requires=[
          "numpy>=1.19.0","scipy>=1.5.0","Cython>=0.29.0"
      ],python_requires='>=3.7',ext_modules=cythonize(extensions,language_level='3',compiler_directives=compiler_directives),include_dirs=[numpy.get_include()])

我使用命令python3.7 setup.py build_ext --inplace编译程序。

分析

通过在num_threads中将变量main.py分别设置为1和8，我在系统上运行了单线程和多线程代码变体。我获得了以下运行时：

num_threads=1: 5.614 seconds
num_threads=8: 50.541 seconds

如您所见，多线程变体的速度慢了大约十倍，而且我不清楚为什么会这样。

这是我系统的规格：

OS: (Arch) Linux
16GB RAM
Intel Core i7-4710MQ CPU (8 cores @ 2.50GHz)

我想知道正在使用哪个版本的Blas和Lapack也将有所帮助。我正在使用scipy 1.5.2打包的版本，但是我不知道确切的版本号。我将不胜感激有关如何了解它们的任何提示。

一些其他测试

我认为可能是由于运行DSYSV例程之前分配了过多的内存而导致出现此问题，从而导致某些其他数据被从某些CPU缓存中推出。为了对此进行调查，我尝试从run_dsysv中的函数test.pyx中删除对DSYSV例程的调用。我用以下代码替换了原始函数，在其中我尝试保持与原始代码的距离尽可能近（特别是，我想保留内存分配）：

cdef double* run_dsysv(double *coefficients,&info)
    lwork = <int>worksize
    cdef double *work = <double*>malloc(lwork * sizeof(double))
    cdef int *ipiv = <int*>malloc(n * sizeof(int))

    # dsysv(uplo,&info)
    cdef double v

    for i in range(n):
        k = (((i + 1) * (i + 2)) // 2) - 1
        v = coefficients[k]
        v = -ordinates[i] / v if v != 0 else 0
        b[i] = v

    free(ipiv)
    free(work)
    free(a)

    return b

运行更新后的代码的单线程和多线程版本会导致以下结果：

num_threads=1: 0.160 seconds
num_threads=8: 0.089 seconds

该程序当然比以前快得多，因为对计算量要求很高的DSYSV例程已被一些线性运行的代码所代替。尽管如此，我们可以看到，与单线程版本相比，现在使用多个线程确实可以提高速度。至少它不会导致明显的减速。因此，我认为问题与DSYSV例程有关，但不一定与为其分配的内存有关。

解决方法

暂无找到可以解决该程序问题的有效方法，小编努力寻找整理中！

如果你已经找到好的解决方法，欢迎将解决方案带上本链接一起发送给小编。

小编邮箱:dio#foxmail.com (将#修改为@）