Thrust Complex Transform de 3 vectores de diferentes tamaños

Question

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

C++ Código

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

Para empezar, algunos comentarios generales. Su 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 norma BLAS gemv operación

donde la matriz se almacena en orden de fila mayor. 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 absolutamente ninguna manera 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 es probablemente va a ser algo como:

1. Uso empuje :: transformar para calcular el (n * n) los productos de los elementos de B y d
2. Uso empuje :: reduce_by_key para reducir los productos en sumas parciales, produciendo Bd
3. uso empuje :: transformada para añadir el producto matriz-vector resultante a v para dar el resultado final.

En el código, definen en primer lugar un funtor como esto:

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

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

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 manera terriblemente ineficiente para hacer el cálculo en comparación con el uso de un código de multiplicación matriz-vector diseñado como el dgemv de CUBLAS.

Answer 2

¿Cuánto difieren sus 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 punto flotante, especialmente aquellas que agregan o multiplican ciertos valores repetidamente, son no asociativo, 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 punto flotante: http://en.wikipedia.org/wiki/Floating_point#Accuracy_problems