Traslazione complessa spinta di 3 vettori di dimensioni diverse

Question

Ciao ho questo ciclo in C +, e stavo cercando di convertirlo in spinta ma senza ottenere gli stessi risultati ... Qualche idea? grazie

codice C++

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

Codice spinta

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

Funzioni di spinta

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

Per cominciare, alcune osservazioni di carattere generale. Il ciclo

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

è l'equivalente dello standard BLAS gemv operazione

dove la matrice è memorizzato nella riga ordine maggiore. Il modo ottimale per farlo sul dispositivo sarebbe utilizzare CUBLAS, non qualcosa costruito con primitive di spinta.

Detto questo, non c'è assolutamente alcun modo in cui il codice di spinta che hai postato farà mai quello che fa il tuo codice seriale. Gli errori che stai vedendo non sono il risultato dell'associatività in virgola mobile. Fondamentalmente, thrust::transform applica il functor fornito a ogni elemento dell'iteratore di input e memorizza il risultato sull'iteratore di output. Per ottenere lo stesso risultato del ciclo che hai postato, la chiamata thrust::transform dovrebbe eseguire (n * n) le operazioni del fmad functor che hai postato. Chiaramente no. Inoltre, non vi è alcuna garanzia che lo thrust::transform esegua l'operazione di sommatoria/riduzione in un modo che sarebbe al sicuro dalle gare della memoria.

La soluzione corretta è destinata probabilmente ad essere qualcosa come:

1. Usa spinta :: trasformare per calcolare la (n * n) i prodotti degli elementi di B e d
2. Usa spinta :: reduce_by_key per ridurre i prodotti nella somme parziali, cedendo Bd
3. Uso spinta :: trasformare aggiungere il prodotto matrice-vettore risultante a v per ottenere il risultato finale.

Nel codice, innanzitutto definire un funtore simili:

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

Quindi fare quanto segue per calcolare la matrice-vettore moltiplicazione

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

Naturalmente, questo è un modo terribilmente inefficiente per eseguire il calcolo rispetto all'utilizzo di un codice di moltiplicazione di matrice vettoriale progettato allo scopo come dgemv da CUBLAS.

Answer 2

Quanto differiscono i risultati? È una risposta completamente diversa o si differenzia solo sulle ultime cifre? Il ciclo viene eseguito solo una volta o è una sorta di processo iterativo?

Le operazioni in virgola mobile, in particolare quelle che sommano ripetutamente o moltiplicano determinati valori, sono non associative, a causa di problemi di precisione. Inoltre, se si utilizzano ottimizzazioni matematiche veloci, le operazioni potrebbero non essere un compilatore IEEE.

Per cominciare, consulta questa sezione wikipedia su numeri in virgola mobile: http://en.wikipedia.org/wiki/Floating_point#Accuracy_problems