É uma boa ideia usar o vetor <vetor <duplo>> para formar uma classe de matriz para código de computação científica de alto desempenho?

É uma boa idéia usar vector<vector<double>>(usando std) para formar uma classe de matriz para código de computação científica de alto desempenho?

Se a resposta for não. Por quê? obrigado

É uma péssima idéia, porque o vetor precisa alocar tantos objetos no espaço quanto houver linhas em sua matriz. A alocação é cara, mas principalmente é uma má idéia, porque os dados de sua matriz agora existem em várias matrizes espalhadas pela memória, em vez de em um único local onde o cache do processador pode acessá-las facilmente.

Também é um formato de armazenamento desnecessário: std :: vector armazena dois ponteiros, um para o início da matriz e outro para o final, porque o comprimento da matriz é flexível. Por outro lado, para que essa seja uma matriz adequada, os comprimentos de todas as linhas devem ser os mesmos e, portanto, seria suficiente armazenar o número de colunas apenas uma vez, em vez de permitir que cada linha armazene seu comprimento independentemente.

Além das razões mencionadas por Wolfgang, se você usar a vector<vector<double> >, precisará desdiferenciá-lo duas vezes toda vez que quiser recuperar um elemento, que é mais caro em termos computacionais do que uma única operação de desreferenciação. Uma abordagem típica é alocar uma única matriz (a vector<double>ou a double *). Também vi pessoas adicionando açúcar sintático às classes matriciais, envolvendo em torno dessa matriz única algumas operações de indexação mais intuitivas, para reduzir a quantidade de "sobrecarga mental" necessária para invocar os índices adequados.

Não, use uma das bibliotecas de álgebra linear disponíveis gratuitamente. Uma discussão sobre diferentes bibliotecas pode ser encontrada aqui: Recomendações para uma biblioteca de matrizes C ++ rápida e utilizável?

É realmente uma coisa tão ruim?

@ Wolfgang: Dependendo do tamanho da matriz densa, dois ponteiros adicionais por linha podem ser desprezíveis. Em relação aos dados dispersos, pode-se pensar em usar um alocador personalizado que garanta que os vetores estejam na memória contígua. Enquanto a memória não for reciclada, mesmo o alocador padrão usará memória contígua com uma diferença de tamanho de dois ponteiros.

@ Geoff: Se você está fazendo acesso aleatório e usa apenas uma matriz, ainda precisa calcular o índice. Pode não ser mais rápido.

Então, vamos fazer um pequeno teste:

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

E agora usando uma matriz:

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

No meu sistema agora existe um vencedor claro (Compilador gcc 4.7 com -O3)

a matriz de vetor de tempo imprime:

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

Também vemos que, enquanto o alocador padrão não reciclar a memória liberada, os dados serão contíguos. (Obviamente, após algumas desalocações, não há garantia para isso.)

a matriz de tempo imprime:

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

Não recomendo, mas não por problemas de desempenho. Será um pouco menos eficiente do que uma matriz tradicional, que geralmente é alocada como uma grande porção de dados contíguos indexados usando uma única desreferência de ponteiro e aritmética de número inteiro. O motivo do impacto no desempenho são principalmente as diferenças de armazenamento em cache, mas quando o tamanho da matriz ficar grande o suficiente, esse efeito será amortizado e se você usar um alocador especial para os vetores internos, para que eles estejam alinhados aos limites do cache, isso atenua ainda mais o problema de armazenamento em cache. .

Isso por si só não é motivo suficiente para não fazê-lo, na minha opinião. A razão para mim é que isso cria muitas dores de cabeça de codificação. Aqui está uma lista de dores de cabeça que isso causará a longo prazo

Uso de bibliotecas HPC

Se você deseja usar a maioria das bibliotecas HPC, precisará iterar sobre seu vetor e colocar todos os dados em um buffer contíguo, porque a maioria das bibliotecas HPC espera esse formato explícito. BLAS e LAPACK vêm à mente, mas também a onipresente biblioteca HPC MPI seria muito mais difícil de usar.

Mais potencial para erro de codificação

std::vectornão sabe nada sobre suas entradas. Se você preencher um std::vectorcom mais std::vectors, é inteiramente seu trabalho garantir que todos tenham o mesmo tamanho, porque lembre-se de que queremos uma matriz e as matrizes não têm número variável de linhas (ou colunas). Portanto, você precisará chamar todos os construtores corretos para cada entrada do seu vetor externo, e qualquer pessoa que use seu código deve resistir à tentação de usar std::vector<T>::push_back()em qualquer um dos vetores internos, o que causaria a quebra de todo o código a seguir. É claro que você pode proibir isso se escrever sua classe corretamente, mas é muito mais fácil aplicar isso simplesmente com uma grande alocação contígua.

Cultura e expectativas da HPC

Os programadores de HPC simplesmente esperam dados de baixo nível. Se você lhes der uma matriz, existe a expectativa de que, se eles pegaram o ponteiro para o primeiro elemento da matriz e um ponteiro para o último elemento da matriz, todos os ponteiros entre esses dois são válidos e apontam para elementos do mesmo matriz. Isso é semelhante ao meu primeiro ponto, mas diferente porque pode não estar muito relacionado às bibliotecas, mas aos membros da equipe ou a qualquer pessoa com quem você compartilhe seu código.

Mais fácil de raciocinar sobre o desempenho de dados de nível inferior

Passar para a representação de nível mais baixo da estrutura de dados desejada facilita sua vida a longo prazo para o HPC. O uso de ferramentas como perfe vtunefornecerá medições de contador de desempenho de nível muito baixo, que você tentará combinar com os resultados de criação de perfil tradicionais para melhorar o desempenho do seu código. Se sua estrutura de dados usa muitos contêineres sofisticados, será difícil entender que as falhas de cache são decorrentes de um problema com o contêiner ou de uma ineficiência no próprio algoritmo. Para contêineres de código mais complicados, são necessários, mas para álgebra matricial eles realmente não são - você pode conviver apenas 1 std::vectorpara armazenar os dados em vez de n std::vectors, então vá em frente.

Eu também escrevo uma referência. Para matrizes de tamanho pequeno (<100 * 100), o desempenho é semelhante para o vetor <vetor <duplo >> e o vetor 1D empacotado. Para matrizes de tamanho grande (~ 1000 * 1000), o vetor 1D empacotado é melhor. A matriz Eigen se comporta pior. Surpreende-me que o Eigen seja o pior.

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

Como outros já apontaram, não tente fazer contas com ela nem faça nada com desempenho.

Dito isso, usei essa estrutura como temporária quando o código precisa montar uma matriz 2-D cujas dimensões serão determinadas no tempo de execução e depois que você começar a armazenar dados. Por exemplo, coletando saídas vetoriais de algum processo caro, em que não é simples calcular exatamente quantos vetores você precisará armazenar na inicialização.

Você pode concatenar todas as suas entradas de vetor em um buffer à medida que elas forem chegando, mas o código será mais durável e mais legível se você usar a vector<vector<T>>.