Как вычислить внутреннее произведение сложных векторов, используя кубы или тяги? - PullRequest
Новая
       23

Как вычислить внутреннее произведение сложных векторов, используя кубы или тяги?

0 голосов
/ 26 апреля 2019

После длительного поиска я все еще не могу решить эту проблему.

У меня есть два вектора: x = [a1, ..., aN], y = [b1, ...,bN].

И я хочу вычислить их внутреннее произведение: = a1 * coni (b1) + ... + aN * con (bN).(Conj (.) означает комплексную сопряженную операцию)

Я пробовал cublasCdotu, и он просто вычисляет a1 * b1 + ... + aN * bN.

И cublasCdotc возвращает con (a1) * Conv (b1) + ... + ConI (aN) * ConI (BN).

Наконец, я попробовал thrust :: inner_product, и он вычисляет a1 * b1 + ... + aN * bNтоже.

Мой код направления такой: 

typedef thrust::complex<float> comThr;
thrust::host_vector< comThr > x( vec_size );
thrust::generate(x.begin(), x.end(), rand);

thrust::host_vector< comThr > y( vec_size );
thrust::generate(y.begin(), y.end(), rand);
comThr z = thrust::inner_product(x.begin(), x.end(), y.begin(), comThr(0.0f,0.0f) );

Не могли бы вы дать мне несколько советов по этой проблеме?Спасибо!

1 Ответ

3 голосов
/ 26 апреля 2019

Вы можете сделать это с помощью thrust::inner_product. Все, что требуется, это пользовательская двоичная функция, которая реализует a * conj(b), где conj - комплексное сопряжение. Библиотека Thrust включает в себя все необходимые сложные операторы, поэтому реализация проста, как оператор, подобный этому: 

  __host__ __device__
  comThr operator()(comThr a, comThr b) 
  { 
    return  a * thrust::conj(b); 
  };

Полный рабочий пример: 

#include <iostream>
#include "thrust/host_vector.h"
#include "thrust/functional.h"
#include "thrust/complex.h"
#include "thrust/inner_product.h"
#include "thrust/random.h"

typedef thrust::complex<float> comThr;

struct a_dot_conj_b : public thrust::binary_function<comThr,comThr,comThr>
{
  __host__ __device__
  comThr operator()(comThr a, comThr b) 
  { 
    return  a * thrust::conj(b); 
  };
};

__host__ static __inline__ comThr rand_comThr()
{
    return comThr((float)rand()/RAND_MAX, (float)rand()/RAND_MAX);
}

int main()
{
    const int vec_size = 16;

    thrust::host_vector< comThr > x( vec_size );
    thrust::generate(x.begin(), x.end(), rand_comThr);

    thrust::host_vector< comThr > y( vec_size );
    thrust::generate(y.begin(), y.end(), rand_comThr);

    comThr z = thrust::inner_product(x.begin(), x.end(), y.begin(), comThr(0.0f,0.0f),
                                     thrust::plus<comThr>(),  a_dot_conj_b());

    comThr zref(0.0,0.0);
    for(int i=0; i<vec_size; i++) {
        comThr val = x[i] * thrust::conj(y[i]);
        std::cout << i << " " << x[i] << " op " << y[i] << " = " << val  << std::endl;
        zref += val;
    }

    std::cout << "z = " << z << " zref = " << zref << std::endl;

    return 0;
}

, который скомпилируется и будет работать так: 

$ nvcc -arch=sm_52 -o dotprod_thrust dotprod_thrust.cu 

$ ./dotprod_thrust 

0 (0.394383,0.840188) op (0.296032,0.61264) = (0.631482,0.00710744)
1 (0.79844,0.783099) op (0.524287,0.637552) = (0.917879,-0.0984784)
2 (0.197551,0.911647) op (0.972775,0.493583) = (0.642147,0.78932)
3 (0.76823,0.335223) op (0.771358,0.292517) = (0.690638,0.0338566)
4 (0.55397,0.277775) op (0.769914,0.526745) = (0.572826,-0.0779383)
5 (0.628871,0.477397) op (0.891529,0.400229) = (0.751725,0.173921)
6 (0.513401,0.364784) op (0.352458,0.283315) = (0.284301,-0.0168827)
7 (0.916195,0.95223) op (0.919026,0.807725) = (1.61115,0.135091)
8 (0.717297,0.635712) op (0.949327,0.0697553) = (0.725294,0.553463)
9 (0.606969,0.141603) op (0.0860558,0.525995) = (0.126716,-0.307077)
10 (0.242887,0.0163006) op (0.663227,0.192214) = (0.164222,-0.0358752)
11 (0.804177,0.137232) op (0.348893,0.890233) = (0.40274,-0.668025)
12 (0.400944,0.156679) op (0.020023,0.0641713) = (0.0180824,-0.0225919)
13 (0.108809,0.12979) op (0.0630958,0.457702) = (0.0662707,-0.0416127)
14 (0.218257,0.998924) op (0.970634,0.23828) = (0.449871,0.917584)
15 (0.839112,0.512932) op (0.85092,0.902208) = (1.17679,-0.32059)
z = (9.23213,1.02127) zref = (9.23213,1.02127)
...