代码之家  ›  专栏  ›  技术社区  ›  malang

cblas-gemm时间依赖于输入矩阵值-Ubuntu 14.04

  •  2
  • malang  · 技术社区  · 8 年前

    这是一个 扩大 属于我的 earlier question ,但我是分开问的,因为我真的很沮丧,所以请不要投反对票!

    问题:cblas_sgemm调用使用 大量的零 与相同的cblassgemm调用稠密矩阵相比?

    我知道gemv是为矩阵向量乘法而设计的,但如果它花费的时间更少,为什么我不能使用gemm进行向量矩阵乘法,特别是对于稀疏矩阵

    下面给出了一个简短的代表代码。它要求输入一个值,然后用该值填充向量。然后用索引替换每32个值。因此,如果我们输入“0”,那么我们会得到一个稀疏向量,但对于其他值,我们会得到稠密向量。

    #include <iostream>
    #include <stdio.h>
    #include <time.h>
    #include <cblas.h>
    
    using namespace std;
    
    int main()
    {
    const int m = 5000; 
    
    timespec blas_start, blas_end;
    long totalnsec; //total nano sec
    double totalsec, totaltime;
    int i, j;
    float *A = new float[m]; // 1 x m
    float *B = new float[m*m]; // m x m
    float *C = new float[m]; // 1 x m
    
    float input;
    cout << "Enter a value to populate the vector (0 for sparse) ";
    cin >> input; // enter 0 for sparse
    
    // input martix A: every 32nd element is non-zero, rest of the values = input
    for(i = 0; i < m; i++)
    {       
    A[i] = input;
    if( i % 32 == 0)    //adjust for sparsity
            A[i] = i;
    }
    
    // input matrix B: identity matrix
    for(i = 0; i < m; i++)
            for(j = 0; j < m; j++)
                B[i*m + j] = (i==j);
    
    clock_gettime(CLOCK_REALTIME, &blas_start);
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, m, m, 1.0f, A, m, B, m, 0.0f, C, m);
    clock_gettime(CLOCK_REALTIME, &blas_end);
    
    /* for(i = 0; i < m; i++)
            printf("%f ", C[i]);
    printf("\n\n");    */
    
    // Print time
    totalsec = (double)blas_end.tv_sec - (double)blas_start.tv_sec;
    totalnsec = blas_end.tv_nsec - blas_start.tv_nsec;
    if(totalnsec < 0)
    {
        totalnsec += 1e9;
        totalsec -= 1;
    }
    totaltime = totalsec + (double)totalnsec*1e-9;
    cout<<"Duration = "<< totaltime << "\n";
    
    return 0;
    }
    

    我用blas 3.0在Ubuntu 14.04中运行如下

    erisp@ubuntu:~/uas/stackoverflow$ g++ gemmcomp.cpp -o gemmcomp.o -lblas
    erisp@ubuntu:~/uas/stackoverflow$ ./gemmcomp.o
    Enter a value to populate the vector (0 for sparse) 5
    Duration = 0.0291558
    erisp@ubuntu:~/uas/stackoverflow$ ./gemmcomp.o
    Enter a value to populate the vector (0 for sparse) 0
    Duration = 0.000959521
    

    编辑

    如果我用下面的gemv调用替换gemm调用,那么矩阵稀疏性无关紧要

    //cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 1, m, m, 1.0f, A, m, B, m, 0.0f, C, m);
    cblas_sgemv(CblasRowMajor, CblasNoTrans, m, m, 1.0f, B, m, A, 1, 0.0f, C, 1);
    

    后果

    erisp@ubuntu:~/uas/stackoverflow$ g++ gemmcomp.cpp -o gemmcomp.o -lblas
    erisp@ubuntu:~/uas/stackoverflow$ ./gemmcomp.o
    Enter a value to populate the vector (0 for sparse) 5
    Duration = 0.0301581
    erisp@ubuntu:~/uas/stackoverflow$ ./gemmcomp.o
    Enter a value to populate the vector (0 for sparse) 0
    Duration = 0.0299282
    

    但问题是,我正在尝试使用cublas优化其他人的代码,他成功地 高效地 使用gemm执行此矢量矩阵乘法。因此,我必须对此进行测试,或者明确证明这一呼吁是 不准确的

    编辑

    今天,我甚至更新了我的blas库,使用

    sudo apt-get install libblas-dev liblapack-dev
    

    EDIT:按照不同参与者的建议执行以下命令

    erisp@ubuntu:~/uas/stackoverflow$ ll -d /usr/lib/libblas* /etc/alternatives/libblas.*
    lrwxrwxrwx 1 root root   26 مارچ  13  2015 /etc/alternatives/libblas.a -> /usr/lib/libblas/libblas.a
    lrwxrwxrwx 1 root root   27 مارچ  13  2015 /etc/alternatives/libblas.so -> /usr/lib/libblas/libblas.so
    lrwxrwxrwx 1 root root   29 مارچ  13  2015 /etc/alternatives/libblas.so.3 -> /usr/lib/libblas/libblas.so.3
    lrwxrwxrwx 1 root root   29 مارچ  13  2015 /etc/alternatives/libblas.so.3gf -> /usr/lib/libblas/libblas.so.3
    drwxr-xr-x 2 root root 4096 مارچ  13  2015 /usr/lib/libblas/
    lrwxrwxrwx 1 root root   27 مارچ  13  2015 /usr/lib/libblas.a -> /etc/alternatives/libblas.a
    lrwxrwxrwx 1 root root   28 مارچ  13  2015 /usr/lib/libblas.so -> /etc/alternatives/libblas.so
    lrwxrwxrwx 1 root root   30 مارچ  13  2015 /usr/lib/libblas.so.3 -> /etc/alternatives/libblas.so.3
    lrwxrwxrwx 1 root root   32 مارچ  13  2015 /usr/lib/libblas.so.3gf -> /etc/alternatives/libblas.so.3gf
    
    erisp@ubuntu:~/uas/stackoverflow$ ldd ./gemmcomp.o
        linux-gate.so.1 =>  (0xb76f6000)
        libblas.so.3 => /usr/lib/libblas.so.3 (0xb765e000)
        libstdc++.so.6 => /usr/lib/i386-linux-gnu/libstdc++.so.6 (0xb7576000)
         libc.so.6 => /lib/i386-linux-gnu/libc.so.6 (0xb73c7000)
         libm.so.6 => /lib/i386-linux-gnu/libm.so.6 (0xb7381000)
    /lib/ld-linux.so.2 (0xb76f7000)
         libgcc_s.so.1 => /lib/i386-linux-gnu/libgcc_s.so.1 (0xb7364000)
    
    1 回复  |  直到 7 年前
        1
  •  6
  •   Robert Crovella    8 年前

    问题:cblas_sgemm调用对包含大量零的矩阵的时间比对密集矩阵的相同cblas_sgemm调用要少得多,其原因可能是什么?

    似乎BLAS实现由 the default libblas-dev package 对于Ubuntu 14.04(可能还有其他Ubuntu发行版),包含了对某些矩阵元素为零的情况的优化。

    对于Ubuntu 14.04,BLAS(和cblas)实现/包的源代码可以从 here .

    解压缩该存档后,我们有一个 cblas/src 包含 cblas API,我们还有另一个 src 目录,其中包含各种blas例程的F77实现。

    对于cblas_sgemm,当参数 CblasRowMajor 已指定,则 cblas/src/cblas_sgemm.c 代码将按如下方式调用底层fortran例程:

    void cblas_sgemm(const enum CBLAS_ORDER Order, const enum CBLAS_TRANSPOSE TransA,
                 const enum CBLAS_TRANSPOSE TransB, const int M, const int N,
                 const int K, const float alpha, const float  *A,
                 const int lda, const float  *B, const int ldb,
                 const float beta, float  *C, const int ldc)
    {
    ...
    } else if (Order == CblasRowMajor)
    ...
      F77_sgemm(F77_TA, F77_TB, &F77_N, &F77_M, &F77_K, &alpha, B, &F77_ldb, A, &F77_lda, &beta, C, &F77_ldc);
    

    请注意,对于此行主调用 A B 当传递给 F77_sgemm 常规这是明智的,但我不会在这里深入探讨原因。值得注意的是 A. 已成为 B fortran调用/代码中,以及 B 已成为 A. .

    当我们在中检查相应的fortran例程时 src/sgemm.f ,我们看到以下代码序列:

    *
    *     Start the operations.
    *
          IF (NOTB) THEN
              IF (NOTA) THEN
    *
    *           Form  C := alpha*A*B + beta*C.
    *
                  DO 90 J = 1,N
                      IF (BETA.EQ.ZERO) THEN
                          DO 50 I = 1,M
                              C(I,J) = ZERO
       50                 CONTINUE
                      ELSE IF (BETA.NE.ONE) THEN
                          DO 60 I = 1,M
                              C(I,J) = BETA*C(I,J)
       60                 CONTINUE
                      END IF
                      DO 80 L = 1,K
                          IF (B(L,J).NE.ZERO) THEN        ***OPTIMIZATION
                              TEMP = ALPHA*B(L,J)
                              DO 70 I = 1,M
                                  C(I,J) = C(I,J) + TEMP*A(I,L)
       70                     CONTINUE
                          END IF
       80             CONTINUE
       90         CONTINUE
    

    上面是处理以下情况的代码段: A. 且无转置 B (对于这个cblas行主要测试用例来说是这样的)。矩阵行/列乘法运算在我添加注释的循环开始处处理 ***OPTIMIZATION 特别是,如果矩阵元素 B(L,J) 为零,则跳过第70行关闭的DO循环。但请记住 B 这里对应于 A. 矩阵传递给 cblas_sgemm 常规

    跳过这个do循环可以使以这种方式实现的sgemm函数在 A. 矩阵传递给 cblas_sgemm公司 当指定了行major时。

    从实验上看,并不是所有的blas实现都有这种优化。在完全相同的平台上测试,但使用 libopenblas-dev 而不是 libblas-dev 不提供这种加速,即当 A. 矩阵大多为零,而不是零。

    请注意,此处的fortran(77)代码似乎与旧版本的 sgemm.f 例行程序,例如 here 。我发现此fortran例程的较新发布版本不包含此优化,例如 here .