Dlaczego mnożenie C++ z dynamiczną tablicą działa lepiej niż std :: vector version

Question

Wdrażam mnożenie C++ dla macierzy z różnymi strukturami danych i technikami (wektorami, tablicami i OpenMP) i znalazłem dziwną sytuację ... Moja dynamika wersja tablica działa lepiej:

razy:

OpenMP mult_1: godzina: 5.882000 s

tablica mult_2: czas: 1,478000 s

Moje kompilacji flagi są:

/usr/bin/g ++ -fopenmp -pthread -std = C++ 1r -O3

wersja C++ wektor

typedef std::vector<std::vector<float>> matrix_f; 
void mult_1 (const matrix_f & matrixOne, const matrix_f & matrixTwo, matrix_f & result) { 
    const int matrixSize = (int)result.size(); 
    #pragma omp parallel for simd 
    for (int rowResult = 0; rowResult < matrixSize; ++rowResult) { 
     for (int colResult = 0; colResult < matrixSize; ++colResult) { 
      for (int k = 0; k < matrixSize; ++k) { 
       result[rowResult][colResult] += matrixOne[rowResult][k] * matrixTwo[k][colResult]; 
      } 
     } 
    } 
} 

Dynamiczna macierz wersji

void mult_2 (float * matrixOne, float * matrixTwo, float * result, int size) { 
    for (int row = 0; row < size; ++row) { 
     for (int col = 0; col < size; ++col) { 
      for (int k = 0; k < size; ++k) { 
       (*(result+(size*row)+col)) += (*(matrixOne+(size*row)+k)) * (*(matrixTwo+(size*k)+col)); 
      } 
     } 
    } 
} 

testy:

C++ wersja wektor

utils::ChronoTimer timer; 
/* set Up simple matrix */ 
utils::matrix::matrix_f matr1 = std::vector<std::vector<float>>(size,std::vector<float>(size)); 
fillRandomMatrix(matr1); 

utils::matrix::matrix_f matr2 = std::vector<std::vector<float>>(size,std::vector<float>(size)); 
fillRandomMatrix(matr2); 

utils::matrix::matrix_f result = std::vector<std::vector<float>>(size,std::vector<float>(size));  
timer.init(); 
utils::matrix::mult_1(matr1,matr2,result); 
std::printf("openmp mult_1: time: %f ms\n",timer.now()/1000); 

Dynamiczna tablica wersja

utils::ChronoTimer timer; 

float *p_matr1 = new float[size*size]; 
float *p_matr2 = new float[size*size]; 
float *p_result = new float[size*size]; 

fillRandomMatrixArray(p_matr1,size); 
fillRandomMatrixArray(p_matr2,size); 

timer.init(); 
utils::matrix::mult_2(p_matr1,p_matr2,p_result,size); 
std::printf("array mult_2: time: %f ms\n",timer.now()/1000); 

delete [] p_matr1; 
delete [] p_matr2; 
delete [] p_result; 

Sprawdzałem kilka poprzednich postów, ale nie mogłem znaleźć żadnych związanych z moim problem link, link2, link3:

UPDATE: I refactorized testy z odpowiedziami, a VectorWorks slighty lepiej:

wektor mult: Godzina: 1.194000 s

tablica mult_2: Godzina: 1,202000 s

C++ wersja wektorowa

void mult (const std::vector<float> & matrixOne, const std::vector<float> & matrixTwo, std::vector<float> & result, int size) { 
    for (int row = 0; row < size; ++row) { 
     for (int col = 0; col < size; ++col) { 
      for (int k = 0; k <size; ++k) { 
       result[(size*row)+col] += matrixOne[(size*row)+k] * matrixTwo[(size*k)+col]; 
      } 
     } 
    } 
} 

Dynamiczna tablica wersja

void mult_2 (float * matrixOne, float * matrixTwo, float * result, int size) { 
    for (int row = 0; row < size; ++row) { 
     for (int col = 0; col < size; ++col) { 
      for (int k = 0; k < size; ++k) { 
       (*(result+(size*row)+col)) += (*(matrixOne+(size*row)+k)) * (*(matrixTwo+(size*k)+col)); 
      } 
     } 
    } 
} 

Również moja vectorized wersja pracuje lepsze (0,803 s);

Answer 1

Wektor wektorów jest analogiczny do tablicy poszarpanej - tablicy, w której każdy wpis jest wskaźnikiem, a każdy wskaźnik wskazuje na inną tablicę zmiennych.

Dla porównania, podstawowa wersja macierzy to jeden blok pamięci, w którym matematyka służy do znajdowania elementów.

Użyj pojedynczego wektora, a nie wektorów wektorów i rób ręcznie matematykę. Lub użyj wektora o rozmiarze rozmiaru std::array s. Lub napisz typ pomocnika, który pobiera (jednowymiarowy) wektor float, i daje dwuwymiarowy widok tego.

Dane w sąsiadującym buforze są bardziej przyjazne dla pamięci podręcznej i optymalizacji niż dane w grupie odłączonych buforów.

template<std::size_t Dim, class T> 
struct multi_dim_array_view_helper { 
    std::size_t const* dims; 
    T* t; 
    std::size_t stride() const { 
    return 
     multi_dim_array_view_helper<Dim-1, T>{dims+1, nullptr}.stride() 
     * *dims; 
    } 
    multi_dim_array_view_helper<Dim-1, T> operator[](std::size_t i)const{ 
    return {dims+1, t+i*stride()}; 
    } 
}; 
template<class T> 
struct multi_dim_array_view_helper<1, T> { 
    std::size_t stride() const{ return 1; } 
    T* t; 
    T& operator[](std::size_t i)const{ 
    return t[i]; 
    } 
    multi_dim_array_view_helper(std::size_t const*, T* p):t(p) {} 
}; 
template<std::size_t Dims> 
using dims_t = std::array<std::size_t, Dims-1>; 
template<std::size_t Dims, class T> 
struct multi_dim_array_view_storage 
{ 
    dims_t<Dims> storage; 
}; 
template<std::size_t Dims, class T> 
struct multi_dim_array_view: 
    multi_dim_array_view_storage<Dims, T>, 
    multi_dim_array_view_helper<Dims, T> 
{ 
    multi_dim_array_view(dims_t<Dims> d, T* t): 
    multi_dim_array_view_storage<Dims, T>{std::move(d)}, 
    multi_dim_array_view_helper<Dims, T>{ 
     this->storage.data(), t 
    } 
    {} 
}; 

teraz można to zrobić:

std::vector<float> blah = { 
    01.f, 02.f, 03.f, 
    11.f, 12.f, 13.f, 
    21.f, 22.f, 23.f, 
}; 

multi_dim_array_view<2, float> view = { {3}, blah.data() }; 
for (std::size_t i = 0; i < 3; ++i) 
{ 
    std::cout << "["; 
    for (std::size_t j = 0; j < 3; ++j) 
    std::cout << view[i][j] << ","; 
    std::cout << "]\n"; 
} 

live example

Żadne dane nie są kopiowane w klasie widoku. Po prostu zapewnia widok płaskiej tablicy, która jest tablicą wielowymiarową.

Answer 2

Twoje podejście jest zupełnie inna:

• W „tablicy dynamicznej” wersji można przydzielić jeden kawałek pamięci dla każdej matrycy i zmapować wiersze macierzy na tej jednowymiarowej zakresie pamięci.

• W wersji "wektorowej" używane są wektory wektorów "rzeczywistych" i "dynamicznie" dwuwymiarowych, co oznacza, że ​​pozycja przechowywania każdego wiersza macierzy nie ma związku z innymi wierszami.

Co prawdopodobnie chcesz zrobić jest:

• Zastosowanie vector<float>(size*size) i wykonywać tę samą mapowanie robisz w „tablicy dynamicznej” przykład ręcznie lub

• Write klasa, która wewnętrznie obsługuje mapowanie dla ciebie i zapewnia dwuwymiarowy interfejs dostępu (T& operator()(size_t, size_t) lub coś w rodzaju row_proxy operator[](size_t), gdzie row_proxy z kolei ma)

Answer 3

Ma to na celu jedynie egzekwowanie teorii (w praktyce) dotyczącej ciągłej pamięci.

Po jakiejś analizę kodu wygenerowanego z g ++ (-O2) źródło można znaleźć na stronie: https://gist.github.com/42be237af8e3e2b1ca03

Odpowiedni kod generowany dla wersji tablicy brzmi:

.L3: 
    lea r9, [r13+0+rbx]    ; <-------- KEEPS THE ADDRESS 
    lea r11, [r12+rbx] 
    xor edx, edx 
.L7: 
    lea r8, [rsi+rdx] 
    movss xmm1, DWORD PTR [r9] 
    xor eax, eax 
.L6: 
    movss xmm0, DWORD PTR [r11+rax*4] 
    add rax, 1 
    mulss xmm0, DWORD PTR [r8] 
    add r8, r10 
    cmp ecx, eax 
    addss xmm1, xmm0 
    movss DWORD PTR [r9], xmm1  ; <------------ ADDRESS IS USED 
    jg .L6 
    add rdx, 4 
    add r9, 4      ; <--- ADDRESS INCREMENTED WITH SIZE OF FLOAT 
    cmp rdx, rdi 
    jne .L7 
    add ebp, 1 
    add rbx, r10 
    cmp ebp, ecx 
    jne .L3 

zobaczyć, jak wykorzystanie o wartości r9 odzwierciedla ciągłą pamięć dla docelowej tablicy i r8 dla jednej z tablic wejściowych.

Na drugim końcu, wektor wektorów generuje kod jak:

.L12: 
    mov r9, QWORD PTR [r12+r11] 
    mov rdi, QWORD PTR [rbx+r11] 
    xor ecx, ecx 
.L16: 
    movss xmm1, DWORD PTR [rdi+rcx] 
    mov rdx, r10 
    xor eax, eax 
    jmp .L15 
.L13: 
    movaps xmm1, xmm0 
.L15: 
    mov rsi, QWORD PTR [rdx] 
    movss xmm0, DWORD PTR [r9+rax] 
    add rax, 4 
    add rdx, 24 
    cmp r8, rax 
    mulss xmm0, DWORD PTR [rsi+rcx] 
    addss xmm0, xmm1 
    movss DWORD PTR [rdi+rcx], xmm0 ; <------------ HERE 
    jne .L13 
    add rcx, 4 
    cmp rcx, r8 
    jne .L16 
    add r11, 24 
    cmp r11, rbp 
    jne .L12 

Nic dziwnego, że kompilator jest na tyle sprytny, aby nie generować kod dla wszystkich połączeń operator [] i robi dobrą robotę inline je, ale zobacz, jak musi śledzić różne adresy poprzez rdi + rcx, gdy przechowuje wartość z powrotem do wektora wyników, a także dodatkowy dostęp do pamięci dla różnych pod-wektorów (mov rsi, QWORD PTR [rdx]), które wszystkie generują pewien narzut.