Pertanyaan vektorisasi eigen dengan array


Saya sedang memproses data cloud point (150k poin per cloud). Saya ingin, untuk setiap titik (x, y), untuk menghitung jarak ke titik referensi O, dan azimut:

for each point p in points
    dx = p.x - ox
    dy = p.y - oy
    d = hypot(dx, dy)
    az = atan2(dy, dx)

Saya memiliki implementasi SSE manual. Saya berharap untuk membuat kode lebih jelas menggunakan eigen:

ArrayXf x(points.size()), y(points.size());
for(unsigned i=0; i<points.size(); ++i) {
    x[i] = points[i].x;
    y[i] = points[i].y;
}
const ArrayXf d = (dx.square() + dy.square()).sqrt();
// implement a polynomial approximation to atan (same as the SSE)

Namun, dari percobaan waktu saya, ini tidak tampak sama sekali karena waktu sama dengan pelaksanaan baseline. Dan saya tahu bahwa SSE2 diaktifkan karena saya menyusun beberapa kode SSE2 dalam file yang sama.

Namun, menurut doc, Eigen memang memanfaatkan SSE2 ketika didukung (dan AVX di 3.3). Apakah hanya untuk operasi vektor dan matrik?

EDIT: Saya mempelajari kode assembly yang diproduksi dan itu berisi beberapa instruksi SSE. Tapi ini masih lambat

EDIT: di sini lebih banyak informasi waktu. Saya mengulang lebih dari 100 frame, sekitar 150k poin per frame.

  • implementasi naif tanpa atan2: 150ms
  • implementasi sse (memproses poin 4 oleh 4 dan membuang beberapa terakhir yang tidak mengisi paket lengkap): 30ms
  • implementasi eigen menggunakan peta eigen: 90ms (diff: 36ms, hypot: 16ms, index: 17ms)

di sini adalah kode eigen saya:

const Eigen::Map<const Eigen::ArrayXf, Eigen::Unaligned, Eigen::InnerStride<4> > px(&(points[0].x), points.size());
const Eigen::Map<const Eigen::ArrayXf, Eigen::Unaligned, Eigen::InnerStride<4> > py(&(points[0].y), points.size());

// difference with the origin (ox and oy are floats)
const Eigen::ArrayXf dx = px - ox, dy = py - oy;

// distance and index
const Eigen::ArrayXf d = sqrt(dx.square() + dy.square());

static const float r_res_mult = 1.0f / r_res; //2x faster than div
const Eigen::ArrayXi didx = (d * r_res_mult).cast<int>();

4
2017-07-07 00:56


asal


Jawaban:


Masalah utama Anda adalah data Anda tidak dalam format yang ramah untuk SIMD. Anda menggunakan larik struct (xyxyxyxyxyxy ...) dan kemudian meng-vectorize kode yang Anda lakukan

for(unsigned i=0; i<points.size(); ++i) {
    x[i] = points[i].x;
    y[i] = points[i].y;
}

yang mengkonversi ke struct of array (xxxxxxxx .... yyyyyyy ...). Konversi ini mahal.

Solusi yang lebih baik adalah sudah memiliki poin Anda disimpan sebagai struct of array. Solusi yang lebih baik adalah menggunakan struct hybrid array alias array struct of array. Untuk SSE, dengan asumsi Anda menggunakan floating point tunggal, Anda kemudian akan melakukan xxxxyyyyxxxxyyyy ....

Selanjutnya saya sarankan Anda menggunakan perpustakaan matematika SIMD. Intel menawarkan SVML yang mahal dan sumber tertutup. AMD menawarkan libm yang gratis tetapi sumber tertutup. Tetapi perpustakaan ini tidak bermain dengan baik di perangkat keras pesaing mereka. Perpustakaan SIMD terbaik adalah Agner Fog Perpustakaan Kelas Vektor (VCL) . Ini open source, gratis, dan dioptimalkan untuk bekerja pada prosesor Intel dan AMD. Ini juga, seperti Eigen, hanya file header dan oleh karena itu, seperti Eigen, Anda tidak perlu mengkompilasi dan menghubungkan pustaka. Anda baru saja menyertakan file header. Di sini adalah bagaimana Anda akan melakukannya untuk SSE atau AVX untuk float (VLC akan meniru AVX pada sistem tanpa AVX).

//    g++ -O3 -Ivectorclass -msse4.2 foo.cpp
// or g++ -O3 -Ivectorclass -mavx foo.cpp
#include <vectorclass.h>
#include <vectormath_trig.h>

struct Point2DBlock {
    float x[8];
    float y[8];
};

int main(void) {
    const int nblocks = 10; //each block contains eight points
    Point2DBlock aosoa[nblocks]; //xxxxxxxxyyyyyyyy xxxxxxxxyyyyyyyy ...
    float ox = 0.0f, oy = 0.0f;
    Vec8f vox = ox, voy = oy;
    for(int i=0; i<nblocks; i++) {
        Vec8f dx = Vec8f().load(aosoa[i].x) - vox;
        Vec8f dy = Vec8f().load(aosoa[i].y) - voy;
        Vec8f d  = sqrt(dx*dx + dy*dy);
        Vec8f az = atan2(dy,dx);
    } 
}

Jika Anda benar-benar membutuhkannya hypot. Anda dapat membangun satu dari VCL menggunakan pseudo-code dari wikipedia.

static inline Vec8f hypot(Vec8f const &x, Vec8f const &y) {
    Vec8f t;
    Vec8f ax = abs(x), ay = abs(y);
    t  = min(ax,ay);
    ax = max(ax,ay);
    t  = t/ax;
    return ax*sqrt(1+t*t);
}

Edit:

Berikut ini adalah metode menggunakan array struct. Ini membutuhkan beberapa menyeret tetapi ini dapat diabaikan dibandingkan dengan perhitungan lainnya. VLC menggunakan pemrograman meta template untuk menentukan metode yang efisien untuk pengacakan.

#include <vectorclass.h>
#include <vectormath_trig.h>

int main(void) {
    const int npoints=80;
    float points[2*npoints]; //xyxyxyxyxyxy...
    float ox = 0.0, oy = 0.0;
    Vec8f vox = ox, voy = oy;
    for(int i=0; i<npoints; i+=16) {
        Vec8f l1 = Vec8f().load(&points[i+0]);
        Vec8f l2 = Vec8f().load(&points[i+8]);
        Vec8f dx = blend8f<0, 2, 4, 6, 8, 10, 12, 14>(l1,l2) - vox;
        Vec8f dy = blend8f<1, 3, 5, 7, 9, 11, 13, 15>(l1,l2) - voy;
        Vec8f d  = sqrt(dx*dx + dy*dy);
        Vec8f az = atan2(dy,dx);
    } 
}

4
2017-07-07 07:39



Salinannya membutuhkan banyak waktu. Sama atau lebih lama dari perhitungan itu sendiri. Anda tidak perlu menyalin data seperti itu. Ini verbose dan mungkin lebih lambat. Anda dapat menggunakan peta sebagai gantinya, atau bahkan menggunakan peta secara langsung untuk menghitung. Saya menulis demo cepat:

int sz = 15000000;
std::vector<Point> points(sz);

Eigen::Map<ArrayXd, Unaligned, InnerStride<2>> mapX(&(points[0].x), sz);
Eigen::Map<ArrayXd, Unaligned, InnerStride<2>> mapY(&(points[0].y), sz);

mapX = ArrayXd::Random(sz);
mapY = ArrayXd::Random(sz);

auto cpstart = std::chrono::high_resolution_clock::now();
ArrayXd x = mapX;
ArrayXd y = mapY;
ArrayXd d;
auto cpend = std::chrono::high_resolution_clock::now();

auto mpSumstart = std::chrono::high_resolution_clock::now();

d = (mapX.square() + mapY.square()).sqrt().eval();

auto mpSumend = std::chrono::high_resolution_clock::now();

std::cout << d.mean() << "\n";

auto arStart = std::chrono::high_resolution_clock::now();

d = (x.square() + y.square()).sqrt().eval();

auto arEnd = std::chrono::high_resolution_clock::now();

std::cout << d.mean() << "\n";

auto elapsed = cpend - cpstart;
std::cout << "Copy: " <<  elapsed.count() << '\n';
std::cout << "Map: " <<  (mpSumend - mpSumstart).count() << '\n';
std::cout << "Array: " <<  (arEnd - arStart).count() << '\n';

Waktu yang saya dapatkan adalah 100 kali panjang larik Anda, saya terlalu malas untuk menulis satu lingkaran untuk menguji lebih baik. Salinan ini membutuhkan sekitar 90ms pada sistem saya (VS2012 / Ox in Release (-DNDEBUG)), versi 185ms yang dipetakan, dan array yang disalin juga sekitar 90ms. Faktor perkiraan dua masuk akal untuk operasi SIMD sebagai versi yang dipetakan melompat setiap ganda lainnya. Jika Anda memiliki struct of array daripada array struct, maka kinerja peta harus sebanding dengan array yang disalin.

EDIT: Saya mendefinisikan EIGEN_DONT_VECTORIZE dan array yang disalin (hampir) menggandakan waktunya (seperti yang diharapkan). Namun, peta itu tetap sama. Ingin tahu. Mungkin harus dilakukan dengan peta yang tidak sejajar. Atau hanya fakta bahwa hanya ada ruang untuk dua ganda dan satu lagi milik peta yang salah.

EDIT 2 Pemikiran bodoh menghantam saya mengenai masalah khusus yang diajukan dalam pertanyaan itu. Anda bisa mengobati x, y nilai sebagai std::complex<double> dan kemudian dimuat sebagai satu blok tanpa salinan memori:

Eigen::Map<ArrayXcd> mapC((std::complex<double>*)(&(points[0].x)), sz);
//...
cd = mapC.cwiseAbs2().sqrt().eval();

Waktunya hanya sedikit lebih lama daripada susunan yang sudah dikopi sebelumnya di komputer saya. Anda juga dapat mengurangi asal sebagai jumlah kompleks yang dilakukan

cd = (mapC - std::complex<double>(ox, oy)).cwiseAbs2().sqrt().eval();

3
2017-07-07 07:02