请注意，下面的所有数据都是 float32，而不是双精度。
UPDATE：
我完全放弃了 gcc，转而使用 Intel 的 icc。当性能至关重要并且您没有时间微调“编译器提示”以强制执行 gcc 向量化时（请参阅例如here http://locklessinc.com/articles/vectorize/)
log_omp.c,
GCC：gcc -o log_omp.so -fopenmp log_omp.c -lm -O3 -fPIC -shared -std=c99
ICC：icc -o log_omp.so -openmp loge_omp.c -lm -O3 -fPIC -shared -std=c99 -vec-report1 -xAVX -I/opt/intel/composer/mkl/include
#include <math.h>
#include "omp.h"
#include "mkl_vml.h"
#define restrict __restrict
inline void log_omp(int m, float * restrict a, float * restrict c);
void log_omp(int m, float * restrict a, float * restrict c)
{
int i;
#pragma omp parallel for default(none) shared(m,a,c) private(i)
for (i=0; i<m; i++) {
a[i] = log(c[i]);
}
}
// VML / icc only:
void log_VML(int m, float * restrict a, float * restrict c)
{
int i;
int split_to = 14;
int iter = m / split_to;
int additional = m % split_to;
// vsLn(m, c, a);
#pragma omp parallel for default(none) shared(m,a,c, additional, iter) private(i) num_threads(split_to)
for (i=0;i < (m-additional); i+=iter)
vsLog10(iter,c+i,a+i);
//vmsLn(iter,c+i,a+i, VML_HA);
if (additional > 0)
vsLog10(additional, c+m-additional, a+m-additional);
//vmsLn(additional, c+m-additional, a+m-additional, VML_HA);
}
在Python中：
from ctypes import CDLL, c_int, c_void_p
def log_omp(xs, out):
lib = CDLL('./log_omp.so')
lib.log_omp.argtypes = [c_int, np.ctypeslib.ndpointer(dtype=np.float32), np.ctypeslib.ndpointer(dtype=np.float32)]
lib.log_omp.restype = c_void_p
n = xs.shape[0]
out = np.empty(n, np.float32)
lib.log_omp(n, out, xs)
return out
Cython 代码（在 ipython 笔记本中，因此有 %% 魔法）：
%%cython --compile-args=-fopenmp --link-args=-fopenmp
import numpy as np
cimport numpy as np
from libc.math cimport log
from cython.parallel cimport prange
import cython
@cython.boundscheck(False)
def cylog(np.ndarray[np.float32_t, ndim=1] a not None,
np.ndarray[np.float32_t, ndim=1] out=None):
if out is None:
out = np.empty((a.shape[0]), dtype=a.dtype)
cdef Py_ssize_t i
with nogil:
for i in prange(a.shape[0]):
out[i] = log(a[i])
return out
Timings:
numexpr.detect_number_of_cores() // 2
28
%env OMP_NUM_THREADS=28
x = np.abs(np.random.randn(50000000).astype('float32'))
y = x.copy()
# GCC
%timeit log_omp(x, y)
10 loops, best of 3: 21.6 ms per loop
# ICC
%timeit log_omp(x, y)
100 loops, best of 3: 9.6 ms per loop
%timeit log_VML(x, y)
100 loops, best of 3: 10 ms per loop
%timeit cylog(x, out=y)
10 loops, best of 3: 21.7 ms per loop
numexpr.set_num_threads(28)
%timeit out = numexpr.evaluate('log(x)')
100 loops, best of 3: 13 ms per loop
所以 numexpr 似乎比（糟糕的）编译的 gcc 代码做得更好，但 icc 获胜。
我发现一些有用且可耻地使用代码的资源：
http://people.duke.edu/~ccc14/sta-663/Optimization_Bakeoff.html http://people.duke.edu/~ccc14/sta-663/Optimization_Bakeoff.html
https://gist.github.com/zed/2051661 https://gist.github.com/zed/2051661
