标签: vectorization

我具有存储在numpy的阵列的3D点的列表A具有形状(N,3)和旋转矩阵R与形状(3,3).我想计算的点积R.x每个点x在A原地.天真我能做到这一点:

for n in xrange(N):
    A[n,:] = dot(R, A[n,:]) 

有没有办法用本地numpy调用来对此进行矢量化？如果重要的话,N就是几千的订单.

我正在做一些图像处理,我从矢量化中受益.我有一个矢量化ok的函数,但是我无法说服编译器输入和输出缓冲区没有重叠,因此不需要进行别名检查.我应该可以使用__restrict__,但如果缓冲区没有定义为__restrict__作为函数参数到达时,没有办法说服编译器我绝对确定2个缓冲区永远不会重叠.

这是功能:

__attribute__((optimize("tree-vectorize","tree-vectorizer-verbose=6")))
void threshold(const cv::Mat& inputRoi, cv::Mat& outputRoi, const unsigned char th) {

    const int height = inputRoi.rows;
    const int width = inputRoi.cols;

    for (int j = 0; j < height; j++) {
        const uint8_t* __restrict in = (const uint8_t* __restrict) inputRoi.ptr(j);
        uint8_t* __restrict out = (uint8_t* __restrict) outputRoi.ptr(j);
        for (int i = 0; i < width; i++) {
           out[i] = (in[i] < valueTh) ? 255 : 0;
        }
    }
}

我可以说服编译器不执行别名检查的唯一方法是将内部循环放在一个单独的函数中,其中指针被定义为__restrict__参数.如果我将此内部函数声明为内联,则再次激活别名检查.

您可以通过此示例查看效果,我认为这是一致的: …

也许这个问题非常愚蠢.

我试图"矢量化"以下循环:

set.seed(0)
x <- round(runif(10), 2)
# [1] 0.90 0.27 0.37 0.57 0.91 0.20 0.90 0.94 0.66 0.63
sig <- sample.int(10)
# [1]  1  2  9  5  3  4  8  6  7 10
for (i in seq_along(sig)) x[i] <- x[sig[i]]
x
# [1] 0.90 0.27 0.66 0.91 0.66 0.91 0.94 0.91 0.94 0.63

我认为这很简单,x[sig]但结果并不匹配.

set.seed(0)
x <- round(runif(10), 2)
x[] <- x[sig]
x
# [1] 0.90 0.27 0.66 0.91 0.37 0.57 0.94 0.20 0.90 0.63

怎么了？

备注 …

在Mathematica中有一个命令Clip[x, {min, max}] ,它给出了xfor min<=x<=max,minfor x<min和maxfor x>max,以及

http://reference.wolfram.com/mathematica/ref/Clip.html (镜像)

在R中实现这一目标的最快方法是什么？理想情况下,它应该是一个可列表的函数,理想情况下应该可以处理单个值,向量,矩阵或数据帧......

汤姆,欢呼声

正在寻找一个如何避免在我的Matlab代码中使用循环的想法,我在SE上的一个问题下找到了以下评论:

自从Matlab ... euhm,R2008a以来,声明"for循环在Matlab中很慢"的说法已不再正确.

和

您是否尝试过针对for循环与已有的循环进行基准测试？有时它比矢量化代码更快......

所以我想问一下,有没有常用的方法在Matlab中测试一个进程的速度？用户是否可以在某个地方看到流程需要多长时间,或者唯一的方法是将流程延长几分钟以便比较彼此之间的时间？

我有一个数组的值,它应该是单调的(比如减少),但是随机区域的值随着索引而增加.

我需要一个数组,其中每个区域都替换为它前面的值,以便对结果数组进行排序.

所以如果给定的数组是:

a = np.array([10.0, 9.5, 8.0, 7.2, 7.8, 8.0, 7.0, 5.0, 3.0, 2.5, 3.0, 2.0])

我想要结果

b = np.array([10.0, 9.5, 8.0, 7.2, 7.2, 7.2, 7.0, 5.0, 3.0, 2.5, 2.5, 2.0])

这是一个图形表示:

我知道如何用Python循环实现它,但有没有办法用NumPy机器来做到这一点？

Python代码为清晰起见:

b = np.array(a)
for i in range(1, b.size):
    if b[i] > b[i-1]:
        b[i] = b[i-1]

我在SSE 4.2和AVX 2的2个向量之间矢量化了点积,如下所示.该代码使用GCC 4.8.4和-O2优化标志进行编译.正如预期的那样,两者的性能都有所提高(和AVX 2比SSE 4.2快),但是当我用PAPI分析代码时,我发现未命中的总数(主要是L1和L2)增加了很多:

没有矢量化:

PAPI_L1_TCM: 784,112,091
PAPI_L2_TCM: 195,315,365
PAPI_L3_TCM: 79,362

使用SSE 4.2:

PAPI_L1_TCM: 1,024,234,171
PAPI_L2_TCM: 311,541,918
PAPI_L3_TCM: 68,842

使用AVX 2:

PAPI_L1_TCM: 2,719,959,741
PAPI_L2_TCM: 1,459,375,105
PAPI_L3_TCM: 108,140

我的代码可能有问题或者这种行为是否正常？

AVX 2代码:

double vec_dotProduct(const vec& vecs, const unsigned int& start_a, const unsigned int& start_b, const int& n) {
    double dot = 0;
    register int i = 0;
    const int loopBound = n-3;

    __m256d vsum, vecPi, vecCi, vecQCi;

    vsum = _mm256_set1_pd(0);

    double * const pA = vecs.x+start_a ;
    double * const …

我想矢量化两个内存对齐数组的乘法.我没有找到任何方法在AVX/AVX2中乘以64*64位,所以我只是循环展开和AVX2加载/存储.有更快的方法吗？

注意:我不想保存每次乘法的高半结果.

void multiply_vex(long *Gi_vec, long q, long *Gj_vec){

    int i;
    __m256i data_j, data_i;
    __uint64_t *ptr_J = (__uint64_t*)&data_j;
    __uint64_t *ptr_I = (__uint64_t*)&data_i;


    for (i=0; i<BASE_VEX_STOP; i+=4) {
        data_i = _mm256_load_si256((__m256i*)&Gi_vec[i]);
        data_j = _mm256_load_si256((__m256i*)&Gj_vec[i]);

        ptr_I[0] -= ptr_J[0] * q;
        ptr_I[1] -= ptr_J[1] * q;
        ptr_I[2] -= ptr_J[2] * q;
        ptr_I[3] -= ptr_J[3] * q;

        _mm256_store_si256((__m256i*)&Gi_vec[i], data_i);
    }


    for (; i<BASE_DIMENSION; i++)
        Gi_vec[i] -= Gj_vec[i] * q;
}

更新: 我正在使用Haswell微体系结构和ICC/GCC编译器.所以AVX和AVX2都很好.我在乘法循环展开后-=用C inrisic 替换_mm256_sub_epi64它,在那里得到一些加速.目前,它是ptr_J[0] *= q; ...

我用, …

我有一个c包含许多不同列的数据框.此外,arr是一个对应于以下子集的数据帧c:arr = c[c['A_D'] == 'A'].

我的代码的主要思想是迭代c-dataframe中的所有行,并搜索arr应该发生某些特定条件的所有可能情况(在数据框中):

• 只需要迭代行c['A_D'] == D和c['Already_linked'] == 0
• 将hour在arr数据帧必须小于hour_aux在c数据帧
• 列Already_linked的的arr数据框必须是零:arr.Already_linked == 0
• c和数据框中的Terminal和Operator需要相同arr

现在,使用布尔索引和groupby get_group存储条件:

• 分组arr数据帧以选择相同的运营商和终端:)g = groups.get_group((row.Operator, row.Terminal)
• 仅选择小时数小于c数据框中小时的到达时间以及Already_linked == 0的位置:vb = g[(g.Already_linked==0) & (g.hour<row.hour_aux)]

对于c验证所有条件的数据框中的每一行,vb都会创建一个数据框.当然,这个数据帧在每次迭代中都有不同的长度.创建vb数据框后,我的目标是选择vb数据框的索引,以最小化vb.STARTc …

我有两个大数组,大约1000行和1000列.我需要比较这些数组的每个元素,如果相应的元素相等,则将1存储在另一个数组中.

我可以使用for循环但这需要很长时间.我怎么能更快地做到这一点？