c++ - 将数字乘以矩阵

Question

谁能告诉我将一系列数字乘以矩阵的最佳方法是什么？

我是说。

我看过矩阵乘法的算法，但要将数字相乘为 matrix1 [4] [4] 和 matrix2 [4] [4]。但是，我想将数字乘以 matrix1 [16] 和 matrix2 [16]。

是否有使用 float 尽可能快地进行此乘法的算法？

非常感谢您的帮助。

编辑

我使用了 cBLAS 并进行了一些速度测试，我对结果感到惊讶。

#include <stdio.h>
#include <stdlib.h>
#include <cblas.h>
#include  <GL/glfw.h>

    void matriz_matriz(float *matriz1,float *matriz2,float *matrizr){
      matrizr[0]  = (matriz1[0]*matriz2[0])+(matriz1[4]*matriz2[1])  +(matriz1[8]*matriz2[2])  +(matriz1[12]*matriz2[3]);
      matrizr[1]  = (matriz1[1]*matriz2[0])+(matriz1[5]*matriz2[1])  +(matriz1[9]*matriz2[2])  +(matriz1[13]*matriz2[3]);
      matrizr[2]  = (matriz1[2]*matriz2[0])+(matriz1[6]*matriz2[1])  +(matriz1[10]*matriz2[2]) +(matriz1[14]*matriz2[3]);
      matrizr[3]  = (matriz1[3]*matriz2[0])+(matriz1[7]*matriz2[1])  +(matriz1[11]*matriz2[2]) +(matriz1[15]*matriz2[3]);

      matrizr[4]  = (matriz1[0]*matriz2[4])+(matriz1[4]*matriz2[5])  +(matriz1[8]*matriz2[6])  +(matriz1[12]*matriz2[7]);
      matrizr[5]  = (matriz1[1]*matriz2[4])+(matriz1[5]*matriz2[5])  +(matriz1[9]*matriz2[6])  +(matriz1[13]*matriz2[7]);
      matrizr[6]  = (matriz1[2]*matriz2[4])+(matriz1[6]*matriz2[5])  +(matriz1[10]*matriz2[6]) +(matriz1[14]*matriz2[7]);
      matrizr[7]  = (matriz1[3]*matriz2[4])+(matriz1[7]*matriz2[5])  +(matriz1[11]*matriz2[6]) +(matriz1[15]*matriz2[7]);

      matrizr[8]  = (matriz1[0]*matriz2[8])+(matriz1[4]*matriz2[9])  +(matriz1[8]*matriz2[10]) +(matriz1[12]*matriz2[11]);
      matrizr[9]  = (matriz1[1]*matriz2[8])+(matriz1[5]*matriz2[9])  +(matriz1[9]*matriz2[10]) +(matriz1[13]*matriz2[11]);
      matrizr[10] = (matriz1[2]*matriz2[8])+(matriz1[6]*matriz2[9])  +(matriz1[10]*matriz2[10])+(matriz1[14]*matriz2[11]);
      matrizr[11] = (matriz1[3]*matriz2[8])+(matriz1[7]*matriz2[9])  +(matriz1[11]*matriz2[10])+(matriz1[15]*matriz2[11]);

      matrizr[12] = (matriz1[0]*matriz2[12])+(matriz1[4]*matriz2[13])+(matriz1[8]*matriz2[14]) +(matriz1[12]*matriz2[15]);
      matrizr[13] = (matriz1[1]*matriz2[12])+(matriz1[5]*matriz2[13])+(matriz1[9]*matriz2[14]) +(matriz1[13]*matriz2[15]);
      matrizr[14] = (matriz1[2]*matriz2[12])+(matriz1[6]*matriz2[13])+(matriz1[10]*matriz2[14])+(matriz1[14]*matriz2[15]);
      matrizr[15] = (matriz1[3]*matriz2[12])+(matriz1[7]*matriz2[13])+(matriz1[11]*matriz2[14])+(matriz1[15]*matriz2[15]);
    }


    int main(){
      int i;
      double tiempo1;
      double tiempo2;

      glfwInit();

      float *mat0 = NULL;
      float *mat1 = NULL;
      float *mat2 = NULL;

      mat0  = (float *)malloc(16 * sizeof(float));
      mat1  = (float *)malloc(16 * sizeof(float));
      mat2  = (float *)malloc(16 * sizeof(float));

      mat0[0]  =  1.0;
      mat0[1]  =  0.0;
      mat0[2]  =  0.0;
      mat0[3]  =  0.0;
      mat0[4]  =  0.0;
      mat0[5]  =  1.0;
      mat0[6]  =  0.0;
      mat0[7]  =  0.0;
      mat0[8]  =  0.0;
      mat0[9]  =  0.0;
      mat0[10] =  1.0;
      mat0[11] =  0.0;
      mat0[12] =  3.281897;
      mat0[13] =  4.714289;
      mat0[14] =  5.124306;
      mat0[15] =  1.0;

      mat1[0]  =  1.0;
      mat1[1]  =  0.0;
      mat1[2]  =  0.0;
      mat1[3]  =  0.0;
      mat1[4]  =  0.0;
      mat1[5]  =  0.924752;
      mat1[6]  =  0.380570;
      mat1[7]  =  0.0;
      mat1[8]  =  0.0;
      mat1[9]  = -0.380570;
      mat1[10] =  0.924752;
      mat1[11] =  0.0;
      mat1[12] =  0.0;
      mat1[13] =  0.0;
      mat1[14] =  0.0;
      mat1[15] =  1.0;

      mat2[0]  =  1.0;
      mat2[1]  =  0.0;
      mat2[2]  =  0.0;
      mat2[3]  =  0.0;
      mat2[4]  =  0.0;
      mat2[5]  =  1.0;
      mat2[6]  =  0.0;
      mat2[7]  =  0.0;
      mat2[8]  =  0.0;
      mat2[9]  =  0.0;
      mat2[10] =  1.0;
      mat2[11] =  0.0;
      mat2[12] =  0.0;
      mat2[13] =  0.0;
      mat2[14] =  0.0;
      mat2[15] =  1.0;

       tiempo1 = glfwGetTime();

       for(i=0;i<100000;i++){
        matriz_matriz(mat0,mat1,mat2);
        //cblas_sgemm(CblasRowMajor,CblasNoTrans,CblasNoTrans,4,4,4,1.0f,mat0,4,mat1,4,0.0f,mat2,4);
       }

      tiempo2 = glfwGetTime();
      printf("Tiempo total: %f\n",tiempo2-tiempo1);

      for(i=0;i<16;i++)printf("valor[%i]: %f\n",i,mat2[i]);

      free(mat0);
      free(mat1);
      free(mat2);

      system("pause");

      glfwTerminate();
      return 0;
    }

如果我使用函数 cblas_sgemm (...) tiempo2 - tiempo1 变量返回值 0.096924，但如果我使用我自己的函数 (matriz_matriz(...)) tiempo2 - tiempo1 返回值 0.046271...会发生什么？我的功能比 Cblas 更快......

此测试是在配备 Pentium 3 处理器的 PC 上进行的。谁能告诉我发生了什么？

非常感谢。

Answer 1

老实说，如果你正在做任何一种线性代数，到目前为止，你最好的选择是使用为此目的而设计的库，例如 BLAS、LAPACK、等等。您将很难用自己的代码接近他们的速度。

矩阵-矩阵运算是 BLAS 级别 3，您想要的特定运算是 SGEMM() for float 和 DGEMM()对于 double。英特尔硬件上最快的 BLAS 实现是 OpenBLAS(源自 GotoBLAS)和英特尔 MKL(数学内核库)中的 BLAS 实现。如果您自己编译，ATLAS 也非常快。