cuda - Transformación compleja de empuje de 3 vectores de diferentes tamaños

thrust

Question

Hola, tengo este ciclo en C +, y estaba tratando de convertirlo en empuje pero sin obtener los mismos resultados ... ¿Alguna idea? gracias

Código C ++

for (i=0;i<n;i++) 
    for (j=0;j<n;j++) 
      values[i]=values[i]+(binv[i*n+j]*d[j]);

Código de empuje

thrust::fill(values.begin(), values.end(), 0);
thrust::transform(make_zip_iterator(make_tuple(
                thrust::make_permutation_iterator(values.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator(0), IndexDivFunctor(n))),
                binv.begin(),
                thrust::make_permutation_iterator(d.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator(0), IndexModFunctor(n))))),
                make_zip_iterator(make_tuple(
                thrust::make_permutation_iterator(values.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator(0), IndexDivFunctor(n))) + n,
                binv.end(),
                thrust::make_permutation_iterator(d.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator(0), IndexModFunctor(n))) + n)),
                thrust::make_permutation_iterator(values.begin(), thrust::make_transform_iterator(thrust::make_counting_iterator(0), IndexDivFunctor(n))),
                function1()
                );

Funciones de empuje

struct IndexDivFunctor: thrust::unary_function<int, int>
{
  int n;

  IndexDivFunctor(int n_) : n(n_) {}

  __host__ __device__
  int operator()(int idx)
  {
    return idx / n;
  }
};

struct IndexModFunctor: thrust::unary_function<int, int>
{
  int n;

  IndexModFunctor(int n_) : n(n_) {}

  __host__ __device__
  int operator()(int idx)
  {
    return idx % n;
  }
};


struct function1
{
  template <typename Tuple>
  __host__ __device__
  double operator()(Tuple v)
  {
    return thrust::get<0>(v) + thrust::get<1>(v) * thrust::get<2>(v);
  }
};

Answer 1

¿Cuánto difieren tus resultados? ¿Es una respuesta completamente diferente, o solo difiere en los últimos dígitos? ¿El ciclo se ejecuta solo una vez o es algún tipo de proceso iterativo?

Las operaciones de coma flotante, especialmente aquellas que agregan o multiplican ciertos valores, no son asociativas debido a problemas de precisión. Además, si usa optimizaciones matemáticas rápidas, las operaciones pueden no ser compiladas por IEEE.

Para empezar, echa un vistazo a esta sección de wikipedia sobre los números de coma flotante: http://en.wikipedia.org/wiki/Floating_point#Accuracy_problems

Answer 2

Para empezar, algunos comentarios generales. Tu bucle

for (i=0;i<n;i++) 
    for (j=0;j<n;j++) 
      v[i]=v[i]+(B[i*n+j]*d[j]);

es el equivalente de la operación estándar de gemv de BLAS

donde la matriz se almacena en orden mayor fila. La forma óptima de hacer esto en el dispositivo sería usar CUBLAS, no algo construido a partir de primitivas de empuje.

Habiendo dicho eso, no hay forma de que el código de empuje que publicas vaya a hacer lo que hace tu código de serie. Los errores que está viendo no son el resultado de la asociatividad de coma flotante. Fundamentalmente thrust::transform aplica el functor suministrado a cada elemento del iterador de entrada y almacena el resultado en el iterador de salida. Para obtener el mismo resultado que el bucle que publicó, la llamada thrust::transform necesitaría realizar (n * n) operaciones del fmad functor que publicó. Claramente no. Además, no hay garantía de que thrust::transform realice la operación de suma / reducción de una manera que estaría a salvo de las carreras de memoria.

La solución correcta probablemente será algo así como:

1. Utilice thrust :: transform para calcular los productos (n * n) de los elementos de B y D
2. Utilice thrust :: reduce_by_key para reducir los productos en sumas parciales, produciendo Bd
3. Utilice thrust :: transform para agregar el producto matriz-vector resultante a v para obtener el resultado final.

En el código, primero defina un functor como este:

struct functor
{
  template <typename Tuple>
  __host__ __device__
  double operator()(Tuple v)
  {
    return thrust::get<0>(v) * thrust::get<1>(v);
  }
};

Luego, haga lo siguiente para calcular la multiplicación matriz-vector

  typedef thrust::device_vector<int> iVec;
  typedef thrust::device_vector<double> dVec;

  typedef thrust::counting_iterator<int> countIt;
  typedef thrust::transform_iterator<IndexDivFunctor, countIt> columnIt;
  typedef thrust::transform_iterator<IndexModFunctor, countIt> rowIt;

  // Assuming the following allocations on the device
  dVec B(n*n), v(n), d(n);

  // transformation iterators mapping to vector rows and columns
  columnIt cv_begin = thrust::make_transform_iterator(thrust::make_counting_iterator(0), IndexDivFunctor(n));
  columnIt cv_end   = cv_begin + (n*n);

  rowIt rv_begin = thrust::make_transform_iterator(thrust::make_counting_iterator(0), IndexModFunctor(n));
  rowIt rv_end   = rv_begin + (n*n);

  dVec temp(n*n);
  thrust::transform(make_zip_iterator(
                      make_tuple(
                        B.begin(),
                        thrust::make_permutation_iterator(d.begin(),rv_begin) ) ),
                    make_zip_iterator(
                      make_tuple(
                        B.end(),
                        thrust::make_permutation_iterator(d.end(),rv_end) ) ),
                    temp.begin(),
                    functor());

  iVec outkey(n);
  dVec Bd(n);
  thrust::reduce_by_key(cv_begin, cv_end, temp.begin(), outkey.begin(), Bd.begin());
  thrust::transform(v.begin(), v.end(), Bd.begin(), v.begin(), thrust::plus<double>());

Por supuesto, esta es una forma terriblemente ineficiente de hacer el cálculo en comparación con el uso de un código de multiplicación matriz-vector diseñado como dgemv de CUBLAS.