умножение чисел как матрицы

Question

Может ли кто-нибудь сказать мне, что является лучшим способом умножить серию чисел в виде матрицы?

Я имею в виду.

Я видел алгоритмы для матричного умножения, но должен умножать числа как матрицы1 [4] [4] и matrix2 [4] [4]. Однако я хочу умножить числа на матрицу 1 [16] и на матрицу2 [16].

Есть ли какой-либо алгоритм для этого умножения как можно быстрее с использованием чисел с плавающей точкой?

Большое спасибо за помощь.

EDIT

Я использовал 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. Может ли кто-нибудь сказать мне, что происходит?

спасибо.

Answer 1

Честно говоря, если вы делаете какой-либо линейной алгебры, безусловно, ваш лучший выбор, чтобы использовать библиотеки, предназначенные для этой цели, такие как BLAS, LaPack и т.д. Вы будете иметь очень трудное время приближаясь к своей скорости с помощью собственного кода.

Matrix-матричные операции BLAS Level 3, и особенно один вы хотите SGEMM() для float с и DGEMM() для double с. Самая быстрая реализация BLAS на оборудовании Intel - OpenBLAS (производная от GotoBLAS) и реализация BLAS в Intel MKL (библиотека математического ядра). ATLAS также очень быстро, если вы скомпилируете его самостоятельно.

Answer 2

Версия для 2 х 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] << ' '; 

}