Czy dobrym pomysłem jest użycie wektora <vector <double>> do utworzenia klasy macierzowej dla wysokowydajnego naukowego kodu komputerowego?

Czy dobrym pomysłem jest użycie vector<vector<double>>(przy użyciu std) do utworzenia klasy macierzy dla wysokowydajnego naukowego kodu komputerowego?

Jeśli odpowiedź brzmi „nie”. Czemu? Dzięki

To zły pomysł, ponieważ wektor musi przydzielić tyle obiektów w przestrzeni, ile macierzy w macierzy. Alokacja jest droga, ale przede wszystkim jest to zły pomysł, ponieważ dane macierzy istnieją teraz w wielu tablicach rozrzuconych wokół pamięci, a nie w jednym miejscu, w którym pamięć podręczna procesora może z łatwością uzyskać do niej dostęp.

Jest to również marnotrawny format przechowywania: std :: vector przechowuje dwa wskaźniki, jeden na początku tablicy i jeden na końcu, ponieważ długość tablicy jest elastyczna. Z drugiej strony, aby była to właściwa matryca, długości wszystkich rzędów muszą być takie same, więc wystarczyłoby przechować liczbę kolumn tylko raz, zamiast pozwolić, aby każdy wiersz przechowywał swoją długość niezależnie.

Oprócz powodów, o których wspomniał Wolfgang, jeśli używasz vector<vector<double> >, będziesz musiał wyrejestrować go dwukrotnie za każdym razem, gdy chcesz odzyskać element, co jest bardziej kosztowne obliczeniowo niż pojedyncza operacja dereferencji. Jednym typowym podejściem jest zamiast tego przydzielenie jednej tablicy (a vector<double>lub a double *). Widziałem także, jak ludzie dodają cukier składniowy do klas macierzy, owijając wokół tej pojedynczej tablicy kilka bardziej intuicyjnych operacji indeksowania, aby zmniejszyć ilość „mentalnego obciążenia” potrzebnego do wywołania właściwych wskaźników.

Nie, użyj jednej z dostępnych bibliotek algebry liniowej. Dyskusję na temat różnych bibliotek można znaleźć tutaj: Zalecenia dotyczące użytecznej, szybkiej biblioteki macierzy C ++?

Czy to naprawdę takie złe?

@Wolfgang: W zależności od wielkości gęstej matrycy dwa dodatkowe wskaźniki na wiersz mogą być nieistotne. Jeśli chodzi o rozproszone dane, można pomyśleć o użyciu niestandardowego alokatora, który zapewnia, że ​​wektory są w ciągłej pamięci. Tak długo, jak pamięć nie zostanie poddana recyklingowi, nawet standardowy alokator będzie nam przylegał do pamięci z odstępem o wielkości dwóch wskaźników.

@Geoff: Jeśli korzystasz z dostępu losowego i używasz tylko jednej tablicy, nadal musisz obliczyć indeks. Nie może być szybciej.

Zróbmy więc mały test:

vectormatrix.cc:

#include<vector>
#include<iostream>
#include<random>
#include <functional>
#include <sys/time.h>

int main()
{
  int N=1000;
  struct timeval start, end;

  std::cout<< "Checking differenz between last entry of previous row and first entry of this row"<<std::endl;
  std::vector<std::vector<double> > matrix(N, std::vector<double>(N, 0.0));
  for(std::size_t i=1; i<N;i++)
    std::cout<< "index "<<i<<": "<<&(matrix[i][0])-&(matrix[i-1][N-1])<<std::endl;
  std::cout<<&(matrix[0][N-1])<<" "<<&(matrix[1][0])<<std::endl;
  gettimeofday(&start, NULL);
  int k=0;

  for(int j=0; j<100; j++)
    for(std::size_t i=0; i<N;i++)
      for(std::size_t j=0; j<N;j++, k++)
        matrix[i][j]=matrix[i][j]*matrix[i][j];
  gettimeofday(&end, NULL);
  double seconds  = end.tv_sec  - start.tv_sec;
  double useconds = end.tv_usec - start.tv_usec;

  double mtime = ((seconds) * 1000 + useconds/1000.0) + 0.5;

  std::cout<<"calc took: "<<mtime<<" k="<<k<<std::endl;

  std::normal_distribution<double> normal_dist(0, 100);
  std::mt19937 engine; // Mersenne twister MT19937
  auto generator = std::bind(normal_dist, engine);
  for(std::size_t i=1; i<N;i++)
    for(std::size_t j=1; j<N;j++)
      matrix[i][j]=generator();
}

A teraz za pomocą jednej tablicy:

arraymatrix.cc

    #include<vector>
#include<iostream>
#include<random>
#include <functional>
#include <sys/time.h>

int main()
{
  int N=1000;
  struct timeval start, end;

  std::cout<< "Checking difference between last entry of previous row and first entry of this row"<<std::endl;
  double* matrix=new double[N*N];
  for(std::size_t i=1; i<N;i++)
    std::cout<< "index "<<i<<": "<<(matrix+(i*N))-(matrix+(i*N-1))<<std::endl;
  std::cout<<(matrix+N-1)<<" "<<(matrix+N)<<std::endl;

  int NN=N*N;
  int k=0;

  gettimeofday(&start, NULL);
  for(int j=0; j<100; j++)
    for(double* entry =matrix, *endEntry=entry+NN;
        entry!=endEntry;++entry, k++)
      *entry=(*entry)*(*entry);
  gettimeofday(&end, NULL);
  double seconds  = end.tv_sec  - start.tv_sec;
  double useconds = end.tv_usec - start.tv_usec;

  double mtime = ((seconds) * 1000 + useconds/1000.0) + 0.5;

  std::cout<<"calc took: "<<mtime<<" k="<<k<<std::endl;

  std::normal_distribution<double> normal_dist(0, 100);
  std::mt19937 engine; // Mersenne twister MT19937
  auto generator = std::bind(normal_dist, engine);
  for(std::size_t i=1; i<N*N;i++)
      matrix[i]=generator();
}

W moim systemie jest teraz wyraźny zwycięzca (kompilator gcc 4.7 z -O3)

odbitki vectormatrix czasu:

index 997: 3
index 998: 3
index 999: 3
0xc7fc68 0xc7fc80
calc took: 185.507 k=100000000

real    0m0.257s
user    0m0.244s
sys     0m0.008s

Widzimy również, że dopóki standardowy alokator nie przetwarza zwolnionej pamięci, dane są ciągłe. (Oczywiście po niektórych zwolnieniach nie ma na to gwarancji.)

odbitki macierzy czasu:

index 997: 1
index 998: 1
index 999: 1
0x7ff41f208f48 0x7ff41f208f50
calc took: 187.349 k=100000000

real    0m0.257s
user    0m0.248s
sys     0m0.004s

Nie polecam tego, ale nie z powodu problemów z wydajnością. Będzie on nieco mniej wydajny niż tradycyjna macierz, która zwykle jest alokowana jako duża część ciągłych danych, które są indeksowane za pomocą pojedynczego wskaźnika dereferencji i arytmetyki liczb całkowitych. Powodem spadku wydajności są głównie różnice w buforowaniu, ale gdy rozmiar macierzy wystarczająco się powiększy, efekt zostanie amortyzowany, a jeśli użyjesz specjalnego alokatora dla wewnętrznych wektorów, aby były one wyrównane do granic pamięci podręcznej, to dodatkowo złagodzi problem buforowania .

Moim zdaniem nie jest to wystarczający powód, aby tego nie robić. Powodem dla mnie jest to, że powoduje to wiele bólów głowy związanych z kodowaniem. Oto lista bólów głowy, które spowodują w dłuższej perspektywie

Korzystanie z bibliotek HPC

Jeśli chcesz korzystać z większości bibliotek HPC, musisz iterować wektor i umieścić wszystkie ich dane w ciągłym buforze, ponieważ większość bibliotek HPC oczekuje tego jawnego formatu. Przychodzą mi na myśl BLAS i LAPACK, ale także wszechobecna biblioteka HPC MPI byłaby znacznie trudniejsza w użyciu.

Większy potencjał błędu kodowania

std::vectornic nie wie o swoich wpisach. Jeśli wypełnisz std::vectorwięcej std::vectors, Twoim zadaniem jest upewnić się, że wszystkie mają ten sam rozmiar, ponieważ pamiętaj, że chcemy, aby macierz i macierze nie miały zmiennej liczby wierszy (lub kolumn). Dlatego będziesz musiał wywoływać wszystkie poprawne konstruktory dla każdego wpisu zewnętrznego wektora, a każdy, kto używa twojego kodu, musi oprzeć się pokusie użycia std::vector<T>::push_back()dowolnego z wewnętrznych wektorów, co spowodowałoby uszkodzenie całego następującego kodu. Oczywiście możesz tego zabronić, jeśli piszesz poprawnie swoją klasę, ale o wiele łatwiej jest egzekwować to po prostu z dużym ciągłym przydziałem.

Kultura i oczekiwania HPC

Programiści HPC po prostu oczekują danych niskiego poziomu. Jeśli dasz im matrycę, oczekuje się, że jeśli złapią wskaźnik do pierwszego elementu macierzy i wskaźnik do ostatniego elementu macierzy, wówczas wszystkie wskaźniki pomiędzy tymi dwoma są prawidłowe i wskazują na elementy tego samego matryca. Jest to podobne do mojego pierwszego punktu, ale inne, ponieważ może nie być tak bardzo powiązane z bibliotekami, ale raczej z członkami zespołu lub kimkolwiek, komu udostępniasz swój kod.

Łatwiej jest uzasadnić wydajność danych niższego poziomu

Zrzucenie na najniższy poziom reprezentacji pożądanej struktury danych ułatwia życie na dłuższą metę dla HPC. Korzystanie z narzędzi takich jak perfi vtunezapewni ci bardzo niski poziom liczników wydajności, które spróbujesz połączyć z tradycyjnymi wynikami profilowania w celu poprawy wydajności twojego kodu. Jeśli struktura danych wykorzystuje wiele fantazyjnych kontenerów, trudno będzie zrozumieć, że brak pamięci podręcznej wynika z problemu z kontenerem lub z nieefektywności samego algorytmu. W przypadku bardziej skomplikowanych kodów potrzebne są pojemniki, ale w przypadku algebry macierzowej tak naprawdę nie są - możesz sobie radzić z 1 std::vectorprzechowywaniem danych zamiast n std::vectors, więc idź z tym.

Piszę również punkt odniesienia. W przypadku matrycy o małym rozmiarze (<100 * 100) wydajność jest podobna dla wektora <wektor <podwójny >> i owiniętego wektora 1D. W przypadku matrycy o dużych rozmiarach (~ 1000 * 1000) lepiej jest owinięty wektor 1D. Macierz własna zachowuje się gorzej. Dziwi mnie, że Eigen jest najgorszy.

#include <iostream>
#include <iomanip>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <map>
#include <vector>
#include <string>
#include <cmath>
#include <numeric>
#include "time.h"
#include <chrono>
#include <cstdlib>
#include <Eigen/Dense>

using namespace std;
using namespace std::chrono;    // namespace for recording running time
using namespace Eigen;

int main()
{
    const int row = 1000;
    const int col = row;
    const int N = 1e8;

    // 2D vector
    auto start = high_resolution_clock::now();
    vector<vector<double>> vec_2D(row,vector<double>(col,0.));
    for (int i = 0; i < N; i++)
    {
        for (int i=0; i<row; i++)
        {
            for (int j=0; j<col; j++)
            {
                vec_2D[i][j] *= vec_2D[i][j];
            }
        }
    }
    auto stop = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(stop - start);
    cout << "2D vector: " << duration.count()/1e6 << " s" << endl;

    // 2D array
    start = high_resolution_clock::now();
    double array_2D[row][col];
    for (int i = 0; i < N; i++)
    {
        for (int i=0; i<row; i++)
        {
            for (int j=0; j<col; j++)
            {
                array_2D[i][j] *= array_2D[i][j];
            }
        }
    }
    stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start);
    cout << "2D array: " << duration.count() / 1e6 << " s" << endl;

    // wrapped 1D vector
    start = high_resolution_clock::now();
    vector<double> vec_1D(row*col, 0.);
    for (int i = 0; i < N; i++)
    {
        for (int i=0; i<row; i++)
        {
            for (int j=0; j<col; j++)
            {
                vec_1D[i*col+j] *= vec_1D[i*col+j];
            }
        }
    }
    stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start);
    cout << "1D vector: " << duration.count() / 1e6 << " s" << endl;

    // eigen 2D matrix
    start = high_resolution_clock::now();
    MatrixXd mat(row, col);
    for (int i = 0; i < N; i++)
    {
        for (int j=0; j<col; j++)
        {
            for (int i=0; i<row; i++)
            {
                mat(i,j) *= mat(i,j);
            }
        }
    }
    stop = high_resolution_clock::now();
    duration = duration_cast<microseconds>(stop - start);
    cout << "2D eigen matrix: " << duration.count() / 1e6 << " s" << endl;
}

Jak zauważyli inni, nie próbuj z tym robić matematyki ani nie rób nic wydajnego.

To powiedziawszy, użyłem tej struktury jako tymczasowej, gdy kod musi złożyć tablicę 2-D, której wymiary zostaną określone w czasie wykonywania i po rozpoczęciu przechowywania danych. Na przykład, zbieranie danych wyjściowych wektora z jakiegoś kosztownego procesu, w którym nie jest łatwo dokładnie obliczyć, ile wektorów trzeba przechowywać przy starcie.

Możesz po prostu połączyć wszystkie swoje dane wektorowe w jednym buforze, gdy tylko się pojawią, ale kod będzie bardziej trwały i bardziej czytelny, jeśli użyjesz vector<vector<T>>.