هل من الممكن أن vectorize myNum += a[b[i]] * c[i];على x86_64?

Question

ما إينترينسيكس وأود أن استخدام vectorize التالية(إذا كان ممكنا أن vectorize) على x86_64?

double myNum = 0;
for(int i=0;i<n;i++){
    myNum += a[b[i]] * c[i]; //b[i] = int, a[b[i]] = double, c[i] = double
}

Answer 1

ها هي ذهابي إليها ، محسّنة واختبارها بالكامل:

#include <emmintrin.h>

__m128d sum = _mm_setzero_pd();
for(int i=0; i<n; i+=2) {
    sum = _mm_add_pd(sum, _mm_mul_pd(
        _mm_loadu_pd(c + i),
        _mm_setr_pd(a[b[i]], a[b[i+1]])
    ));
}

if(n & 1) {
    sum = _mm_add_pd(sum, _mm_set_sd(a[b[n-1]] * c[n-1]));
}

double finalSum = _mm_cvtsd_f64(_mm_add_pd(
    sum, _mm_shuffle_pd(sum, sum, _MM_SHUFFLE2(0, 1))
));

هذا ينتج رمز تجميع جميل جدا باستخدام gcc -O2 -msse2 (4.4.1).

كما يمكنك أن تقول ، وجود حتى n سوف تجعل هذه الحلقة أسرع وكذلك محاذاة c. إذا كنت تستطيع المحاذاة c, ، يتغيرون _mm_loadu_pd ل _mm_load_pd لأوقات تنفيذ أسرع.

Answer 2

أود أن أبدأ بإلغاء حلقة. شيء مثل

double myNum1 = 0, myNum2=0;
for(int i=0;i<n;i+=2)
{
    myNum1 += a[b[ i ]] * c[ i ];
    myNum2 += a[b[i+1]] * c[i+1];
}
// ...extra code to handle the remainder when n isn't a multiple of 2...
double myNum = myNum1 + myNum2;

نأمل أن يتيح أن يسمح للمترجم بالتواصل مع الأحمال مع الحساب ؛ الملف الشخصي وانظر إلى التجمع لمعرفة ما إذا كان هناك تحسن. من الناحية المثالية ، سيقوم برنامج التحويل البرمجي بإنشاء تعليمات SSE ، لكنني لست إذا حدث ذلك في الممارسة العملية.

قد يتيح لك عدم التوصل إلى القيام بذلك:

__m128d sum0, sum1;
// ...initialize to zero...
for(int i=0;i<n;i+=4)
{
    double temp0 = a[b[ i ]] * c[ i ];
    double temp1 = a[b[i+1]] * c[i+1];
    double temp2 = a[b[i+2]] * c[i+2];
    double temp3 = a[b[i+3]] * c[i+3];
    __m128d pair0 = _mm_set_pd(temp0, temp1);
    __m128d pair1 = _mm_set_pd(temp2, temp3);
    sum0 = _mm_add_pd(sum0, pair0);
    sum1 = _mm_add_pd(sum1, pair1);
}
// ...extra code to handle the remainder when n isn't a multiple of 4...
// ...add sum0 and sum1, then add the result's components...

(الاعتذار عن الرمز الكاذب في البداية والنهاية ، أعتقد أن الجزء المهم هو الحلقة). لا أعرف على وجه اليقين ما إذا كان ذلك سيكون أسرع ؛ يعتمد ذلك على اختلافات الكامنات ومدى جودة المترجمات إعادة ترتيب كل شيء. تأكد من أنك ملف تعريف قبل وبعد معرفة ما إذا كان هناك تحسن فعلي.

امل ان يساعد.

Answer 3

معالجات إنتل يمكن إصدار اثنين من عمليات النقطة العائمة ولكن حمولة واحدة في كل دورة ، حتى يصل الذاكرة هي أشد القيد.مع هذا في الاعتبار, تهدف أنا أول من استخدم معبأة الأحمال للحد من عدد من الحمل والتعليمات المستخدمة معبأة الحساب فقط لأنها كانت مريحة.لقد أدركت منذ تشبع عرض النطاق الترددي الذاكرة قد تكون أكبر مشكلة ، كل العبث مع تعليمات SSE قد يكون من السابق لأوانه الأمثل إذا كانت نقطة لجعل رمز تذهب بسرعة بدلا من تعلم vectorize.

SSE

أقل عدد ممكن الأحمال مع أي افتراض على المؤشرات في b يتطلب الفتح الحلقة أربع مرات.واحد 128 بت تحميل يحصل على أربعة مؤشرات من b, اثنين 128 بت الأحمال الحصول على كل زوج المجاورة الزوجي من c, و جمع a مطلوب المستقلة 64 بت الأحمال.هذه الكلمة من 7 دورات في أربعة تكرارات بالنسبة رمز المسلسل.(يكفي أن تشبع ذاكرتي عرض النطاق الترددي إذا كان الوصول إلى a لا ذاكرة التخزين المؤقت بشكل جيد).لقد تركت بعض الأشياء المزعجة مثل التعامل مع عدد من التكرارات التي ليست متعددة من 4.

entry: ; (rdi,rsi,rdx,rcx) are (n,a,b,c)
  xorpd xmm0, xmm0
  xor r8, r8
loop:
  movdqa xmm1, [rdx+4*r8]
  movapd xmm2, [rcx+8*r8]
  movapd xmm3, [rcx+8*r8+8]
  movd   r9,   xmm1
  movq   r10,  xmm1
  movsd  xmm4, [rsi+8*r9]
  shr    r10,  32
  movhpd xmm4, [rsi+8*r10]
  punpckhqdq xmm1, xmm1
  movd   r9,   xmm1
  movq   r10,  xmm1
  movsd  xmm5, [rsi+8*r9]
  shr    r10,  32
  movhpd xmm5, [rsi+8*r10]
  add    r8,   4
  cmp    r8,   rdi
  mulpd  xmm2, xmm4
  mulpd  xmm3, xmm5
  addpd  xmm0, xmm2
  addpd  xmm0, xmm3
  jl loop

الحصول على المؤشرات هو الجزء الأكثر تعقيدا. movdqa الأحمال 128 بت من البيانات عدد صحيح من 16 بايت محاذاة عنوان (Nehalem وقد الكمون العقوبات على خلط "عدد صحيح" و "تعويم" تعليمات SSE). punpckhqdq يتحرك عالية 64 بت منخفضة 64 بت ولكن في صحيح الوضع على عكس يدعى ببساطة movhlpd.32 بت التحولات تتم في سجلات للأغراض العامة. movhpd الأحمال مزدوجة واحدة في الجزء العلوي من xmm التسجيل دون إزعاج الجزء السفلي - يستخدم هذا إلى تحميل عناصر a مباشرة إلى معبأة السجلات.

وهذا رمز واضح أسرع من رمز أعلاه ، والذي هو بدوره أسرع من كود بسيط, و على كل نمط الوصول ولكن الحالة بسيطة B[i] = i حيث السذاجة حلقة هو في الواقع أسرع.كما أنني حاولت بعض الشيء مثل وظيفة في جميع أنحاء SUM(A(B(:)),C(:)) في Fortran التي انتهت أساسا ما يعادل حلقة بسيطة.

أنا اختبرت على Q6600 (65 نانومتر Core 2 في 2.4 Ghz) 4GB DDR2-667 الذاكرة في 4 وحدات.اختبار الذاكرة وعرض النطاق الترددي يعطي عن 5333 MB/s لذلك يبدو أنا أرى سوى قناة واحدة.أنا أجمع مع دبيان في دول مجلس التعاون 4.3.2-1.1, -O3 -Ffast-الرياضيات -msse2 -Ftree-vectorize -std=gnu99.

لاختبار أنا السماح n مليون, تهيئة المصفوفات حتى a[b[i]] و c[i] سواء متساوية 1.0/(i+1), مع عدد قليل من أنماط مختلفة من المؤشرات.واحد يخصص a مع مليون من عناصر مجموعات b عشوائي التقليب ، تخصص آخر a مع 10M عناصر يستخدم كل 10, وآخر يخصص a مع 10M عناصر مجموعات b[i+1] عن طريق إضافة رقم عشوائي من 1 إلى 9 b[i].أنا توقيت متى مكالمة واحدة يأخذ مع gettimeofday, مسح مخابئ عن طريق الاتصال clflush على المصفوفات ، وقياس 1000 التجارب من كل وظيفة.لقد تآمر ممهدة وقت التوزيعات باستخدام بعض التعليمات البرمجية من الشجاعة المعيار (على وجه الخصوص ، نواة الكثافة مقدر في statistics حزمة).

عرض النطاق الترددي

الآن الفعلية ملاحظة هامة حول عرض النطاق الترددي.5333MB/s مع 2.4 Ghz على مدار الساعة هي فقط أكثر من اثنين بايت لكل دورة.بياناتي طويلة بما فيه الكفاية أن لا شيء يجب أن يكون cacheable و ضرب التشغيل من حلقة (16+2*16+4*64) بايت تحميل في التكرار إذا كان كل شيء يفتقد يعطيني بالضبط تقريبا ~5333MB/s من عرض النطاق الترددي النظام لديه.وينبغي أن يكون من السهل جدا أن تشبع هذا النطاق الترددي دون SSE.حتى مع افتراض a تماما مؤقتا فقط القراءة b و c واحد التكرار التحركات 12 بايت من البيانات ، ومن السذاجة أن تبدأ من جديد التكرار من أي وقت مضى الدورة الثالثة مع pipelining.

على افتراض أي شيء أقل من كاملة التخزين المؤقت على a يجعل الحساب و عدد التعليمات حتى أقل من عنق الزجاجة.لن أكون مندهشا إذا كان معظم تسريع في قانون بلدي يأتي من إصدار أقل الأحمال b و c حتى أكثر مساحة مجانية لتتبع التكهن الماضي ذاكرة التخزين المؤقت يخطئ على a.

أوسع الأجهزة قد جعل أكثر الفرق.أ Nehalem نظام تشغيل ثلاث قنوات DDR3-1333 سوف تحتاج إلى نقل 3*10667/2.66 = 12.6 بايت لكل دورة تشبع عرض النطاق الترددي الذاكرة.التي سيكون من المستحيل بالنسبة خيط واحد إذا a يناسب في ذاكرة التخزين المؤقت ولكن في 64 بايت خط ذاكرة التخزين المؤقت يخطئ على ناقلات تضيف ما يصل بسرعة - واحد فقط من أربعة الأحمال في حلقة مفقودة في مخابئ يجلب ما يصل متوسط النطاق الترددي المطلوب إلى 16 بايت/دورة.

Answer 4

إجابة قصيرة لا. إجابة طويلة نعم ، ولكن ليس بكفاءة. سوف تتحمل عقوبة القيام بأحمال غير محددة والتي ستنفي أي نوع من الفوائد. ما لم تتمكن

إذا كنت تعرف مسبقًا ما هي المؤشرات ، فإن أفضل ما لديك هو إلغاء تحديد مؤشرات واضحة وتحديدها. فعلت شيئًا مشابهًا باستخدام تخصص القالب وتوليد الكود. إذا كنت مهتمًا ، يمكنني المشاركة

للإجابة على تعليقك ، يجب عليك التركيز بشكل أساسي على صفيف. أسهل شيء في المحاولة على الفور هو حظر الحلقة بعامل اثنين ، وتحميل منخفضة وعالية على حدة ، ثم استخدام مم*_pd مثل عادة. كود مزيف:

__m128d a, result;
for(i = 0; i < n; i +=2) {
  ((double*)(&a))[0] = A[B[i]];
  ((double*)(&a))[1] = A[B[i+1]];
  // you may also load B using packed integer instruction
  result = _mm_add_pd(result, _mm_mul_pd(a, (__m128d)(C[i])));
}

لا أتذكر أسماء الوظائف بالضبط ، قد أرغب في التحقق. أيضا ، استخدم قصر الكلمة الرئيسية مع المؤشرات إذا كنت تعلم أنه لا يمكن أن تكون هناك مشكلات مستعملة. هذا سيسمح للمترجم أن يكون أكثر عدوانية.

Answer 5

هذا لن يتجهت كما هو ، بسبب عدم التوجيه المزدوج لمؤشرات الصفيف. نظرًا لأنك تعمل مع الزوجي ، لا يوجد شيء أو لا شيء يتم الحصول عليه من SSE ، لا سيما لأن معظم وحدات المعالجة المركزية الحديثة لديها 2 fpus على أي حال.