标签: vectorization

我有一堆时间序列,每个时间序列由两个组件描述,一个时间戳向量(以秒为单位),以及一个测量值向量.时间向量是不均匀的(即以非常规间隔采样)

我试图计算每个1分钟间隔值的平均值/ SD(采用X分钟间隔,计算其平均值,采用下一个间隔,......).

我当前的实现使用循环.这是我到目前为止的样本:

t = (100:999)' + rand(900,1);       %' non-uniform time
x = 5*rand(900,1) + 10;             % x(i) is the value at time t(i)

interval = 1;         % 1-min interval
tt = ( floor(t(1)):interval*60:ceil(t(end)) )';  %' stopping points of each interval
N = length(tt)-1;

mu = zeros(N,1);
sd = zeros(N,1);

for i=1:N
    indices = ( tt(i) <= t & t < tt(i+1) ); % find t between tt(i) and tt(i+1)
    mu(i) = mean( x(indices) );
    sd(i) = std( x(indices) …

我试图对包含在数学库中使用'pow'函数的循环进行矢量化.我知道英特尔编译器支持使用'pow'作为sse指令 - 但我似乎无法使用gcc运行(我认为).这是我正在使用的情况:

int main(){
        int i=0;
        float a[256],
        b[256];

        float x= 2.3;


        for  (i =0 ; i<256; i++){
                a[i]=1.5;
        }

        for (i=0; i<256; i++){
                b[i]=pow(a[i],x);
        }

        for (i=0; i<256; i++){
                b[i]=a[i]*a[i];
        }
    return 0;

}

我正在编译以下内容:

gcc -O3 -Wall -ftree-vectorize -msse2 -ftree-vectorizer-verbose=5 code.c -o runthis

这是在使用gcc版本4.2的os X 10.5.8上(我也使用4.5并且无法判断它是否已经向量化了 - 因为它根本没有输出任何内容).似乎没有一个循环矢量化 - 是否存在一个对齐问题或者我需要使用限制的其他问题？如果我将其中一个循环写为函数,我会得到更详细的输出(代码):

void pow2(float *a, float * b, int n) {
        int i;
        for (i=0; i<n; i++){
                b[i]=a[i]*a[i];
        }
}

输出(使用7级详细输出):

note: not vectorized: can't determine dependence between *D.2878_13 …

在用基于范围的循环替换我的许多"旧"for循环之前,我使用visual studio 2013进行了一些测试:

std::vector<int> numbers;

for (int i = 0; i < 50; ++i) numbers.push_back(i);

int sum = 0;

//vectorization
for (auto number = numbers.begin(); number != numbers.end(); ++number) sum += *number;

//vectorization
for (auto number = numbers.begin(); number != numbers.end(); ++number) {
    auto && ref = *number;
    sum += ref;
}

//definition of range based for loops from http://en.cppreference.com/w/cpp/language/range-for
//vectorization
for (auto __begin = numbers.begin(),
    __end = numbers.end();
    __begin != __end; ++__begin) {
    auto && ref = *__begin;
    sum …

根据英特尔内部指南,

• vxorpd ymm, ymm, ymm:计算a和b中打包的双精度(64位)浮点元素的按位XOR,并将结果存储在dst中.
• vpxor ymm, ymm, ymm:计算a和b中256位(表示整数数据)的按位XOR,并将结果存储在dst中.

两者有什么区别？在我看来,两个指令都会对ymm寄存器的所有256位执行按位异或.如果我使用vxorpd整数数据会有任何性能损失(反之亦然)？

我有两个函数用于在Julia中以数字方式确定pi.第二个函数(我认为是矢量化的)比第一个函数慢.为什么矢量化速度较慢？是否有规则何时进行矢量化以及何时不进行？

function determine_pi(n)
    area = zeros(Float64, n);
    sum = 0;
    for i=1:n
        if ((rand()^2+rand()^2) <=1)
            sum = sum + 1;
        end
            area[i] = sum*1.0/i;
    end
    return area
end

和另一个功能

function determine_pi_vec(n)
    res = cumsum(map(x -> x<=1?1:0, rand(n).^2+rand(n).^2))./[1:n]
    return res
end

当运行n = 10 ^ 7时,以下是执行时间(运行几次后)

n=10^7
@time returnArray = determine_pi(n)
#output elapsed time: 0.183211324 seconds (80000128 bytes allocated)
@time returnArray2 = determine_pi_vec(n);
#elapsed time: 2.436501454 seconds (880001336 bytes allocated, 30.71% gc time)

一个可重现的代码示例我正在尝试矢量化.

cutOffs <- seq(1,10,0.2)

plotOutput <- matrix(nrow=length(cutOffs), ncol=2)
colnames(plotOutput) <- c("x","y")
plotOutput[,"y"] <- cutOffs

for(plotPoint in 1:length(cutOffs))
{
  plotOutput[plotPoint, "x"] <-
    nrow(iris[ which(iris$Sepal.Length > cutOffs[plotPoint] &
                   iris$Sepal.Width > cutOffs[plotPoint]), ])
}

plotOutput

特别是我想要找到的是,如果有一种方法来矢量化这部分.

nrow(iris[ which(iris$Sepal.Length > cutOffs[plotPoint] &
                   iris$Sepal.Width > cutOffs[plotPoint]), ])

假设我是使用plyr库或某种形式的应用程序,可能没有太多加速,这正是我正在寻找的.从根本上说,我试图看看是否存在一些我在搜索时忽略或设法错过的矢量化技术.

更新:

Unit: milliseconds
  expr         min          lq        mean      median          uq         max neval
  op() 33663.39700 33663.39700 33663.39700 33663.39700 33663.39700 33663.39700     1
  jr()  3976.53088  3976.53088  3976.53088  3976.53088  3976.53088  3976.53088     1
  dd()  4253.21050  4253.21050  4253.21050  4253.21050  4253.21050  4253.21050     1
 exp()  5085.45331 …

我需要检查所有向量元素是否为非零.到目前为止,我找到了以下解 有一个更好的方法吗？我在Linux/x86_64上使用gcc 4.8.2,指令直到SSE4.2.

typedef char ChrVect __attribute__((vector_size(16), aligned(16)));

inline bool testNonzero(ChrVect vect)
{
    const ChrVect vzero = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
    return (0 == (__int128_t)(vzero == vect));
}

更新:上面的代码被编译为以下汇编代码(当编译为非内联函数时):

movdqa  %xmm0, -24(%rsp)
pxor    %xmm0, %xmm0
pcmpeqb -24(%rsp), %xmm0
movdqa  %xmm0, -24(%rsp)
movq    -24(%rsp), %rax
orq -16(%rsp), %rax
sete    %al
ret

有一个相对简单的代码块,它循环遍历两个数组,相乘并累加:

import numpy as np

a = np.array([1, 2, 4, 6, 7, 8, 9, 11])
b = np.array([0.01, 0.2, 0.03, 0.1, 0.1, 0.6, 0.5, 0.9])
c = []
d = 0
for i, val in enumerate(a):
    d += val
    c.append(d)
    d *= b[i]

有没有办法在没有迭代的情况下做到这一点？我想可以使用cumsum/cumprod,但我无法弄清楚如何.当你逐步分解正在发生的事情时,它看起来像这样:

# 0: 0 + a[0]
# 1: ((0 + a[0]) * b[0]) + a[1]
# 2: ((((0 + a[0]) * b[0]) + a[1]) * b[1]) + a[2]

编辑澄清:我对列表(或数组)感兴趣c.

我想将大小的矩阵分成大小的k x l块,n x n考虑到一个问题o(就像Mathematica的分区函数那样).

例如,给定矩阵之A类的

A <- matrix(seq(1:16), nrow = 4, ncol = 4)

     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

和块大小= 3,偏移= 1,我想要输出我得到的四个子矩阵

A[1:3, 1:3]
A[1:3, 2:4]
A[2:4, 1:3]
A[2:4, 2:4]

如果offset等于2或3,则此示例的输出应该只是我得到的子矩阵

A[1:3, 1:3]

我该如何对此进行矢量化？

我正在使用 Matlab。

我有一个二元方阵。对于每一行，有一个或多个 1 的条目。我想遍历此矩阵的每一行并返回这些 1 的索引并将它们存储在单元格的条目中。

我想知道是否有一种方法可以在不循环此矩阵的所有行的情况下执行此操作，因为在 Matlab 中 for 循环非常慢。

例如，我的矩阵

M = 0 1 0
    1 0 1
    1 1 1 

然后最终，我想要类似的东西

A = [2]
    [1,3]
    [1,2,3]

A细胞也是如此。

有没有办法在不使用 for 循环的情况下实现这个目标，目的是更快地计算结果？