關于上一節(jié)讀者某些疑問：為什么你用進程并行不是線程并行？

回答：由于Python解釋器有GIL（全局解釋器鎖），在單進程的解釋器上有線程安全鎖，也就是說每次只能一個線程訪問解釋器，因此Python在語法上的多線程（multithreads）實現是不會提高并行性能的。
這一點和C\C++上的編譯級別的并行是不一樣的，Python能做到的極限是多進程的解釋級別并行。（上一節(jié)我實現的是多進程并行，和老師課上MPI的多線程是不一樣的?。。?/p>

0. 引言

OpenMP是一種C/C++的并行編譯標準方案，嚴謹地說，在Python上使用OpenMP是不可能的，因為本來Python是一門解釋語言。

但如果在某個解釋流程實現并行優(yōu)化是可行的，也有一些方案。

如下所示C和Python在實現OpenMP接口的可能性對比：

如上兩個地方可以實現mp的優(yōu)化，分別為：在可交互級別的轉化，在解釋階段級別的轉化。

• 在可交互級別的轉化：一些第三方包pymp，pyopenmp，使用fork脫離GIL

• 在解釋階段級別的轉化：使用Cython直接編寫C code，這樣寫出來的module在解釋的時候會被解釋器做并行優(yōu)化。

早期Cython是不支持openmp的，但如今已支持了，我們可以非常方便地使用Cython語法編寫支持openmp并行優(yōu)化編譯的代碼。

因此，本文通過第一種實現方式對上一節(jié)的矩陣乘法優(yōu)化做探究。

生成示例矩陣向量、計算代碼都使用上一節(jié)內容。

鏈接：高性能計算的矩陣乘法優(yōu)化 - Python +MPI的實現

import numpy as np
from functools import wraps
import time

def generate_example_matrix(h, w):
    _vs = (-1, 2, -1)
    _i = -1  # shift bits
    example_data = np.zeros([h, w], dtype=np.int32)
    for i in range(3):
        example_data_eye_mask = np.eye(h, w, _i + i,
                                       dtype=np.bool_)  # build eyes and shift
        example_data[example_data_eye_mask == True] = _vs[i]
    return example_data


def generate_example_vector(w):
    _rest_dict = {
        1: [1],
        2: [1, 2],
        3: [1, 2, 3],
    }
    rest_bits = int(w % 3)
    repeat_times = int(w // 3)
    example_vector = np.repeat([[1, 2, 3]], repeat_times, axis=0).ravel()

    if rest_bits > 0:
        tail_vec = _rest_dict[rest_bits]
        tail_vec = np.array(tail_vec, dtype=np.int32)
        example_vector = np.concatenate([example_vector, tail_vec], axis=0)
    return example_vector
    

計算的naive code如下：

def naive_method(example_matrix, example_vector):
    result = []
    h, w = example_matrix.shape
    for hi in range(h):
        colv = example_matrix[hi, :]
        temp = 0
        for wi in range(w):
            temp += colv[wi] * example_vector[wi]
        result.append(temp)
    return np.array(result)

單進程單線程執(zhí)行：

from utils import generate_example_matrix, generate_example_vector
import pymp
import time
import numpy as np

def naive_method(example_matrix, example_vector):
    result = []
    h, w = example_matrix.shape
    for hi in range(h):
        colv = example_matrix[hi, :]
        temp = 0
        for wi in range(w):
            temp += colv[wi] * example_vector[wi]
        result.append(temp)
    return np.array(result)

def main():
    start_time = time.perf_counter()
    h = 5000
    w = 5000
    print('--- Current matrix scale is: {} x {} ---'.format(h,w))
    example_matrix = generate_example_matrix(h, w)
    example_vector = generate_example_vector(w)
    result = naive_method(example_matrix, example_vector)
    end_time = time.perf_counter()
    print('single method used time is: {:.2f}s\n'.format(end_time - start_time))

if __name__ == '__main__':
    main()

1. 在可交互級別的轉化

一個實現上是openmp stye的第三方庫：pymp-pypi.

注意，該方法不能在Windows上實現，因為他使用fork來繞過GIL的。

操作系統的fork方法可以繞過GIL（全局解釋器鎖），從而實現多線程，這比Python原生的multithreads會更加并行高效。

題外話：當然還有一堆其他類似的第三方庫，看看別人的教程也可以輕松實現。我這里只挑一個

安裝：pip install pymp-pypi

基本用法：

優(yōu)化前：

ex_array = np.zeros((100,), dtype='uint8')

for index in range(0, 100):
    ex_array[index] = 1
    print('Yay! {} done!'.format(index))
    

優(yōu)化后：

import pymp

ex_array = pymp.shared.array((100,), dtype='uint8')

with pymp.Parallel(4) as p:
    for index in p.range(0, 100):
        ex_array[index] = 1
        # The parallel print function takes care of asynchronous output.
        p.print('Yay! {} done!'.format(index))

1.1 OpenMP Style的改寫

將該用法改寫到我們原來的運算代碼中：

優(yōu)化后：

from utils import generate_example_matrix, generate_example_vector
import pymp
import time
import numpy as np

def naive_multi_method(example_matrix, example_vector, threads):
    # result = pymp.shared.list()
    result = []
    h, w = example_matrix.shape
    with pymp.Parallel(num_threads = threads) as p:
        for hi in p.range(0, h):
            colv = example_matrix[hi, :]
            temp = 0
            for wi in p.range(0, w):
                temp += colv[wi] * example_vector[wi]
            result.append(temp)

def main():
    start_time = time.perf_counter()
    h = 50000
    w = 50000
    threads_num = 200 # 250, 200, 100, 50, 25, 16, 8
    print('--- Current matrix scale is: {} x {} ---'.format(h,w))
    print('<- Threads num is: {} ->'.format(threads_num))
    example_matrix = generate_example_matrix(h, w)
    example_vector = generate_example_vector(w)
    result = naive_multi_method(example_matrix, example_vector, threads = threads_num)
    end_time = time.perf_counter()
    print('multi-thread method used time is: {:.2f}s\n'.format(end_time - start_time))
    return np.array(result)
if __name__ == '__main__':
    main()

2.2 實驗對比：

相較于上一次的實驗，我換了一臺設備，CPU配置如下：

11th Gen Intel? Core? i5-11600K @ 3.90GHz，此外，我超頻到了4.5GHz

這一次直接測試50000規(guī)模的運算。

單線程運行：

--- Current matrix scale is: 50000 x 50000 ---
single method used time is: 475.64s

Baseline：475秒

多線程優(yōu)化之后：

8線程：

--- Current matrix scale is: 50000 x 50000 ---
<- Threads num is: 8 ->
multi-thread method used time is: 19.46s

T8：19秒

16線程：13秒

--- Current matrix scale is: 50000 x 50000 ---
<- Threads num is: 16 ->
multi-thread method used time is: 13.28s

T16：13秒

25線程：

--- Current matrix scale is: 50000 x 50000 ---
<- Threads num is: 25 ->
multi-thread method used time is: 11.12s

T25：11秒

50線程：

--- Current matrix scale is: 50000 x 50000 ---
<- Threads num is: 50 ->
multi-thread method used time is: 9.18s

T50：9.18秒

100線程：

--- Current matrix scale is: 50000 x 50000 ---
<- Threads num is: 100 ->
multi-thread method used time is: 8.27s

T100：8.27秒

200線程：

--- Current matrix scale is: 50000 x 50000 ---
<- Threads num is: 200 ->
multi-thread method used time is: 8.23s

T200：8.23秒

250線程：

--- Current matrix scale is: 50000 x 50000 ---
<- Threads num is: 250 ->
multi-thread method used time is: 8.47s

T250：8.47秒

從以上數據可見，線程瓶頸數在200到250之間，接下來可以使用二分查找測試出最佳的線程數文章來源地址http://www.zghlxwxcb.cn/news/detail-553868.html

到了這里，關于高性能計算的矩陣乘法優(yōu)化 - Python + OpenMP實現的文章就介紹完了。如果您還想了解更多內容，請在右上角搜索TOY模板網以前的文章或繼續(xù)瀏覽下面的相關文章，希望大家以后多多支持TOY模板網！