Как наиболее эффективно применить функтор к подмножеству массива устройств? - PullRequest
Новая
       71

Как наиболее эффективно применить функтор к подмножеству массива устройств?

0 голосов
/ 13 февраля 2019

Я переписываю библиотеку, которая выполняет вычисления и другие операции с данными, хранящимися в непрерывных порциях памяти, чтобы она могла работать на графических процессорах с использованием инфраструктуры CUDA.Данные представляют информацию, которая живет в 4-мерной сетке.Общий размер сетки может варьироваться от 1000 до миллионов точек сетки.Вдоль каждого направления сетка может иметь от 8 до 100 точек.Мой вопрос о том, каков наилучший способ реализации операций над подмножеством сетки.Например, предположим, что моя сетка имеет вид [0, nx) x [0, ny) x [0, nz) x [0, nq), и я хочу реализовать преобразование, которое умножает все точки, индексы которых принадлежат [1, nx-1) x [1, ny-1) x [1, nz-1) x [0, nq-1) на минус 1.

Сейчас я делаю через вложенные циклы.Это скелет кода 

{ 
int nx,ny,nz,nq;
nx=10,ny=10,nz=10,nq=10;
typedef thrust::device_vector<double> Array;
Array A(nx*ny*nz*nq);
thrust::fill(A.begin(), A.end(), (double) 1);

for (auto q=1; q<nq-1; ++q){
for (auto k=1; k<nz-1; ++k){
for (auto j=1; j<ny-1; ++j){
int offset1=+1+j*nx+k*nx*ny+q*nx*ny*nz;
int offset2=offset1+nx-2;
thrust::transform(A.begin()+offset1, 
                  A.begin()+offset2, 
                  thrust::negate<double>());
      }
    }
  }
}

Однако мне интересно, является ли это наиболее эффективным способом, потому что мне кажется, что в этом случае одновременно может быть запущено не более nx-2 потоков.Поэтому я подумал, что, возможно, лучшим способом было бы создать итератор последовательности (возвращающий линейную позицию вдоль массива), сжать его в массив с помощью итератора zip и определить функтор, который проверяет второй элемент кортежа (значение позиции), и если это значение попадает в допустимый диапазон, измените первый элемент кортежа.Тем не менее, может быть лучший способ сделать это.Я новичок в CUDA, и, что еще хуже, я действительно порезал себе зубы с Фортраном, поэтому мне трудно думать вне рамок для цикла ...

1 Ответ

0 голосов
/ 13 февраля 2019

Я не уверен, что является наиболее эффективным способом.Я могу предложить то, что, на мой взгляд, будет более эффективным, чем ваш скелетный код.

Ваше предложение в тексте направлено в правильном направлении.Вместо того, чтобы использовать набор вложенных циклов for, которые будут повторяться потенциально несколько раз, мы должны стремиться выполнить все одним вызовом.Но нам все еще нужно, чтобы один вызов Thrust изменял только значения массива в индексах внутри «кубического» тома, с которым нужно работать.

Мы не хотим использовать метод, включающий тестирование сгенерированного индексапротив действительного объема индекса, однако, как вы, кажется, предлагаете.Это потребует от нас запуска сетки размером с наш массив, даже если мы хотим изменить только ее небольшой объем.

Вместо этого мы запускаем операцию, достаточно большую, чтобы покрыть необходимое количествоэлементы для изменения, и мы создаем функтор, который выполняет линейный индекс -> 4D индекс -> скорректированное преобразование линейного индекса.Затем этот функтор работает в итераторе преобразования, чтобы преобразовать обычную линейную последовательность, начинающуюся с 0, 1, 2 и т. Д., В последовательность, которая начинается и остается в объеме, который нужно изменить.Итератор перестановки затем используется с этой измененной последовательностью для выбора значений массива для изменения.

Вот пример, показывающий разницу во времени для вашего метода вложенного цикла (1) по сравнению с моим (2) длямассив размером 64x64x64x64 и модифицированный том 62x62x62x62: 

$ cat t39.cu
#include <thrust/device_vector.h>
#include <thrust/transform.h>
#include <thrust/iterator/permutation_iterator.h>
#include <thrust/iterator/zip_iterator.h>
#include <thrust/iterator/counting_iterator.h>
#include <thrust/functional.h>
#include <thrust/equal.h>
#include <cassert>
#include <iostream>

struct my_idx
{
  int nx, ny, nz, nq, lx, ly, lz, lq, dx, dy, dz, dq;
  my_idx(int _nx, int _ny, int _nz, int _nq, int _lx, int _ly, int _lz, int _lq, int _hx, int _hy, int _hz, int _hq) {
    nx = _nx;
    ny = _ny;
    nz = _nz;
    nq = _nq;
    lx = _lx;
    ly = _ly;
    lz = _lz;
    lq = _lq;
    dx = _hx - lx;
    dy = _hy - ly;
    dz = _hz - lz;
    dq = _hq - lq;
    // could do a lot of assert checking here
  }

  __host__ __device__
  int operator()(int idx){
    int rx = idx / dx;
    int ix = idx - (rx * dx);
    int ry = rx / dy;
    int iy = rx - (ry * dy);
    int rz = ry / dz;
    int iz = ry - (rz * dz);
    int rq = rz / dq;
    int iq = rz - (rq * dq);
    return (((iq+lq)*nz+iz+lz)*ny+iy+ly)*nx+ix+lx;
  }
};

#include <time.h>
#include <sys/time.h>
#define USECPSEC 1000000ULL

unsigned long long dtime_usec(unsigned long long start){

  timeval tv;
  gettimeofday(&tv, 0);
  return ((tv.tv_sec*USECPSEC)+tv.tv_usec)-start;
}


int main()
{
  int nx,ny,nz,nq,lx,ly,lz,lq,hx,hy,hz,hq;
  nx=64,ny=64,nz=64,nq=64;
  lx=1,ly=1,lz=1,lq=1;
  hx=nx-1,hy=ny-1,hz=nz-1,hq=nq-1;
  thrust::device_vector<double> A(nx*ny*nz*nq);
  thrust::device_vector<double> B(nx*ny*nz*nq);
  thrust::fill(A.begin(), A.end(), (double) 1);
  thrust::fill(B.begin(), B.end(), (double) 1);
  // method 1
  unsigned long long m1_time = dtime_usec(0);
  for (auto q=lq; q<hq; ++q){
    for (auto k=lz; k<hz; ++k){
      for (auto j=ly; j<hy; ++j){
        int offset1=lx+j*nx+k*nx*ny+q*nx*ny*nz;
        int offset2=offset1+(hx-lx);
        thrust::transform(A.begin()+offset1,
                  A.begin()+offset2, A.begin()+offset1,
                  thrust::negate<double>());
      }
    }
  }
  cudaDeviceSynchronize();
  m1_time = dtime_usec(m1_time);

  // method 2
  unsigned long long m2_time = dtime_usec(0);
  auto p = thrust::make_permutation_iterator(B.begin(), thrust::make_transform_iterator(thrust::counting_iterator<int>(0), my_idx(nx, ny, nz, nq, lx, ly, lz, lq, hx, hy, hz, hq)));
  thrust::transform(p, p+(hx-lx)*(hy-ly)*(hz-lz)*(hq-lq), p, thrust::negate<double>());
  cudaDeviceSynchronize();
  m2_time = dtime_usec(m2_time);
  if (thrust::equal(A.begin(), A.end(), B.begin()))
    std::cout << "method 1 time: " << m1_time/(float)USECPSEC << "s method 2 time: " << m2_time/(float)USECPSEC << "s" << std::endl;
  else
    std::cout << "mismatch error" << std::endl;
}
$ nvcc -std=c++11 t39.cu -o t39
$ ./t39
method 1 time: 1.6005s method 2 time: 0.013182s
$
...