¿Es una buena idea usar el vector <vector <double>> para formar una clase de matriz para el código informático científico de alto rendimiento?

¿Es una buena idea usar vector<vector<double>>(usando std) para formar una clase de matriz para el código de computación científica de alto rendimiento?

Si la respuesta es no. ¿Por qué? Gracias

Es una mala idea porque el vector necesita asignar tantos objetos en el espacio como filas en su matriz. La asignación es costosa, pero principalmente es una mala idea porque los datos de su matriz ahora existen en una serie de arreglos dispersos por la memoria, en lugar de todos en un lugar donde el caché del procesador puede acceder fácilmente.

También es un formato de almacenamiento derrochador: std :: vector almacena dos punteros, uno al principio de la matriz y otro al final porque la longitud de la matriz es flexible. Por otro lado, para que esta sea una matriz adecuada, las longitudes de todas las filas deben ser las mismas y, por lo tanto, sería suficiente almacenar el número de columnas solo una vez, en lugar de permitir que cada fila almacene su longitud de forma independiente.

Además de las razones que mencionó Wolfgang, si usa un vector<vector<double> >, deberá desreferenciarlo dos veces cada vez que desee recuperar un elemento, que es más costoso desde el punto de vista computacional que una sola operación de desreferenciación. Un enfoque típico es asignar una matriz única (ao vector<double>a double *) en su lugar. También he visto a personas agregar azúcar sintáctico a las clases de matriz al envolver alrededor de esta matriz única algunas operaciones de indexación más intuitivas, para reducir la cantidad de "sobrecarga mental" necesaria para invocar los índices adecuados.

No, use una de las bibliotecas de álgebra lineal disponibles gratuitamente. Aquí se puede encontrar una discusión sobre las diferentes bibliotecas: ¿ Recomendaciones para una biblioteca de matriz C ++ rápida y utilizable?

¿Es realmente tan malo?

@Wolfgang: Dependiendo del tamaño de la matriz densa, dos punteros adicionales por fila pueden ser insignificantes. Con respecto a los datos dispersos, se podría pensar en usar un asignador personalizado que se asegure de que los vectores estén en memoria contigua. Siempre que la memoria no se recicle, incluso el asignador estándar utilizará memoria contigua con un espacio de dos punteros.

@Geoff: si está haciendo acceso aleatorio y usa solo una matriz, aún tiene que calcular el índice. Puede que no sea más rápido.

Entonces, hagamos una pequeña prueba:

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

Y ahora usando una 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();
}

En mi sistema ahora hay un claro ganador (Compilador gcc 4.7 con -O3)

impresiones de vectormatrix de tiempo:

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

También vemos que, mientras el asignador estándar no recicle la memoria liberada, los datos son contiguos. (Por supuesto, después de algunas desasignaciones no hay garantía para esto).

impresiones de matriz de tiempo:

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

No lo recomiendo, pero no por problemas de rendimiento. Será un poco menos eficaz que una matriz tradicional, que generalmente se asigna como una gran porción de datos contiguos que se indexan utilizando una única referencia de puntero y aritmética de enteros. La razón del impacto en el rendimiento es principalmente las diferencias de almacenamiento en caché, pero una vez que el tamaño de su matriz sea lo suficientemente grande, este efecto se amortizará y si usa un asignador especial para los vectores internos para que estén alineados con los límites del caché, esto mitiga aún más el problema de almacenamiento en caché .

Eso en sí mismo no es razón suficiente para no hacerlo, en mi opinión. La razón para mí es que crea muchos dolores de cabeza de codificación. Aquí hay una lista de dolores de cabeza que esto causará a largo plazo

Uso de bibliotecas HPC

Si desea utilizar la mayoría de las bibliotecas de HPC, deberá iterar sobre su vector y colocar todos sus datos en un búfer contiguo, porque la mayoría de las bibliotecas de HPC esperan este formato explícito. Me vienen a la mente BLAS y LAPACK, pero también la ubicua biblioteca HPC MPI sería mucho más difícil de usar.

Más potencial para codificar errores

std::vectorno sabe nada de sus entradas. Si llena una std::vectorcon más std::vectors, entonces es completamente su trabajo asegurarse de que todas tengan el mismo tamaño, porque recuerde que queremos una matriz y las matrices no tienen un número variable de filas (o columnas). Por lo tanto, tendrá que llamar a todos los constructores correctos para cada entrada de su vector externo, y cualquier otra persona que use su código debe resistir la tentación de usar std::vector<T>::push_back()cualquiera de los vectores internos, lo que provocaría que se rompa todo el código siguiente. Por supuesto, puede rechazar esto si escribe su clase correctamente, pero es mucho más fácil aplicar esto simplemente con una gran asignación contigua.

Cultura y expectativas de HPC

Los programadores de HPC simplemente esperan datos de bajo nivel. Si les da una matriz, existe la expectativa de que si agarran el puntero al primer elemento de la matriz y un puntero al último elemento de la matriz, todos los punteros entre estos dos son válidos y apuntan a elementos de ese mismo matriz. Esto es similar a mi primer punto, pero diferente porque puede no estar relacionado tanto con las bibliotecas sino con los miembros del equipo o cualquier persona con la que comparta su código.

Es más fácil razonar sobre el rendimiento de los datos de nivel inferior

Bajar al nivel más bajo de representación de su estructura de datos deseada hace que su vida sea más fácil a largo plazo para HPC. El uso de herramientas como perfy vtunele proporcionará mediciones de contador de rendimiento de muy bajo nivel que intentará combinar con los resultados de creación de perfiles tradicionales para mejorar el rendimiento de su código. Si su estructura de datos utiliza muchos contenedores sofisticados, será difícil entender que las fallas de caché provienen de un problema con el contenedor o de una ineficiencia en el algoritmo mismo. Para los contenedores de código más complicados son necesarios, pero para el álgebra matricial realmente no lo son: puede arreglárselas solo 1 std::vectorpara almacenar los datos en lugar de n std::vectors, así que vaya con eso.

También escribo un punto de referencia. Para la matriz de tamaño pequeño (<100 * 100), el rendimiento es similar para el vector <vector <double >> y el vector 1D envuelto. Para una matriz de gran tamaño (~ 1000 * 1000), el vector 1D envuelto es mejor. La matriz de Eigen se comporta peor. Me sorprende que el Eigen sea el peor.

#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 otros han señalado, no intentes hacer cálculos matemáticos con él ni hacer nada performante.

Dicho esto, he usado esta estructura como temporal cuando el código necesita ensamblar una matriz 2-D cuyas dimensiones se determinarán en tiempo de ejecución y después de que haya comenzado a almacenar datos. Por ejemplo, recopilar salidas de vectores de algún proceso costoso en el que no es simple calcular exactamente cuántos vectores necesitará almacenar al inicio.

Podrías concatenar todas tus entradas de vectores en un búfer a medida que entran, pero el código será más duradero y más legible si utilizas a vector<vector<T>>.