将数字乘以矩阵

Question

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

我是说。

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

有没有使用浮点数尽可能快地进行这种乘法的算法？

非常感谢您的帮助。

编辑

我使用了 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 (...) tiemp2 - tiempo1 变量返回值 0.096924，但如果我使用自己的函数 (matriz_matriz(...)) tiempo2 - tiempo1 返回值 0.046271...发生什么了？我的函数比 Cblas 快...

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

非常感谢。

Answer 1

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

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

Answer 2

2 x 2 矩阵的版本（基于link）：

#include<iostream>
using namespace std;

int main()
{
    const int rows = 2;
    const int cols = 2;

    float a[4]={1,2,3,4};
    float b[4]={1,2,3,4};
    float c[4]={0,0,0,0};

    for (int i = 0; i <rows; i++) {
        for (int j = 0; j <cols; j++) 
        {   
            float sum = 0.0;
            for (int k = 0; k < rows; k++)
                sum = sum + a[i * cols + k] * b[k * cols + j]; 
            c[i * cols + j] = sum;
        }   
    }   
    for (int ix =0; ix <4; ++ix)
            cout << c[ix] << ' ';

}