Apakah ide yang baik untuk menggunakan vektor <vektor <double>> untuk membentuk kelas matriks untuk kode komputasi ilmiah kinerja tinggi?

Apakah ide yang baik untuk menggunakan vector<vector<double>>(menggunakan std) untuk membentuk kelas matriks untuk kode komputasi ilmiah kinerja tinggi?

Jika jawabannya tidak. Mengapa? Terima kasih

Itu ide yang buruk karena vektor perlu mengalokasikan objek sebanyak mungkin di ruang karena ada baris dalam matriks Anda. Alokasi mahal, tetapi terutama itu adalah ide yang buruk karena data matriks Anda sekarang ada di sejumlah array yang tersebar di sekitar memori, daripada semua di satu tempat di mana cache prosesor dapat dengan mudah mengaksesnya.

Ini juga merupakan format penyimpanan yang boros: std :: vector menyimpan dua pointer, satu ke awal array dan satu lagi ke akhir karena panjang array fleksibel. Di sisi lain, agar ini menjadi matriks yang tepat, panjang semua baris harus sama sehingga cukup untuk menyimpan jumlah kolom hanya sekali, daripada membiarkan setiap baris menyimpan panjangnya secara independen.

Selain alasan yang disebutkan Wolfgang, jika Anda menggunakan a vector<vector<double> >, Anda harus melakukan dereferensi dua kali setiap kali Anda ingin mengambil elemen, yang secara komputasi lebih mahal daripada operasi dereferencing tunggal. Salah satu pendekatan yang umum adalah mengalokasikan satu array (a vector<double>atau a double *) sebagai gantinya. Saya juga melihat orang menambahkan gula sintaksis ke kelas matriks dengan membungkus array tunggal ini beberapa operasi pengindeksan yang lebih intuitif, untuk mengurangi jumlah "overhead mental" yang diperlukan untuk memohon indeks yang tepat.

Tidak, gunakan salah satu perpustakaan aljabar linear gratis yang tersedia. Diskusi tentang berbagai pustaka dapat ditemukan di sini: Rekomendasi untuk pustaka matriks C ++ yang dapat digunakan dan cepat?

Apakah ini benar-benar hal yang buruk?

@ Wolfgang: Tergantung pada ukuran matriks padat, dua pointer tambahan per baris mungkin dapat diabaikan. Mengenai data yang tersebar, seseorang dapat berpikir untuk menggunakan pengalokasi kustom yang memastikan bahwa vektor berada dalam memori yang berdekatan. Selama memori tidak didaur ulang, bahkan pengalokasi standar akan menggunakan memori yang berdekatan dengan celah ukuran dua penunjuk.

@ Geoff: Jika Anda melakukan akses acak dan menggunakan hanya satu array Anda masih harus menghitung indeks. Mungkin tidak akan lebih cepat.

Jadi mari kita lakukan tes kecil:

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();
}

Dan sekarang menggunakan satu array:

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();
}

Di sistem saya sekarang ada pemenang yang jelas (Compiler gcc 4.7 dengan -O3)

cetakan waktu vectormatrix:

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

Kita juga melihat, bahwa selama pengalokasi standar tidak mendaur ulang memori yang dibebaskan, data tersebut bersebelahan. (Tentu saja setelah beberapa deallocations tidak ada jaminan untuk ini.)

cetakan arraymatrix:

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

Saya tidak merekomendasikannya, tetapi bukan karena masalah kinerja. Ini akan menjadi sedikit kurang berkinerja daripada matriks tradisional, yang biasanya dialokasikan sebagai bagian besar dari data yang berdekatan yang diindeks menggunakan dereference pointer tunggal dan aritmatika integer. Alasan untuk hit kinerja sebagian besar perbedaan caching, tetapi setelah ukuran matriks Anda cukup besar efek ini akan diamortisasi dan jika Anda menggunakan pengalokasi khusus untuk vektor bagian dalam sehingga mereka selaras dengan batas-batas cache maka ini lebih lanjut mengurangi masalah caching .

Dengan sendirinya itu bukan alasan yang cukup untuk tidak melakukannya, menurut saya. Alasan saya adalah membuat banyak sakit kepala kode. Berikut daftar sakit kepala yang akan ditimbulkan dalam jangka panjang

Penggunaan perpustakaan HPC

Jika Anda ingin menggunakan sebagian besar pustaka HPC Anda harus mengulangi vektor Anda dan menempatkan semua datanya dalam buffer yang bersebelahan, karena sebagian besar pustaka HPC mengharapkan format eksplisit ini. BLAS dan LAPACK datang ke pikiran, tetapi juga MPI perpustakaan HPC mana-mana akan jauh lebih sulit untuk digunakan.

Lebih banyak potensi kesalahan pengkodean

std::vectortidak tahu apa-apa tentang entri mereka. Jika Anda mengisi std::vectordengan lebih dari std::vectors maka itu sepenuhnya tugas Anda untuk memastikan bahwa mereka semua memiliki ukuran yang sama, karena ingat bahwa kami ingin matriks dan matriks tidak memiliki jumlah baris (atau kolom) yang bervariasi. Dengan demikian Anda harus memanggil semua konstruktor yang benar untuk setiap entri vektor luar Anda, dan siapa pun yang menggunakan kode Anda harus menahan godaan untuk digunakan std::vector<T>::push_back()pada salah satu vektor bagian dalam, yang akan menyebabkan semua kode berikut rusak. Tentu saja Anda dapat melarang ini jika Anda menulis kelas dengan benar, tetapi jauh lebih mudah untuk menegakkan ini hanya dengan alokasi bersebelahan besar.

Budaya dan harapan HPC

Programer HPC hanya mengharapkan data level rendah. Jika Anda memberi mereka sebuah matriks, ada harapan bahwa jika mereka meraih pointer ke elemen pertama dari matriks dan sebuah pointer ke elemen terakhir dari matriks, maka semua pointer di antara keduanya valid dan arahkan ke elemen yang sama. matriks. Ini mirip dengan poin pertama saya, tetapi berbeda karena mungkin tidak terkait banyak dengan perpustakaan tetapi anggota tim atau siapa pun yang Anda bagikan kode Anda.

Lebih mudah untuk alasan tentang kinerja data tingkat yang lebih rendah

Menurunkan ke tingkat representasi terendah dari struktur data yang Anda inginkan menjadikan hidup Anda lebih mudah dalam jangka panjang untuk HPC. Menggunakan alat-alat seperti perfdan vtuneakan memberi Anda pengukuran penghitung kinerja tingkat sangat rendah yang akan Anda coba gabungkan dengan hasil profil tradisional untuk meningkatkan kinerja kode Anda. Jika struktur data Anda menggunakan banyak wadah mewah, akan sulit untuk memahami bahwa kesalahan cache berasal dari masalah dengan wadah atau ketidakefisienan dalam algoritma itu sendiri. Diperlukan wadah kode yang lebih rumit, tetapi untuk aljabar matriks sebenarnya tidak - Anda bisa bertahan hanya 1 std::vectordengan menyimpan data daripada n std::vectors, jadi ikuti saja.

Saya juga menulis patokan. Untuk matriks ukuran kecil (<100 * 100), kinerjanya mirip untuk vektor <vektor <ganda >> dan membungkus vektor 1D. Untuk matriks ukuran besar (~ 1000 * 1000), vektor 1D yang dibungkus lebih baik. Matriks Eigen berperilaku lebih buruk. Sangat mengejutkan bagi saya bahwa Eigen adalah yang terburuk.

#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;
}

Seperti yang telah ditunjukkan orang lain, jangan mencoba melakukan matematika dengan itu atau melakukan pemain apa pun.

Yang mengatakan, saya telah menggunakan struktur ini sebagai sementara ketika kode perlu merakit array 2-D yang dimensinya akan ditentukan pada saat runtime dan setelah Anda mulai menyimpan data. Misalnya, mengumpulkan keluaran vektor dari beberapa proses mahal di mana tidak mudah untuk menghitung dengan tepat berapa banyak vektor yang perlu Anda simpan saat startup.

Anda bisa menggabungkan semua input vektor Anda menjadi satu buffer saat mereka masuk, tetapi kode akan lebih tahan lama dan lebih mudah dibaca jika Anda menggunakan a vector<vector<T>>.