👨🏿‍🍳 😢 👎🏽 एन निकायों का कार्य या रसोई छोड़ने के बिना आकाशगंगा को कैसे उड़ाया जाए 🥣 👆🏼 🛡️

बहुत पहले नहीं मैंने लियू किक्सिन द्वारा विज्ञान कथा उपन्यास "द थ्री-बॉडी प्रॉब्लम" पढ़ा। इसमें, कुछ एलियंस को एक समस्या थी - वे नहीं जानते थे कि कैसे, उनके लिए पर्याप्त सटीकता के साथ, अपने घर ग्रह के प्रक्षेपवक्र की गणना करने के लिए। हमारे विपरीत, वे एक तीन सितारा प्रणाली में रहते थे, और ग्रह पर "मौसम" अपने रिश्तेदार स्थान पर दृढ़ता से निर्भर करता था - गर्मी की बढ़ती गर्मी से। और मैंने यह जांचने का फैसला किया कि क्या हम ऐसी समस्याओं को हल कर सकते हैं।

घटना का भौतिकी

समस्या को समझने के लिए, घटना के भौतिकी से निपटना आवश्यक है। शास्त्रीय सिद्धांत के ढांचे में, दो निकायों के आकर्षक बल को न्यूटन के नियम द्वारा निर्धारित किया गया है:

ज ी

$\ vec {F} (\ vec {r} _1, \ vec {r} _2) = - जी m_1 m_2 \ frac {\ vec {r} _1 - \ vec {r} _2} {\ _ vec {r} | _1 - \ vec {r} _2 | ^ 3},$

जहाँ

$\ vec {r} _1, \ vec {r} _2$ - अंतरिक्ष में पिंडों की स्थिति,

$m_1, m_2$ - निकायों का द्रव्यमान,

ज ी

$जी$ - गुरुत्वाकर्षण स्थिर।
की व्यवस्था में

ए न

$एन$ उनमें से प्रत्येक पर निकाय बाकी से आकर्षण बल द्वारा प्रभावित होगा, जो समीकरण द्वारा व्यक्त किया गया है:

$\ vec {F} _n = -G \ sum_ {k \ neq n} m_n m_k \ frac {\ vec {r} _n - \ vec {r} _k} {| \ vec [r} _n - \ vec {r } _k | ^ 3} |$

न्यूटन के दूसरे नियम का उपयोग करते हुए , हम प्रत्येक कण के लिए त्वरण लिखते हैं:

$\ vec {a} _n = \ vec {F} _n / m_n = -G \ sum_ {k \ neq n} m_k \ frac {\ vec {r} _n - \ vec {r_ _c} {\ _cc r} _n - \ vec {r} _k | ^ 3}।$

उस त्वरण को याद करते हुए समन्वय का दूसरी बार व्युत्पन्न होना, हम एक दूसरे क्रम का आंशिक अंतर समीकरण प्राप्त करते हैं, जिसे प्रत्येक शरीर के प्रक्षेपवक्र को प्राप्त करने के लिए हल किया जाना चाहिए:

$\ frac {\ आंशिक ^ 2 \ vec {r} _n} {\ आंशिक t ^ 2} = f_n = -G \ sum_ {k \ neq n} m_k \ frac {\ vec {r} _n_ \ vec {r } _k} {| \ vec {r} _n - \ vec {r} _k | ^ 3} |$

यहां यह ध्यान रखना महत्वपूर्ण है कि किसी फ़ंक्शन की गणना करने की जटिलता

$f_n$ के बराबर है

$O (N ^ 2)$ और अंतःक्रियात्मक निकायों की बढ़ती संख्या के साथ काफी बढ़ जाती है।

गणित

अंतर समीकरणों को हल करने के लिए पहली और सरल विधि ईलर विधि है , जिसे फॉर्म के समीकरणों को हल करने के लिए डिज़ाइन किया गया है:

$\ frac {डाई} {dx} = f (x, y)।$

असतत क्षेत्र में संक्रमण होने पर, हम प्राप्त करते हैं:

$y_ {i} = y_ {i-1} + h f (x_ {i-1}, y_ {i-1}), \ quad i = 1,2,3, \ dots, m,$

जहाँ

$ज$ एकीकरण कदम है, और

$एम$ - एकीकरण चरणों की संख्या। इस प्रकार, अगर हमें एक समय में निकायों की स्थिति की गणना करने की आवश्यकता है

$T$ तो हमें करना चाहिए

$m = T / h$ एकीकरण के कदम। यहां पहली समस्या दिखाई देती है - यदि

$T$ बड़े, तो हमें बड़ी संख्या में एकीकरण कदम उठाने की आवश्यकता है।

हमारी समस्या के लिए यूलर विधि को लागू करने के लिए, इसे एक प्रथम-क्रम प्रणाली में कम किया जाना चाहिए। ऐसा करने के लिए, हम एक अतिरिक्त चर - कण वेग का परिचय देते हैं:

$\ frac {\ आंशिक \ vec {v} _n} {\ आंशिक t} = f_n, \\ \ frac {\ आंशिक \ vec {r} _n} {\ आंशिक t} = \ vec {v} _n।$

विभेदक समीकरणों के सिस्टम को हल करने में दूसरी समस्या समाधान की सटीकता और इसके नियंत्रण की है। सटीकता में दो तरीकों से सुधार किया जा सकता है: एकीकरण चरण को कम करके और सटीकता के उच्च क्रम के साथ एक विधि का चयन करके। दोनों तरीकों से कम्प्यूटेशनल जटिलता बढ़ जाती है, लेकिन अलग-अलग तरीकों से। उदाहरण के लिए, आप शास्त्रीय चौथे क्रम रन-कुट्टा विधि का उपयोग कर सकते हैं, इसके लिए चार फ़ंक्शन गणनाओं की आवश्यकता होती है

$f_n$ हर कदम पर, लेकिन सटीकता का क्रम है

$O (h ^ 4)$ (तुलना के लिए, यूलर विधि में सटीकता का क्रम है

$O (h)$ और एक गणना की आवश्यकता है

$f_n$ )। समाधान की सटीकता को भी कई तरीकों से नियंत्रित किया जा सकता है: विश्लेषणात्मक समाधान के साथ तुलना करें, विभिन्न तरीकों का उपयोग करके या विभिन्न चरणों के साथ हल करें और परिणामों की तुलना करें, समाधान का अनुपालन करना आवश्यक है।
इसके अलावा, इन तरीकों में से प्रत्येक में इसकी कमियां हैं। विश्लेषणात्मक फैसले अनुपस्थित हो सकते हैं, या, यहां तक कि ज्यादातर मामलों में, पूरी तरह से अनुपस्थित हैं। उदाहरण के लिए, हमारे कार्य के लिए

$एन$ टेली विश्लेषणात्मक समाधान केवल के लिए है

$एन = 2$ , लेकिन यह भी तरीकों की सटीकता का परीक्षण करने के लिए पर्याप्त है। किसी समस्या को दो तरीकों से या अलग-अलग चरणों में हल करने से अभिकलन समय बढ़ जाता है, लेकिन यह दृष्टिकोण लगभग किसी भी कार्य पर लागू किया जा सकता है। प्रत्येक कार्य की सीमाएँ नहीं हैं, लेकिन हमारे लिए उनके पास है: एकीकरण के प्रत्येक चरण पर हम संरक्षण कानूनों के कार्यान्वयन को नियंत्रित कर सकते हैं । यह दृष्टिकोण गणना की जटिलता को भी बढ़ाता है, लेकिन इसमें से चुनने के लिए बहुत कुछ है, सभी कणों के संवेग या कोणीय संवेग की गणना से क्रम की जटिलता होती है

$O (N)$ , जबकि एक प्रणाली की कुल ऊर्जा की गणना के लिए आदेश जटिलता है

$O (N ^ 2)$

कुल ऊर्जा की गणना पर ध्यान दें

हमारे मामले में, सिस्टम की कुल ऊर्जा में दो भाग होते हैं - गतिज और संभावित ऊर्जा। काइनेटिक ऊर्जा में सभी निकायों की गतिज ऊर्जा का योग होता है। संभावित ऊर्जा की गणना करने के लिए, हमें शेष कणों के गुरुत्वाकर्षण क्षेत्र में प्रत्येक कण की संभावित ऊर्जा को जोड़ना होगा, इसलिए हमें जोड़ने की आवश्यकता है

$O (N ^ 2)$ शर्तों। कठिनाई यह है कि सभी शब्द बहुत अलग क्रम के हैं, और यहां तक कि दोहरे-सटीक गणनाओं के साथ भी इस मूल्य की गणना अलग-अलग चरणों में तुलना के लिए पर्याप्त सटीकता के साथ संभव नहीं है। इस समस्या को दूर करने के लिए , कहन एल्गोरिथ्म के अनुसार सारांश लागू करना आवश्यक था।

अंजीर। 1. एक अण्डाकार प्रक्षेपवक्र का एक उदाहरण।

पृथ्वी के चारों ओर एक अण्डाकार कक्षा में एक उपग्रह के सरल मामले पर विचार करें। जैसे-जैसे उपग्रह पृथ्वी के करीब आता है, उसकी गति बढ़ती जाती है, और पृथ्वी से दूर जाने पर यह घटता जाता है, तदनुसार, कक्षा के हरे भाग में समय में एकीकरण के चरण को कम करने की संभावना पैदा होती है, और समाधान की सटीकता को बदले बिना इसे लाल और नीले रंग में बढ़ा दिया जाता है। आइए अधिक विस्तार से तुलना करने का प्रयास करें।

तालिका 1. अंतर समीकरणों को हल करने के लिए खोजी गई विधियां

पद	क्रम	विवरण
एडम्स	1-5	एडम्स-बैशफोर्ट विधि
यूलर	1	यूलर विधि
rk4	4	क्लासिक रन-कुट्टा विधि
rkck	5	काश-कार्प विधि
rkdp	5	डोरेमॉन-प्रिंस विधि
rkdverk	6	वर्नर विधि 1) पी। 181
RKF	7	फेलबर्ग विधि 1) पी। 180
rkgl	6	निहित गॉस-लीजेंडरी विधि
rklc	4	लब्धप्रतिष्ठित विधि
trapeze	2	ट्रेपेज़ॉइड विधि

अपने कार्य के लिए सर्वोत्तम विधि का चयन करने के लिए, हम कई ज्ञात विधियों की तुलना करेंगे। ऐसा करने के लिए, हम निकायों की दो प्रणालियों की टक्कर का अनुकरण करते हैं

$N = 512$ और सिमुलेशन के अंत में कुल ऊर्जा , गति और उसके क्षण में सापेक्ष परिवर्तन को मापें (अधिकतम सिमुलेशन समय

$T_ {अधिकतम} = 2.5$ )। इस मामले में, हम एकीकरण विधियों के चरण और मापदंडों को अलग-अलग करेंगे और फ़ंक्शन कॉल की संख्या को मापेंगे

$f_n$ , क्रमशः, उन तरीकों से जिन्हें कम संख्या में कॉल के साथ कम नुकसान होगा, हम अधिक स्वीकार्य मानेंगे।


क)	ख)
चित्रा 2. ऊर्जा में सापेक्ष परिवर्तन ए), गति बी), सिस्टम सिमुलेशन के अंत में $N = 512$ फ़ंक्शन गणना की संख्या के आधार पर विभिन्न विधियों द्वारा निकाय $f_n$ दोहरी सटीकता

चित्रा 2 के रेखांकन से यह देखा जा सकता है कि फ़ंक्शन गणना की मात्रा का सबसे अच्छा अनुपात

$f_n$ और पांचवें क्रम के एडम्स और डॉर्मन-प्रिंस विधियों में निकायों की प्रणाली की ऊर्जा में सापेक्ष परिवर्तन। यह भी देखा जाता है कि सभी तरीकों के लिए संगणना की संख्या में वृद्धि

$f_n$ प्रणाली की गति में सापेक्ष परिवर्तन बढ़ जाता है। ऊर्जा में एक सापेक्ष परिवर्तन के लिए, यह भी ध्यान देने योग्य है, लेकिन केवल कुछ तरीकों के लिए जो दहलीज तक पहुंच सकता है

$dE / E_0 <10 ^ {- 12}$ । इस प्रभाव को अब तरीकों की त्रुटि के लिए जिम्मेदार नहीं ठहराया जा सकता है, लेकिन गणना की त्रुटि के लिए, और सटीकता में और वृद्धि केवल फ़्लोटिंग पॉइंट के साथ गणना की सटीकता में वृद्धि के साथ ही संभव है।


क)	ख)
अंजीर। 3. ऊर्जा में सापेक्ष परिवर्तन ए), गति बी), सिस्टम सिमुलेशन के अंत में $N = 512$ फ़ंक्शन गणना की संख्या के आधार पर विभिन्न विधियों द्वारा निकाय $f_n$ चौगुनी सटीकता के साथ (__float128)

आंकड़े 3 ए और 3 बी बताते हैं कि चौगुनी सटीकता के साथ गणना का उपयोग सापेक्ष ऊर्जा के नुकसान को कम कर सकता है

$10 ^ {- 23}$ , लेकिन आपको यह समझने की आवश्यकता है कि गणना समय डबल परिशुद्धता की तुलना में परिमाण के दो आदेशों द्वारा बढ़ाया जाता है। जैसा कि डबल-सटीक गणनाओं के मामले में, सटीकता का सबसे अच्छा अनुपात और फ़ंक्शन गणनाओं की संख्या

$f_n$ पांचवें क्रम के एडम्स और डॉरमन-प्रिंस के पास।

डॉर्मन-प्रिंस और वर्नर विधियां नेस्टेड विधियों के वर्ग से संबंधित हैं और एक साथ सटीकता के उच्च और निम्न क्रम के लिए दो समाधानों की गणना कर सकते हैं (डॉर्मन-प्रिंस विधि के लिए, 5 और 4 के आदेश, और वर्नर विधि के लिए, 6 और 5 का आदेश देते हैं)। यदि ये दो समाधान बहुत अलग हैं, तो हम वर्तमान एकीकरण कदम को छोटे लोगों में तोड़ सकते हैं। यह हमें एकीकरण चरण को गतिशील रूप से बदलने और इसे केवल उन क्षेत्रों में कम करने की अनुमति देता है जहां इसकी आवश्यकता होती है।

हमारे सिस्टम के एक लंबे सिमुलेशन अंतराल पर पांचवे क्रम के डॉर्मन-प्रिंस, वर्नर और एडम्स विधियों की अधिक विस्तार से तुलना करें।

$T_ {अधिकतम} = 300$ )।

क) मॉडलिंग की प्रक्रिया में ऊर्जा, गति और उसके क्षण में सापेक्ष परिवर्तन

ख) फ़ंक्शन की गणना की संख्या

$f_n$ सिमुलेशन अंतराल पर

$\ Delta T = 0.3$

अंजीर। 4. सिस्टम के भौतिक मापदंडों और फ़ंक्शन गणना की संख्या में सापेक्ष परिवर्तन की निर्भरता

$f_n$ डॉर्मन-प्रिंस चर-चरण विधि के साथ समय-समय पर मॉडलिंग करते हैं

$h = 10 ^ {- 5} \ dots 10 ^ {- 9}$

क) मॉडलिंग की प्रक्रिया में ऊर्जा, गति और उसके क्षण में सापेक्ष परिवर्तन

ख) फ़ंक्शन की गणना की संख्या

$f_n$ सिमुलेशन अंतराल पर

$\ Delta T = 0.3$

अंजीर। 5. सिस्टम के भौतिक मापदंडों और फ़ंक्शन गणना की संख्या में सापेक्ष परिवर्तन की निर्भरता

$f_n$ एक चर चरण के साथ वर्नर विधि द्वारा मॉडलिंग बनाम समय

$h = 10 ^ {- 5} \ dots 10 ^ {- 9}$

अंजीर। 6. एक चरण में पांचवें क्रम के एडम्स-बाशफोर्ट विधि द्वारा मॉडलिंग की प्रक्रिया में ऊर्जा, गति और उसके क्षण में सापेक्ष परिवर्तन

$h = 10 ^ {- 6}$

अंजीर 7. फ़ंक्शन गणना की संख्या पर निर्भरता

$f_n$ पाँचवें क्रम के एडम्स, डॉर्मन-प्रिंस और वर्नर सिमुलेशन समय से

यह देखा जा सकता है कि हमारी प्रणाली के विकास के प्रारंभिक चरण में (

$T <25$ ) सभी तीन विधियां समान विशेषताएं दिखाती हैं, लेकिन बाद के चरणों में सिस्टम में कुछ घटनाएं होती हैं, जिसके परिणामस्वरूप सिस्टम के मुख्य मापदंडों में त्रुटियां (कुल ऊर्जा, गति और इसकी गति) तेजी से कूदती हैं। लेकिन डॉर्मन-प्रिंस और वर्नर विधियां इन परिवर्तनों के साथ "जटिल" खंडों में एकीकरण कदम को कम करने की क्षमता के कारण बेहतर तरीके से सामना करती हैं, जिसके परिणामस्वरूप फ़ंक्शन की संख्या बढ़ जाती है

$f_n$ जैसा कि आंकड़े 4 बी और 5 बी में देखा गया है, लेकिन गणना की कुल संख्या

$f_n$ नेस्टेड तरीके एडम्स-बाशफोर्ट विधि से छोटे हैं, जैसा कि चित्र 7 में देखा जा सकता है।

मुझे आश्चर्य है कि इन क्षणों में सिस्टम का क्या हुआ

वीडियो 1. 512 निकायों की एक प्रणाली मॉडलिंग। डोरेमॉन-प्रिंस विधि। गतिशील कदम

$h = 10 ^ {- 5} \ dots 10 ^ {- 9}$

वीडियो समय में उस बिंदु तक प्रदर्शित करता है

$T = 25$ गति अपेक्षाकृत शांत है, और "आकाशगंगाओं" के केंद्रों के टकराने के बाद होता है, जो प्रक्षेपवक्र में तेज बदलाव और कुछ कणों के वेग में एक मजबूत वृद्धि की ओर जाता है। इसके अलावा, समाधान की सटीकता बनाए रखने के लिए, एकीकरण चरण को कम करना आवश्यक है। नेस्टेड विधियाँ स्वचालित रूप से कर सकती हैं; रेखांकन बताते हैं कि सिस्टम के विकास के कुछ हिस्सों में, एकीकरण कदम को परिमाण के लगभग दो आदेशों से कम कर दिया गया था

$10 ^ {- 5}$ को

$h = 10 ^ {- 7}$ । सिस्टम के विकास के पूरे अंतराल पर एडम्स पद्धति और इस तरह के एक कदम का उपयोग करते समय, हम एक समाधान के लिए इंतजार नहीं करेंगे।

परिणाम

समाधान के लिए, नेस्टेड विधियों का उपयोग करना बेहतर है जो आपको एकीकरण चरण को गतिशील रूप से नियंत्रित करने की अनुमति देता है, और इसे केवल प्रक्षेपवक्र के "जटिल" खंडों पर कम करता है।

उच्चतम आदेश विधियों का पीछा न करें। यहां तक कि डेटा प्रकार 'डबल' का उपयोग करते समय, वे अपनी संभावित क्षमताओं तक नहीं पहुंचते हैं, डेटा प्रकारों का अधिक सटीकता के साथ उपयोग करने से समस्या को हल करने में लगने वाला समय बहुत बढ़ जाता है।

सीपीयू कार्यान्वयन

अब जबकि समीकरणों को हल करने के लिए एक विधि का विकल्प परिभाषित किया गया है, आइए हम प्रत्येक कण के लिए बातचीत बल की गणना का पता लगाने की कोशिश करें। हम सभी कणों के लिए एक दोहरा चक्र प्राप्त करते हैं:

कार्यान्वयन कोड 'सरल'

for(size_t body1 = 0; body1 < count; ++body1) { const nbvertex_t v1(rx[ body1 ], ry[ body1 ], rz[ body1 ]); nbvertex_t total_force; for(size_t body2 = 0; body2 != count; ++body2) { if(body1 == body2) { continue; } const nbvertex_t v2(rx[ body2 ], ry[ body2 ], rz[ body2 ]); const nbvertex_t force(m_data->force(v1, v2, mass[body1], mass[body2])); total_force += force; } frx[body1] = vx[body1]; fry[body1] = vy[body1]; frz[body1] = vz[body1]; fvx[body1] = total_force.x / mass[body1]; fvy[body1] = total_force.y / mass[body1]; fvz[body1] = total_force.z / mass[body1]; }

प्रत्येक शरीर के लिए गुरुत्वाकर्षण बलों की स्वतंत्र रूप से गणना की जाती है, और सभी प्रोसेसर कोर का उपयोग करने के लिए, पहले चक्र से पहले ओपनएमपी निर्देश लिखना पर्याप्त है:

'Openmp' के कार्यान्वयन से कोड का एक टुकड़ा

 #pragma omp parallel for for(size_t body1 = 0; body1 < count; ++body1)

क्योंकि चूंकि प्रत्येक शरीर रैम के साथ प्रोसेसर की बातचीत की संख्या को कम करने और कैश के उपयोग में सुधार करने के लिए प्रत्येक के साथ बातचीत करता है, हमारे पास डेटा के भाग को कैश में लोड करने और इसे बार-बार उपयोग करने की क्षमता है:

कार्यान्वयन कोड 'ओपनैंप + ब्लॉक'

 #pragma omp parallel for for(size_t n1 = 0; n1 < count; n1 += BLOCK_SIZE) { nbcoord_t x1[BLOCK_SIZE]; nbcoord_t y1[BLOCK_SIZE]; nbcoord_t z1[BLOCK_SIZE]; nbcoord_t m1[BLOCK_SIZE]; nbvertex_t total_force[BLOCK_SIZE]; for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { size_t local_n1 = b1 + n1; x1[b1] = rx[local_n1]; y1[b1] = ry[local_n1]; z1[b1] = rz[local_n1]; m1[b1] = mass[local_n1]; } for(size_t n2 = 0; n2 < count; n2 += BLOCK_SIZE) { nbcoord_t x2[BLOCK_SIZE]; nbcoord_t y2[BLOCK_SIZE]; nbcoord_t z2[BLOCK_SIZE]; nbcoord_t m2[BLOCK_SIZE]; for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2) { size_t local_n2 = b2 + n2; x2[b2] = rx[local_n2]; y2[b2] = ry[local_n2]; z2[b2] = rz[local_n2]; m2[b2] = mass[n2 + b2]; } for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { const nbvertex_t v1(x1[ b1 ], y1[ b1 ], z1[ b1 ]); for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2) { const nbvertex_t v2(x2[ b2 ], y2[ b2 ], z2[ b2 ]); const nbvertex_t force(m_data->force(v1, v2, m1[b1], m2[b2])); total_force[b1] += force; } } } for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { size_t local_n1 = b1 + n1; frx[local_n1] = vx[local_n1]; fry[local_n1] = vy[local_n1]; frz[local_n1] = vz[local_n1]; fvx[local_n1] = total_force[b1].x / m1[b1]; fvy[local_n1] = total_force[b1].y / m1[b1]; fvz[local_n1] = total_force[b1].z / m1[b1]; } }

आगे के अनुकूलन में मुख्य चक्र में बल की गणना करने और शरीर द्रव्यमान m1 [b1] द्वारा विभाजन और गुणन को समाप्त करने के कार्य की सामग्री को शामिल किया गया है। इस तथ्य के अलावा कि हमने गणनाओं को थोड़ा कम कर दिया है, कंपाइलर इस तरह के एक विस्तारित चक्र पर एसएसई और एवीएक्स प्रोसेसर के वेक्टर निर्देशों को लागू करने में सक्षम होगा।

कार्यान्वयन कोड 'ओपनैंप + ब्लॉक + ऑप्टिमाइज़ेशन'

 #pragma omp parallel for for(size_t n1 = 0; n1 < count; n1 += BLOCK_SIZE) { nbcoord_t x1[BLOCK_SIZE]; nbcoord_t y1[BLOCK_SIZE]; nbcoord_t z1[BLOCK_SIZE]; nbcoord_t total_force_x[BLOCK_SIZE]; nbcoord_t total_force_y[BLOCK_SIZE]; nbcoord_t total_force_z[BLOCK_SIZE]; for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { size_t local_n1 = b1 + n1; x1[b1] = rx[local_n1]; y1[b1] = ry[local_n1]; z1[b1] = rz[local_n1]; total_force_x[b1] = 0; total_force_y[b1] = 0; total_force_z[b1] = 0; } for(size_t n2 = 0; n2 < count; n2 += BLOCK_SIZE) { nbcoord_t x2[BLOCK_SIZE]; nbcoord_t y2[BLOCK_SIZE]; nbcoord_t z2[BLOCK_SIZE]; nbcoord_t m2[BLOCK_SIZE]; for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2) { size_t local_n2 = b2 + n2; x2[b2] = rx[local_n2]; y2[b2] = ry[local_n2]; z2[b2] = rz[local_n2]; m2[b2] = mass[n2 + b2]; } for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { for(size_t b2 = 0; b2 != BLOCK_SIZE; ++b2) { nbcoord_t dx = x1[b1] - x2[b2]; nbcoord_t dy = y1[b1] - y2[b2]; nbcoord_t dz = z1[b1] - z2[b2]; nbcoord_t r2(dx * dx + dy * dy + dz * dz); if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = (m2[b2]) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; total_force_x[b1] -= dx; total_force_y[b1] -= dy; total_force_z[b1] -= dz; } } } for(size_t b1 = 0; b1 != BLOCK_SIZE; ++b1) { size_t local_n1 = b1 + n1; frx[local_n1] = vx[local_n1]; fry[local_n1] = vy[local_n1]; frz[local_n1] = vz[local_n1]; fvx[local_n1] = total_force_x[b1]; fvy[local_n1] = total_force_y[b1]; fvz[local_n1] = total_force_z[b1]; } }

तालिका 2. फ़ंक्शन की गणना समय (सेकंड में) की निर्भरता

$f_n$ बातचीत करने वाले निकायों की संख्या पर

$एन$ विभिन्न CPU कार्यान्वयन के लिए

$एन$	2048	4096	8192	16384	32768
सरल	0.0425	0.1651	0.6594	2.65	10.52
OpenMP	0.0078	0.0260	0.1079	0.417	1.655
Openmp + ब्लॉक + अनुकूलन	0.0037	0.0128	0.0495	0.194	0.774

सिस्टम पैरामीटर:

प्रणाली: डेबियन 9, इंटेल कोर i7-5820K (6 कोर)
संकलक: जीसीसी 6.3.0

यह स्पष्ट रूप से देखा गया है कि ओपनएमपी समर्थन वाले संस्करण को कोर की संख्या के संदर्भ में छह बार त्वरित किया जाता है, और अनुकूलित संस्करण दो बार से थोड़ा अधिक तेजी से होता है। इसलिए, अनुकूलन के साथ, आपको केवल संगामिति पर भरोसा नहीं करना चाहिए। दिलचस्प बात यह है कि एकल-स्ट्रीम गणना में, प्रोसेसर ने 3.6 गीगाहर्ट्ज की आवृत्ति पर काम किया, समानांतर संस्करण (ओपनम्प) में इसने आवृत्ति को 3.4 गीगाहर्ट्ज पर गिरा दिया, और समानांतर और ऑप्टिमाइज़्ड (ओपनएम्प + ब्लॉक + ऑप्टिमाइज़ेशन) में यह घटकर 3.3 गीगाहर्ट्ज़ हो गया, लेकिन यह इसने उसे 13.6 गुना तेजी से काम करने से नहीं रोका। यह भी देखा जाता है कि समस्या के आकार में वृद्धि के साथ गणना समय में वृद्धि द्विघात है, और आगे की वृद्धि

$एन$ उचित समय में कार्य को अस्वीकार्य बनाता है।

GPU कार्यान्वयन

लेकिन गणनाओं को और भी तेज करने की इच्छा है। त्वरण के लिए कई दिशाएँ उपलब्ध हैं: GPU अभिकलन, फ़ंक्शन सन्निकटन

$f_n$ । सबसे पहले, GPU कंप्यूटिंग के लिए, मैंने OpenCL तकनीक को चुना। अधिक सुविधाजनक काम के लिए, CLHPP पुस्तकालय का उपयोग किया गया था। ओपनसीएल का मुख्य लाभ यह है कि कोड प्रोसेसर और जीपीयू पर चलाया जा सकता है, जो लेखन और डिबगिंग को सरल बनाता है, साथ ही चलाने के लिए हार्डवेयर की सूची का विस्तार करता है। Oclgrind टूल डिबगिंग में मदद करता है, जो रनटाइम में गलत OpenCL API कॉल और मेमोरी एक्सेस समस्याओं को दिखाता है।

OpenCL

ओपनसीएल के साथ शुरुआत करने के लिए, आपको उपलब्ध प्लेटफार्मों की एक सूची प्राप्त करने की आवश्यकता है। सबसे आम प्लेटफॉर्म एएमडी, इंटेल और एनवीडिया हैं।

कोड

 std::vector<cl::Platform> platforms; cl::Platform::get(&platforms);

अगला, एक प्लेटफ़ॉर्म चुनने के बाद, आपको उस कंप्यूटिंग डिवाइस का चयन करना होगा जो यह प्लेटफ़ॉर्म दर्शाता है:

कोड

 const cl::Platform& platform(platforms[platform_n]); std::vector<cl::Device> all_devices; platform.getDevices(CL_DEVICE_TYPE_ALL, &all_devices);

और तैयारी के चरण के अंत में, एक संदर्भ और कतार बनाना आवश्यक है, जिसके भीतर स्मृति आवंटित की जाएगी और गणना की जाएगी। उदाहरण के लिए, एक संदर्भ जो चयनित मंच के सभी कंप्यूटिंग उपकरणों को जोड़ता है, इस प्रकार बनाया गया है:

प्रसंग और कतार निर्माण संहिता

 cl::Context context(all_devices); std::vector<cl::CommandQueue> queues; for(cl::Device device: all_devices) queues.push_back(cl::CommandQueue(context, device));

एक कंप्यूटिंग डिवाइस के लिए स्रोत कोड को डाउनलोड करने के लिए, इसे संकलित किया जाना चाहिए, इसके लिए cl :: कार्यक्रम वर्ग का इरादा है।

कर्नेल संकलन कोड

 std::vector< std::string > source_data; cl::Program::Sources sources; for(int i = 0; i != files.size(); ++i) { source_data.push_back(load_program(files[i]));//   sources.push_back(std::make_pair(source_data.back().data(), source_data.back().size())); } cl::Program prog(context, sources); devices.push_back(all_devices); prog.build(devices, options);

किसी फ़ंक्शन (कर्नेल) के मापदंडों का वर्णन करने के लिए जो एक कंप्यूटिंग डिवाइस पर निष्पादित किया जाता है, एक सीएल टेम्पलेट है: make_kernel।

इंटरएक्शन ताकत गणना कर्नेल उदाहरण

 typedef cl::make_kernel< cl_int, cl_int, //Block offset cl::Buffer, //mass cl::Buffer, //y cl::Buffer, //f cl_int, cl_int, //yoff,foff cl_int, cl_int //points_count,stride > ComputeBlock;

इसके अलावा, सब कुछ सरल है: हम कर्नेल के प्रकार के साथ एक चर घोषित करते हैं, संकलित कार्यक्रम और कम्प्यूटेशनल कर्नेल के नाम को पास करते हैं, हम कर्नेल को लगभग एक सामान्य फ़ंक्शन की तरह शुरू कर सकते हैं।

कर्नेल लॉन्च कोड

 ComputeBlock fcompute(prog, "ComputeBlockLocal"); cl::NDRange global_range(device_data_size); cl::NDRange local_range(block_size); cl::EnqueueArgs eargs(ctx.m_queue, global_range, local_range); fcompute(eargs, ...  ); //,   .

ओपनसीएल के लिए कम्प्यूटेशनल कोर सीपीयू संस्करण के विपरीत, केवल सीपीयू संस्करण के लिए 'ओपनएम्प + ब्लॉक + ऑप्टिमाइज़ेशन' विकल्प के समान है, पहला चक्र ओपनसीएल का उपयोग करके नियंत्रित किया जाता है (साइकिल रेंज वैश्विक स्टार्टअप प्रणाली से निर्धारित होती है, जो कि स्टार्टअप स्टार्टअप कोड से होती है, और करंट इट्रिशन नंबर कर्नेल से उपलब्ध होता है। फ़ंक्शन का उपयोग करना get_global_id (0))। सबसे पहले, शरीर के डेटा का हिस्सा स्थानीय मेमोरी से लोड किया जाता है, फिर संसाधित किया जाता है। स्थानीय मेमोरी समूह में सभी थ्रेड्स के लिए आम है, इसलिए समूह के लिए डाउनलोड एक बार होता है, और समूह में प्रत्येक थ्रेड द्वारा संसाधित किया जाता है, क्योंकि स्थानीय मेमोरी वैश्विक की तुलना में बहुत तेज है, गणना बहुत तेज है।

कर्नेल कोड

 __kernel void ComputeBlockLocal(int offset_n1, int offset_n2, __global const nbcoord_t* mass, __global const nbcoord_t* y, __global nbcoord_t* f, int yoff, int foff, int points_count, int stride) { int n1 = get_global_id(0) + offset_n1; __global const nbcoord_t* rx = y + yoff; __global const nbcoord_t* ry = rx + stride; __global const nbcoord_t* rz = rx + 2 * stride; __global const nbcoord_t* vx = rx + 3 * stride; __global const nbcoord_t* vy = rx + 4 * stride; __global const nbcoord_t* vz = rx + 5 * stride; __global nbcoord_t* frx = f + foff; __global nbcoord_t* fry = frx + stride; __global nbcoord_t* frz = frx + 2 * stride; __global nbcoord_t* fvx = frx + 3 * stride; __global nbcoord_t* fvy = frx + 4 * stride; __global nbcoord_t* fvz = frx + 5 * stride; nbcoord_t x1 = rx[n1]; nbcoord_t y1 = ry[n1]; nbcoord_t z1 = rz[n1]; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; __local nbcoord_t x2[NBODY_DATA_BLOCK_SIZE]; __local nbcoord_t y2[NBODY_DATA_BLOCK_SIZE]; __local nbcoord_t z2[NBODY_DATA_BLOCK_SIZE]; __local nbcoord_t m2[NBODY_DATA_BLOCK_SIZE]; // NB! get_local_size(0) == NBODY_DATA_BLOCK_SIZE for(int b2 = 0; b2 < points_count; b2 += NBODY_DATA_BLOCK_SIZE) { int n2 = b2 + offset_n2 + get_local_id(0); // Copy data block to local memory x2[ get_local_id(0) ] = rx[n2]; y2[ get_local_id(0) ] = ry[n2]; z2[ get_local_id(0) ] = rz[n2]; m2[ get_local_id(0) ] = mass[n2]; // Synchronize local work-items copy operations barrier(CLK_LOCAL_MEM_FENCE); nbcoord_t local_res_x = 0.0; nbcoord_t local_res_y = 0.0; nbcoord_t local_res_z = 0.0; for(int local_n2 = 0; local_n2 != NBODY_DATA_BLOCK_SIZE; ++local_n2) { nbcoord_t dx = x1 - x2[local_n2]; nbcoord_t dy = y1 - y2[local_n2]; nbcoord_t dz = z1 - z2[local_n2]; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = (m2[local_n2]) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; local_res_x -= dx; local_res_y -= dy; local_res_z -= dz; } // Synchronize local work-items computations barrier(CLK_LOCAL_MEM_FENCE); res_x += local_res_x; res_y += local_res_y; res_z += local_res_z; } frx[n1] = vx[n1]; fry[n1] = vy[n1]; frz[n1] = vz[n1]; fvx[n1] = res_x; fvy[n1] = res_y; fvz[n1] = res_z; }

CUDA

NVidia CUDA प्लेटफ़ॉर्म के लिए कार्यान्वयन OpenCL की तुलना में थोड़ा सरल है, हमें डिवाइस संदर्भ बनाने और निष्पादन कतार का प्रबंधन करने की आवश्यकता नहीं है (कम से कम तब तक जब तक हम बहु-जीपीयू कार्यान्वयन नहीं करना चाहते)। OpenCL के मामले में, हमें GPU पर मेमोरी आवंटित करने की आवश्यकता है, हमारे डेटा को इसमें कॉपी करें, और फिर हम कंप्यूटिंग कोर शुरू कर सकते हैं।

आप यहां CUDA के साथ काम करने के बारे में अधिक पढ़ सकते हैं ।

CUDA कर्नेल लॉन्च कोड

 dim3 grid(count / block_size); dim3 block(block_size); size_t shared_size(4 * sizeof(nbcoord_t) * block_size); kfcompute <<< grid, block, shared_size >>> (... ...);

OpenCL के विपरीत, CUDA पुनरावृत्तियों की पूरी श्रृंखला को निर्दिष्ट नहीं करता है (OpenCL कार्यान्वयन में यह global_range है), लेकिन ग्रिड में ग्रिड आकार और ब्लॉक आकार सेट करता है, तदनुसार वर्तमान शरीर संख्या की गणना थोड़ा बदल जाती है, अन्यथा कर्नेल अन्य नामों के अपवाद के साथ OpenCL के समान है। साझा मेमोरी के लिए सिंक्रोनाइज़ेशन और स्पेसियर फ़ंक्शंस। CUDA की एक और उपयोगी विशिष्ट विशेषता यह है कि कर्नेल शुरू होने पर हम साझा मेमोरी के आवश्यक आकार को निर्दिष्ट कर सकते हैं। OpenCL कार्यान्वयन के रूप में, प्रत्येक पुनरावृत्ति खंड की शुरुआत में हम डेटा के हिस्से को साझा मेमोरी में कॉपी करते हैं और फिर ब्लॉक के सभी थ्रेड्स से इस मेमोरी के साथ काम करते हैं।

CUDA कर्नेल कोड

 __global__ void kfcompute(int offset_n2, const nbcoord_t* y, int yoff, nbcoord_t* f, int foff, const nbcoord_t* mass, int points_count, int stride) { int n1 = blockDim.x * blockIdx.x + threadIdx.x; const nbcoord_t* rx = y + yoff; const nbcoord_t* ry = rx + stride; const nbcoord_t* rz = rx + 2 * stride; nbcoord_t x1 = rx[n1]; nbcoord_t y1 = ry[n1]; nbcoord_t z1 = rz[n1]; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; extern __shared__ nbcoord_t shared_xyzm_buf[]; nbcoord_t* x2 = shared_xyzm_buf; nbcoord_t* y2 = x2 + blockDim.x; nbcoord_t* z2 = y2 + blockDim.x; nbcoord_t* m2 = z2 + blockDim.x; for(int b2 = 0; b2 < points_count; b2 += blockDim.x) { int n2 = b2 + offset_n2 + threadIdx.x; // Copy data block to local memory x2[ threadIdx.x ] = rx[n2]; y2[ threadIdx.x ] = ry[n2]; z2[ threadIdx.x ] = rz[n2]; m2[ threadIdx.x ] = mass[n2]; // Synchronize local work-items copy operations __syncthreads(); nbcoord_t local_res_x = 0.0; nbcoord_t local_res_y = 0.0; nbcoord_t local_res_z = 0.0; for(int n2 = 0; n2 != blockDim.x; ++n2) { nbcoord_t dx = x1 - x2[n2]; nbcoord_t dy = y1 - y2[n2]; nbcoord_t dz = z1 - z2[n2]; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = (m2[n2]) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; local_res_x -= dx; local_res_y -= dy; local_res_z -= dz; } // Synchronize local work-items computations __syncthreads(); res_x += local_res_x; res_y += local_res_y; res_z += local_res_z; } n1 += foff; f[n1 + 3 * stride] = res_x; f[n1 + 4 * stride] = res_y; f[n1 + 5 * stride] = res_z; }

तालिका 3. गणना समय (सेकंड में) फ़ंक्शन की निर्भरता

$f_n$ बातचीत करने वाले निकायों की संख्या पर

$N$ विभिन्न GPU कार्यान्वयन और बेहतर CPU कार्यान्वयन के लिए

$N$	4096	8192	16384	32768	65536	131,072
ओपन + ब्लॉक + अनुकूलन	0.0128	0.0495	0.194	0.774	---	---
ओपनसीएल + आधा एनवीडिया के 80	0.004	0.008	0.026	0.134	0.322	1.18
क्यूडा + आधा एनवीडिया K80	0.004	0.008	0.0245	0.115	0.291	1.13

एनवीडिया K80 कहां से लाएं

OpenCL CUDA Colab Google, NVidia K80. .

सामान्य तौर पर, परिणाम निराशाजनक है, सीपीयू कार्यान्वयन की तुलना में केवल 5-6 गुना तेजी से। यहां तक कि अगर हम पूरे K80 पर गणना करते हैं, तो भी हमें 12 गुना तक त्वरण मिलेगा, लेकिन तब से चूंकि कार्य की जटिलता द्विघात है, तो एक उचित समय में हम 32768 अंतःक्रियात्मक निकायों को संसाधित नहीं कर पाएंगे, लेकिन 131072, जो केवल 4 गुना अधिक है।

समारोह स्वीकृति $f_n$

यदि आप फ़ंक्शन को करीब से देखते हैं, जो दो निकायों के आकर्षक बल द्वारा निर्धारित किया गया है, तो यह स्पष्ट है कि यह दूरी के साथ द्विघात रूप से घटता है। इसलिए, हम करीबी निकायों के बीच और लगभग दूर के लोगों के बीच बातचीत के बल की सही गणना कर सकते हैं। एक प्रसिद्ध दृष्टिकोण डी। बार्न्स और पी। हैट द्वारा प्रस्तावित ट्रेब्लॉस
एल्गोरिथ्म है । ऑक्टो- ट्री को सिम्युलेटेड स्पेस में बनाया जा रहा हैयुक्त शरीर के निर्देशांक और द्रव्यमानों को इसके पत्तों में समाहित करता है। माता-पिता के नोड्स में द्रव्यमान का केंद्र होता है, बच्चे के नोड्स का कुल द्रव्यमान और बच्चे के शरीर के चारों ओर वर्णित गोले की त्रिज्या। पेड़ की जड़ में सभी निकायों के द्रव्यमान का केंद्र, उनका कुल द्रव्यमान और उनके चारों ओर वर्णित गोले की त्रिज्या होती है। अंतःक्रिया बल की गणना करते समय, पेड़ की जड़ से दूरी को पहले माना जाता है यदि नोड से उसकी त्रिज्या की दूरी का अनुपात कुछ स्थिर से अधिक हो

$\lambda_{crit}$ , फिर रूट को एक निकाय माना जाता है, जिसमें समरूप दिशाओं के केंद्र के समतुल्य समतुल्य पिंड होते हैं, जिसमें सम्‍मिलित पिंडों के द्रव्यमान के योग के बराबर द्रव्यमान होता है, यदि अनुपात के बराबर या उससे कम होता है

$\lambda_{crit}$ , फिर प्रक्रिया प्रत्येक बच्चे के नोड के लिए पुनरावृत्ति दोहराई जाती है। यदि एक पेड़ का एक पत्ता तक पहुंच जाता है, तो बातचीत बल को सामान्य तरीके से माना जाता है। इस प्रकार, यदि एक शरीर को अन्य निकायों के कॉम्पैक्ट समूह से बहुत दूर कर दिया जाता है, तो इस समूह को एक शरीर के रूप में प्रस्तुत किया जाता है, और बातचीत बल की गणना प्रत्येक शरीर तक नहीं, बल्कि केवल एक शरीर तक की जाती है। इसके कारण, एल्गोरिथ्म की जटिलता कम हो जाती है

$O(N^2)$ को

$O(N \log(N))$ सटीकता के कुछ नुकसान की कीमत पर।

अपने दृष्टिकोण में, मैंने ऑक्ट्री ट्री का उपयोग नहीं करने का फैसला किया , लेकिन केडी ट्री, क्योंकि यह उपयोग करने के लिए आसान है और ऑक्ट्री की तुलना में कम स्टोरेज ओवरहेड है।

चलो सीपीयू पर कार्यान्वयन पर वापस जाते हैं। एक केडी-ट्री नोड को एक वर्ग के रूप में दर्शाया जा सकता है जिसमें बाएं और दाएं वंश के लिए निर्देशांक और निर्देशांक और द्रव्यमान के बारे में जानकारी है:

केडी वृक्ष नोड

 class node { node* m_left; //!<   node* m_right; //!<   nbvertex_t m_mass_center; //!<     nbcoord_t m_mass; //!<   nbcoord_t m_radius_sqr; //!<    ,   lambda_crit nbvertex_t m_bmin; //!<     nbvertex_t m_bmax; //!<     size_t m_body_n; //!<  ,    };

पेड़ को संचय करने की इस पद्धति के साथ, हमारे पास पेड़ को हटाने के लिए दो संभावित विकल्प हैं: या तो स्पष्ट पुनरावर्तन का उपयोग करें, या स्वयं स्टैक का उपयोग करें। मैं दूसरे विकल्प पर बस गया।

एक पेड़ पर चढ़कर बातचीत बल की गणना

 nbvertex_t force_compute(const nbvertex_t& v1, const nbcoord_t mass1) { nbvertex_t total_force; node* stack_data[MAX_STACK_SIZE] = {}; node** stack = stack_data; node** stack_head = stack; *stack++ = m_root; while(stack != stack_head) { node* curr = *--stack; const nbcoord_t distance_sqr((v1 - curr->m_mass_center).norm()); if(distance_sqr > curr->m_radius_sqr) {//  ,       total_force += force(v1, curr->m_mass_center, mass1, curr->m_mass); } else {//   ,      if(curr->m_right != NULL) { *stack++ = curr->m_right; } if(curr->m_left != NULL) { *stack++ = curr->m_left; } } } return total_force; }

जैसा कि "सटीक" सीपीयू कार्यान्वयन के मामले में, प्रत्येक शरीर के लिए बल गणना फ़ंक्शन को कहा जाता है। OpenMP निर्देशों का उपयोग करके सभी निकायों में चक्र को आसानी से समानांतर किया जा सकता है।

लेकिन पड़ोसी लूप पुनरावृत्तियों इस मामले में पेड़ के पूरी तरह से अलग हिस्सों को संदर्भित करेगा, जो प्रोसेसर कैश का कुशल उपयोग करने की अनुमति नहीं देता है। इस समस्या को दूर करने के लिए, सभी निकायों को मूल क्रम में नहीं, बल्कि उस क्रम में रखा जाना चाहिए, जिस क्रम में निकाय kd-tree की पत्तियों में स्थित होते हैं, तब पड़ोसी पुनरावृत्तियां उन निकायों के लिए घटित होंगी जो अंतरिक्ष में करीब हैं, और लगभग समान पथों के साथ पेड़ को पीछे छोड़ देंगे।

पेड़ की पत्ती का फंदा

 template<class Visitor> void traverse(Visitor visit) { node* stack_data[MAX_STACK_SIZE] = {}; node** stack = stack_data; node** stack_head = stack; *stack++ = m_root; while(stack != stack_head) { node* curr = *--stack; if(curr->m_radius_sqr > 0) {//  .    . if(curr->m_left != NULL) { *stack++ = curr->m_left; } if(curr->m_right != NULL) { *stack++ = curr->m_right; } } else {//   .  . visit(curr->m_body_n, curr->m_mass_center, curr->m_mass); } } }

इस कार्यान्वयन के पास एक और समस्या है - इस तरह के पेड़ के ट्रैवर्सल को समानांतर करने का कोई सार्वभौमिक तरीका नहीं है। लेकिन हम पूरी तरह से पेड़ को स्मृति में संग्रहीत करने के तरीके को बदल सकते हैं - हम सभी नोड्स को एक रैखिक सरणी में संग्रहीत कर सकते हैं और पूरी तरह से बाइनरी हीप के निर्माण के साथ सादृश्य द्वारा, वंश के भंडारण का त्याग कर सकते हैं । जब नोड्स की संख्या के साथ शुरू करना

$1$ यदि नोड सेल नंबर में है

$i$ , तब उसका बाया बच्चा सेल में है

$2i$ , सेल में सही बच्चा है

$2i + 1$ , और सेल में जनक

$i/2$ । संख्या के साथ बाईं ओर दाएं नोड

$i$ एक नंबर होगा

$i + 1$ । सरणी में ही एक लंबाई होगी

$2N$ , और सभी पत्ती नोड्स आखिरी में स्थित होंगे

$N$ सेल, इसके अलावा, बारीकी से स्थित नोड्स सरणी के करीबी कोशिकाओं में स्थित होंगे। ट्री ट्रैवर्सल फ़ंक्शन थोड़ा बदल जाएगा, लेकिन फ्रेमवर्क समान रहेगा:

एक सरणी में एक पेड़ को काटकर बल की गणना

 nbvertex_t force_compute(const nbvertex_t& v1, const nbcoord_t mass1) { nbvertex_t total_force; size_t stack_data[MAX_STACK_SIZE] = {}; size_t* stack = stack_data; size_t* stack_head = stack; *stack++ = NBODY_HEAP_ROOT_INDEX; while(stack != stack_head) { size_t curr = *--stack; const nbcoord_t distance_sqr((v1 - m_mass_center[curr]).norm()); if(distance_sqr > m_radius_sqr[curr]) { total_force += force(v1, m_mass_center[curr], mass1, m_mass[curr]); } else { size_t left(left_idx(curr)); size_t rght(rght_idx(curr)); if(rght < m_body_n.size()) { *stack++ = rght; } if(left < m_body_n.size()) { *stack++ = left; } } } return total_force; }

लेकिन यह सभी संभावनाएं नहीं है कि सरणी में नोड्स को संग्रहीत करना हमारे लिए खुलता है - हम रेंगने पर स्टैक को मना कर सकते हैं। ऐसा करने के लिए, कोड शाखा जिसमें हम नोड के बच्चों के पास जाते हैं, हम अगले नोड की गणना करने के लिए फ़ंक्शन जोड़ते हैं (

$next_{up}$ ), और उस शाखा में, जिसमें हम इंटरैक्शन फोर्स की गणना करते हैं, हम अगले नोड की गणना को वर्तमान सबट्री स्कीप के साथ जोड़ते हैं (

$skip_{down}$ )।

वर्तमान सबट्री को छोड़ने के लिए, हमें रूट (दिशा) में जाने की आवश्यकता है

$D$ ), जबकि हम दाएं वंश में हैं, जैसे ही हम बाईं ओर पहुंचते हैं, हम संबंधित दाएं उप-वर्ग (दिशा) में जाते हैं

$R$ ), यदि हम जड़ तक पहुंचते हैं, तो पेड़ का चक्कर पूरा हो जाता है।

सबट्री स्किप फंक्शन कोड

 index_t skip_down(index_t idx) { // While index is 'right' -> go down while(is_right(idx)) { index_t parent = parent_idx(idx); // We at root again. Stop traverse. if(parent == NBODY_HEAP_ROOT_INDEX) { return NBODY_HEAP_ROOT_INDEX; } idx = parent; } return left2right(idx); }

चित्रा 8. एक उपप्रकार लंघन

$skip_{down}$ ।

अगले नोड पर जाने के लिए, यदि संभव हो तो, बाएं बच्चे पर जाएं (दिशा)

$यू$ ), और अगर कोई वंशज नहीं है, तो फ़ंक्शन का उपयोग करके 'नीचे' से अगले नोड पर जाएं

$skip_{down}$ ।

अगले नोड फ़ंक्शन कोड पर जाएं

 index_t next_up(index_t idx, index_t tree_size) { index_t left = left_idx(idx); if(left < tree_size) { return left; } return skip_down(idx); }

चित्रा 9. अगले नोड के लिए संक्रमण

$next_{up}$ ।

ऐसा लग सकता है कि हमने रिकर्सन को लूप से बदल दिया है

$while$ समारोह में

$skip_{down}$ , और ऐसा प्रतिस्थापन कुछ भी नहीं देता है, लेकिन आइए देखें कि किसी संख्या के साथ नोड कैसे निर्धारित करें

$i$ सही वंशज। यह बस इसकी विषम संख्या (सही बच्चा सेल में है) की जाँच करके किया जा सकता है

$2i + 1$ ), इसके लिए यह गणना करने के लिए पर्याप्त है

$i\&1$ । यानी हम संख्या को विभाजित करते हैं

$i$ दो अगर कम से कम महत्वपूर्ण बिट एक पर सेट है। लेकिन यह एक लूप के बिना किया जा सकता है, कई प्रोसेसर में एक पहला सेट निर्देश होता है जो पहले सेट बिट की स्थिति देता है, इसका उपयोग करते हुए हम लूप को चार प्रोसेसर निर्देशों में बदलते हैं।

ऑप्टिमाइज़्ड स्किप सबट्री फंक्शन कोड

 index_t skip_down(index_t idx) { idx = idx >> (__builtin_ffs(~idx) - 1); return left2right(idx); }

उसके बाद, हम ट्री ट्रैवर्सल फ़ंक्शन से स्टैक को बाहर कर सकते हैं और इसे एक जोड़ी के साथ बदल सकते हैं

$skip_{down} + next_{up}$ , उसके बाद फ़ंक्शन और भी सरल दिखता है:

एक ढेर का उपयोग किए बिना एक सरणी में एक पेड़ को पीछे हटाने के द्वारा बल की गणना

 nbvertex_t force_compute(const nbvertex_t& v1, const nbcoord_t mass1) const { nbvertex_t total_force; size_t curr = NBODY_HEAP_ROOT_INDEX; size_t tree_size = m_mass_center.size(); do { const nbcoord_t distance_sqr((v1 - m_mass_center[curr]).norm()); if(distance_sqr > m_radius_sqr[curr]) { total_force += force(v1, m_mass_center[curr], mass1, m_mass[curr]); curr = skip_down(curr); } else { curr = next_up(curr, tree_size); } } while(curr != NBODY_HEAP_ROOT_INDEX); return total_force; }

कुल में, हमें पेड़ की छल और बल गणना के छह संयोजन मिले। गणना समय और गुणवत्ता के संदर्भ में इन दृष्टिकोणों की तुलना करें। हम 100 पुनरावृत्तियों के बाद प्रणाली की कुल ऊर्जा में सापेक्ष परिवर्तन की गुणवत्ता के उपाय के रूप में लेते हैं। एक मॉडल प्रणाली के रूप में, हम दो परस्पर क्रिया "आकाशगंगाओं" से मिलकर बनाते हैं

$16384$ प्रत्येक निकाय।

तालिका 4. पेड़ की त्रैमासिक विधि और बल गणना का संयोजन

ट्री ट्रैवर्सल / बल गणना प्रकार	ढेर के साथ पेड़	ढेर के साथ 'ढेर'	ढेर के बिना 'हीप'
शरीर की संख्या द्वारा परिवर्तन	चक्र + वृक्ष	चक्र + ढेर	चक्र + ढेर
पत्ती बाईपास	nestedtree + वृक्ष	नेस्टट्री + ढेर	नेस्टट्री + हीपस्टैकलेस


ए) फ़ंक्शन गणना समय निर्भरता $f_n$ पेड़ के नोड के लिए इसकी त्रिज्या के लिए महत्वपूर्ण दूरी के अनुपात से ( $\lambda_{crit}$ )।	ख) ट्री नोड के त्रिज्या के लिए महत्वपूर्ण दूरी के अनुपात पर सिस्टम में ऊर्जा में सापेक्ष परिवर्तन की निर्भरता ( $\lambda_{crit}$ )।
चित्रा 10. समारोह गणना परिणाम $f_n$ सीपीयू पर विभिन्न तरीकों से (निकायों की संख्या) $N=32768$ )

यह देखा जा सकता है कि 'नेस्टट्री + ट्री' के क्रियान्वयन की उम्मीद तेजी से पीछे है, क्योंकि इसमें संक्षिप्तता का अभाव है। सरणी में ट्री नोड्स के स्थान के साथ कार्यान्वयन और एक बाइनरी हीप में अनुक्रमण प्रमुख हैं। ऊर्जा में सापेक्षिक परिवर्तन सभी प्रकारों के लिए नगण्य है

$\lambda_{crit} > 1$ ।अंजीर में भी। 10a से पता चलता है कि (के लिए

$\lambda_{crit} < 10$ ) फ़ंक्शन गणना के सभी प्रकार

$f_n$ आगे की गणना के साथ सटीक गणना ('ओपनम्प + ब्लॉक + ऑप्टिमाइज़ेशन') के सबसे अनुकूलित संस्करण को गति से आगे ले जाएं

$\lambda_{crit}$ एक पेड़ के साथ कार्यान्वयन सटीक संस्करण खो देते हैं।

GPU पेड़ ट्रैवर्सल

मैंने OpenCL तकनीक और CUDA दोनों का उपयोग करके GPU पर पेड़ को बायपास करने की कोशिश की। एक पेड़ के रूप में नोड्स को संग्रहीत करने का विकल्प तुरंत खारिज कर दिया गया था, और केवल विकल्पों को एक बाइनरी हीप में अनुक्रमण के साथ एक सरणी में पेड़ को संग्रहीत करने के साथ छोड़ दिया गया था। सामान्य तौर पर, कंप्यूटिंग कोर का कार्यान्वयन सीपीयू संस्करण से बहुत अलग नहीं है।

OpenCL कर्नेल एक पेड़ (नंबरिंग निकायों के क्रम में ट्रावर्सल) का पता लगाने के लिए कंप्यूटिंग शक्ति के लिए

 __kernel void ComputeTreeBH(int offset_n1, int points_count, int tree_size, __global const nbcoord_t* y, __global nbcoord_t* f, __global const nbcoord_t* tree_cmx, __global const nbcoord_t* tree_cmy, __global const nbcoord_t* tree_cmz, __global const nbcoord_t* tree_mass, __global const nbcoord_t* tree_crit_r2) { int n1 = get_global_id(0) + offset_n1; int stride = points_count; __global const nbcoord_t* rx = y; __global const nbcoord_t* ry = rx + stride; __global const nbcoord_t* rz = rx + 2 * stride; __global const nbcoord_t* vx = rx + 3 * stride; __global const nbcoord_t* vy = rx + 4 * stride; __global const nbcoord_t* vz = rx + 5 * stride; __global nbcoord_t* frx = f; __global nbcoord_t* fry = frx + stride; __global nbcoord_t* frz = frx + 2 * stride; __global nbcoord_t* fvx = frx + 3 * stride; __global nbcoord_t* fvy = frx + 4 * stride; __global nbcoord_t* fvz = frx + 5 * stride; nbcoord_t x1 = rx[n1]; nbcoord_t y1 = ry[n1]; nbcoord_t z1 = rz[n1]; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; int stack_data[MAX_STACK_SIZE] = {}; int stack = 0; int stack_head = stack; stack_data[stack++] = NBODY_HEAP_ROOT_INDEX; while(stack != stack_head) { int curr = stack_data[--stack]; nbcoord_t dx = x1 - tree_cmx[curr]; nbcoord_t dy = y1 - tree_cmy[curr]; nbcoord_t dz = z1 - tree_cmz[curr]; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 > tree_crit_r2[curr]) { if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = tree_mass[curr] / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; res_x -= dx; res_y -= dy; res_z -= dz; } else { int left = left_idx(curr); int rght = rght_idx(curr); if(left < tree_size) { stack_data[stack++] = left; } if(rght < tree_size) { stack_data[stack++] = rght; } } } frx[n1] = vx[n1]; fry[n1] = vy[n1]; frz[n1] = vz[n1]; fvx[n1] = res_x; fvy[n1] = res_y; fvz[n1] = res_z; }

पहले संस्करण में, वृक्षों की कटाई मूल सरणी में शवों को क्रमबद्ध करने के क्रम में शुरू हुई, इसलिए पड़ोसी धागे ने पेड़ के पूरी तरह से अलग हिस्सों को दरकिनार कर दिया, जिसने GPU मेमोरी कैश के प्रदर्शन को नकारात्मक रूप से प्रभावित किया। इसलिए, दूसरे अवतार में, पेड़ के ऊपर से शुरू होने वाला एक ट्रैवर्सल लगाया गया था, इस मामले में, पड़ोसी धागे उसी रास्ते से पेड़ को पीछे हटाना शुरू कर देते हैं, क्योंकि पड़ोसी पेड़ के शीर्ष पास और अंतरिक्ष में हैं। यह भी महत्वपूर्ण है कि हमने स्क्रैच से नहीं बल्कि ट्री नोड्स के एरे में नंबरिंग को चुना है, लेकिन एक से, इस मामले में ट्री की पत्तियों को सरणी के दूसरे भाग में संग्रहीत किया जाता है, और दो की शक्ति के बराबर निकायों की संख्या के साथ, हम सूचकांक tn1 द्वारा मेमोरी के बराबर उपयोग करेंगे।

OpenCL कर्नेल एक पेड़ की कटाई करके ताकत की गणना के लिए (एक पेड़ के नंबरिंग नोड्स का पता लगाने)

 __kernel void ComputeHeapBH(int offset_n1, int points_count, int tree_size, __global const nbcoord_t* y, __global nbcoord_t* f, __global const nbcoord_t* tree_cmx, __global const nbcoord_t* tree_cmy, __global const nbcoord_t* tree_cmz, __global const nbcoord_t* tree_mass, __global const nbcoord_t* tree_crit_r2, __global const int* body_n) { int tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX; int stride = points_count; int tn1 = get_global_id(0) + offset_n1 + tree_offset; int n1 = body_n[tn1]; nbcoord_t x1 = tree_cmx[tn1]; nbcoord_t y1 = tree_cmy[tn1]; nbcoord_t z1 = tree_cmz[tn1]; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; int stack_data[MAX_STACK_SIZE] = {}; int stack = 0; int stack_head = stack; stack_data[stack++] = NBODY_HEAP_ROOT_INDEX; while(stack != stack_head) { int curr = stack_data[--stack]; nbcoord_t dx = x1 - tree_cmx[curr]; nbcoord_t dy = y1 - tree_cmy[curr]; nbcoord_t dz = z1 - tree_cmz[curr]; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 > tree_crit_r2[curr]) { if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = tree_mass[curr] / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; res_x -= dx; res_y -= dy; res_z -= dz; } else { int left = left_idx(curr); int rght = rght_idx(curr); if(left < tree_size) { stack_data[stack++] = left; } if(rght < tree_size) { stack_data[stack++] = rght; } } } __global const nbcoord_t* vx = y + 3 * stride; __global const nbcoord_t* vy = y + 4 * stride; __global const nbcoord_t* vz = y + 5 * stride; __global nbcoord_t* frx = f; __global nbcoord_t* fry = frx + stride; __global nbcoord_t* frz = frx + 2 * stride; __global nbcoord_t* fvx = frx + 3 * stride; __global nbcoord_t* fvy = frx + 4 * stride; __global nbcoord_t* fvz = frx + 5 * stride; frx[n1] = vx[n1]; fry[n1] = vy[n1]; frz[n1] = vz[n1]; fvx[n1] = res_x; fvy[n1] = res_y; fvz[n1] = res_z; }

ट्री नोड्स की संख्या के क्रम को ट्रेस करने पर, हमें एक प्रदर्शन में वृद्धि हुई। लेकिन इस विकल्प में सुधार भी किया जा सकता है। वैश्विक स्मृति जिसमें ट्री नोड्स वर्तमान में स्थित हैं , साझा पहुंच के लिए अनुकूलित है , अर्थात। एक समूह के थ्रेड्स को एक मेमोरी ब्लॉक में स्थित शब्दों को पढ़ना चाहिए। हमारे मामले में, ट्री ट्रैवर्सल एक ही रास्तों के साथ शुरू होता है, और हम समूह के सभी थ्रेड्स के साथ एक ही डेटा का अनुरोध करते हैं, लेकिन आगे हम पेड़ में गहराई तक जाते हैं, पड़ोसी धागे के रास्ते अधिक से अधिक मोड़ते हैं, और हमें अलग-अलग डेटा का अनुरोध करना पड़ता है, जो कई बार सबसिस्टम के प्रदर्शन को कम कर देता है। स्मृति। लेकिन एक ही स्तर से संबंधित प्रत्येक उपप्रकार के नोड अपेक्षाकृत करीब स्मृति कोशिकाओं में स्थित हैं। यानीजब पेड़ के बाकी हिस्सों को ट्रेस किया जाता है, तो कंप्यूटिंग कोर के पड़ोसी धागे पेड़ के समान नोड्स तक नहीं पहुंचते हैं, लेकिन स्मृति में बारीकी से स्थित होते हैं। इस मेमोरी एक्सेस को ऑप्टिमाइज़ करने के लिए टेक्सचर मेमोरी का इस्तेमाल किया जा सकता है। लेकिन एक रोड़ा है। फिलहाल, बनावट दोहरे सटीक डेटा के साथ काम करने का समर्थन नहीं करती है (हम सही गणना करना चाहते हैं)। लेकिन CUDA में एक __hiloint2double फ़ंक्शन है जो दो पूर्णांक से एक डबल-सटीक संख्या एकत्र करता है।

एक बनावट भंडारण पूर्णांक से डबल सटीक संख्या अनुरोध कोड

 template<> struct nb1Dfetch<double> { typedef double4 vec4; static __device__ double fetch(cudaTextureObject_t tex, int i) { int2 p(tex1Dfetch<int2>(tex, i)); return __hiloint2double(py, px); } static __device__ vec4 fetch4(cudaTextureObject_t tex, int i) { int ii(2 * i); int4 p1(tex1Dfetch<int4>(tex, ii)); int4 p2(tex1Dfetch<int4>(tex, ii + 1)); vec4 d4 = {__hiloint2double(p1.y, p1.x), __hiloint2double(p1.w, p1.z), __hiloint2double(p2.y, p2.x), __hiloint2double(p2.w, p2.z) }; return d4; } };

इस स्थिति में, दो कार्यान्वयन किए गए थे, पेड़ के प्रत्येक तत्व (एक्स, वाई, जेड, ट्री_क्रिट_आर 2) में स्वतंत्र रूप से अनुरोध किया गया था, और दूसरे कार्यान्वयन में इन अनुरोधों को संयुक्त किया गया था। नोड के द्रव्यमान के लिए अनुरोध बहुत कम बार होता है, केवल अगर हालत r2> tree_crit_r2 [वक्र] पूरा हो गया है , तो यह दूसरों के साथ इस अनुरोध को संयोजित करने के लिए कोई मतलब नहीं है। CUDA फ्रेमवर्क की एक अन्य उपयोगी विशेषता L1 कैश के आकार और साझा मेमोरी ( cudaFuncSetCacheConfig ) के आकार के अनुपात को नियंत्रित करने की क्षमता है । ट्री ट्रैवर्सल के मामले में, हम साझा मेमोरी का उपयोग नहीं करते हैं, इसलिए हम L1 कैश को इसके नुकसान के लिए बढ़ा सकते हैं।

CUDA कर्नेल एक पेड़ की कटाई करके ताकत की गणना करने के लिए (पेड़ के नोड्स के क्रम क्रम का पता लगाता है)

 __global__ void kfcompute_heap_bh_tex(int offset_n1, int points_count, int tree_size, nbcoord_t* f, cudaTextureObject_t tree_xyzr, cudaTextureObject_t tree_mass, const int* body_n) { nb1Dfetch<nbcoord_t> tex; int tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX; int stride = points_count; int tn1 = blockDim.x * blockIdx.x + threadIdx.x + offset_n1 + tree_offset; int n1 = body_n[tn1]; nbvec4_t xyzr = tex.fetch4(tree_xyzr, tn1); nbcoord_t x1 = xyzr.x; nbcoord_t y1 = xyzr.y; nbcoord_t z1 = xyzr.z; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; int stack_data[MAX_STACK_SIZE] = {}; int stack = 0; int stack_head = stack; stack_data[stack++] = NBODY_HEAP_ROOT_INDEX; while(stack != stack_head) { int curr = stack_data[--stack]; nbvec4_t xyzr2 = tex.fetch4(tree_xyzr, curr); nbcoord_t dx = x1 - xyzr2.x; nbcoord_t dy = y1 - xyzr2.y; nbcoord_t dz = z1 - xyzr2.z; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 > xyzr2.w) { if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = tex.fetch(tree_mass, curr) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; res_x -= dx; res_y -= dy; res_z -= dz; } else { int left = nbody_heap_func<int>::left_idx(curr); int rght = nbody_heap_func<int>::rght_idx(curr); if(left < tree_size) { stack_data[stack++] = left; } if(rght < tree_size) { stack_data[stack++] = rght; } } } f[n1 + 3 * stride] = res_x; f[n1 + 4 * stride] = res_y; f[n1 + 5 * stride] = res_z; }

एनवीटीएफ प्रोफाइलर में कार्यक्रम के विश्लेषण से पता चला कि पेड़ को स्टोर करने के लिए बनावट मेमोरी का उपयोग करते समय, वैश्विक मेमोरी पर अभी भी बहुत अधिक भार है।

वास्तव में, CUDA में, 'गणना' पतों को संबोधित की जाने वाली सभी कर्नेल मेमोरी को वैश्विक मेमोरी में संग्रहीत किया जाता है, और तदनुसार स्टैक जो पेड़ को पार करने के लिए आवश्यक है, वैश्विक मेमोरी में स्थित है और मेमोरी चिप बैंडविड्थ के एक महत्वपूर्ण हिस्से को खा जाता है, क्योंकि स्टैक है प्रत्येक रन थ्रेड, और बहुत सारे थ्रेड हैं।

लेकिन, सौभाग्य से, हम पहले से ही जानते हैं कि एक ढेर का उपयोग किए बिना किसी पेड़ को कैसे पार करना है। अगले पेड़ के नोड की गणना के कार्यों के साथ पिछले कम्प्यूटेशनल कोर को लागू करते हुए, हमें एक नया कर्नेल, इसके अलावा, अधिक कॉम्पैक्ट मिलता है।

एक स्टैक का उपयोग किए बिना किसी पेड़ को ट्रेस करके कंप्यूटिंग शक्ति के लिए CUDA कोर

 __global__ void kfcompute_heap_bh_stackless(int offset_n1, int points_count, int tree_size, nbcoord_t* f, cudaTextureObject_t tree_xyzr, cudaTextureObject_t tree_mass, const int* body_n) { nb1Dfetch<nbcoord_t> tex; int tree_offset = points_count - 1 + NBODY_HEAP_ROOT_INDEX; int stride = points_count; int tn1 = blockDim.x * blockIdx.x + threadIdx.x + offset_n1 + tree_offset; int n1 = body_n[tn1]; nbvec4_t xyzr = tex.fetch4(tree_xyzr, tn1); nbcoord_t x1 = xyzr.x; nbcoord_t y1 = xyzr.y; nbcoord_t z1 = xyzr.z; nbcoord_t res_x = 0.0; nbcoord_t res_y = 0.0; nbcoord_t res_z = 0.0; int curr = NBODY_HEAP_ROOT_INDEX; do { nbvec4_t xyzr2 = tex.fetch4(tree_xyzr, curr); nbcoord_t dx = x1 - xyzr2.x; nbcoord_t dy = y1 - xyzr2.y; nbcoord_t dz = z1 - xyzr2.z; nbcoord_t r2 = (dx * dx + dy * dy + dz * dz); if(r2 > xyzr2.w) { if(r2 < NBODY_MIN_R) { r2 = NBODY_MIN_R; } nbcoord_t r = sqrt(r2); nbcoord_t coeff = tex.fetch(tree_mass, curr) / (r * r2); dx *= coeff; dy *= coeff; dz *= coeff; res_x -= dx; res_y -= dy; res_z -= dz; curr = nbody_heap_func<int>::skip_idx(curr); } else { curr = nbody_heap_func<int>::next_up(curr, tree_size); } } while(curr != NBODY_HEAP_ROOT_INDEX); f[n1 + 3 * stride] = res_x; f[n1 + 4 * stride] = res_y; f[n1 + 5 * stride] = res_z; }

GPU पर चलने वाले कोर का प्रदर्शन उन ब्लॉकों के आकार पर अत्यधिक निर्भर करता है जिनमें हम कार्य को विभाजित करते हैं। इस आकार पर निर्भर करता है कि प्रत्येक कंप्यूटिंग थ्रेड के लिए कितने रजिस्टर, स्थानीय मेमोरी और अन्य संसाधन उपलब्ध होंगे। आपको यह भी ध्यान रखना होगा कि एक थ्रेड में मेमोरी एक्सेस के लिए प्रतीक्षा करते समय, दूसरा थ्रेड shader प्रोसेसर पर गणना कर सकता है, इस प्रकार, एक प्रोसेसर पर पर्याप्त संख्या में एक साथ चलने वाले थ्रेड्स के साथ, मेमोरी एक्सेस टाइम गणनाओं के पीछे छिपा होगा। इसलिए, हमारे कोर के प्रदर्शन की तुलना करने से पहले, हमें उनमें से प्रत्येक के लिए इष्टतम ब्लॉक आकार की गणना करने की आवश्यकता है। आइए हम एनवीडिया K80 से आधे पर उपलब्ध तुलना करें।

तालिका 5. गणना समय (सेकंड में) फ़ंक्शन की निर्भरता

$f_n$ GPU

$N=131072$

$\lambda_{crit}=10$

/	8	16	32	64	128	256	512	1024
opencl+dense	5.77	2.84	1.46	1.13	1.15	1.14	1.14	1.13
cuda+dense	5.44	2.55	1.27	0.96	0.97	0.99	0.99	-
opencl+heap+cycle	5.88	5.65	5.74	5.96	5.37	5.38	5.35	5.38
opencl+heap+nested	4.54	3.68	3.98	5.25	5.46	5.41	5.48	5.31
cuda+heap+nested	3.62	2.81	2.68	4.26	4.84	4.75	4.8	4.67
cuda+heap+nested+tex	2.6	1.51	0.912	0.7	1.85	1.75	1.69	1.61
cuda+heap+nested+tex+stackless	2.3	1.29	0.773	0.5	0.51	0.52	0.52	0.52

6. ( )

$f_n$ GPU

$N=1M$

$\lambda_{crit}=4$

/	8	16	32	64	128	256	512	1024
opencl+dense	366	179	89.9	69.3	70.3	69.1	68.9	68.0
cuda+dense	346	162	79.6	60.8	60.8	60.7	59.6	-
opencl+heap+cycle	16.2	18.2	20.1	21.2	21.2	21.3	21.2	21.1
opencl+heap+nested	10.5	7.63	6.38	8.23	9.95	9.89	9.65	9.66
cuda+heap+nested	8.67	6.44	5.39	5.93	8.65	8.61	8.41	8.27
cuda+heap+nested+tex	6.38	3.57	2.13	1.44	3.56	3.46	3.30	3.29
cuda+heap+nested+tex+stackless	5.78	3.19	1.83	1.21	1.11	1.10	1.11	1.13

एक कठिन स्थिति, लेकिन, ट्री ट्रैवर्सल के सीपीयू संस्करण के विपरीत, यह स्पष्ट है कि प्रत्येक अनुकूलन कदम मूर्त परिणाम लाता है। पूर्ण समारोह गणना के साथ एक सटीक समाधान की तुलना में 'ओपेंसेल + हीप + चक्र' का कार्यान्वयन लगभग 6 गुना धीमा है

$f_n$ ।'ऑपनकल + हीप + नेस्टेड' का क्रियान्वयन, जिसमें एक पेड़ की कटाई पड़ोसी नोड्स से शुरू होती है, पिछले एक की तुलना में 1.4 गुना तेज, क्योंकि कैश मेमोरी का बेहतर उपयोग किया जाता है। + क्यूडा + हीप + नेस्टेड ’के कार्यान्वयन में, एल 1 कैश को साझा मेमोरी की वृद्धि के लिए बढ़ाया गया था, जिससे गति में 1.4 गुना की वृद्धि हुई थी, हालांकि यह संभव है कि कोडा कार्यान्वयन में संकलन कर्नेल अधिक आशावादी संकलित हो (enc opencl + dense’ और uda cuda के संस्करणों में) + घने 'कोर समान हैं, और कोडा का संस्करण प्रदर्शन ~ 1.2 गुना अधिक है)। कम्प्यूटेशनल गति (3.8 गुना) में सबसे महत्वपूर्ण वृद्धि तब प्राप्त होती है जब पेड़ बनावट की स्मृति में स्थित होता है और ट्री नोड के तत्वों के लिए प्रश्न संयुक्त होते हैं। 'कोडा + ढेर + नेस्टेड + टेक्स + स्टैकलेस' स्टैक का उपयोग किए बिना ट्री ट्रैवर्सल के साथ कार्यान्वयन 1.4 गुना तेज है।इसमें, मेमोरी बस के पूरे बैंडविड्थ का उपयोग केवल पेड़ के नोड्स के डेटा तक पहुंचने के लिए किया जाता है और स्टैक पर खर्च नहीं किया जाता है। इस प्रकार, के साथ

$\lambda_{crit}=10$ फ़ंक्शन की पूर्ण गणना की तुलना में दो बार त्वरण प्राप्त करना संभव था

$f_n$ । लेकिन

$\lambda_{crit}=10$ पैरामीटर के एक अत्यधिक बड़े मूल्य, सीपीयू कार्यान्वयन के लिए पेड़ के नोड से इसकी त्रिज्या के लिए महत्वपूर्ण दूरी के अनुपात से सिस्टम में ऊर्जा में सापेक्ष परिवर्तन के ग्राफ पर, यह देखा जा सकता है कि निम्न मानों का उपयोग किया जा सकता है

$\lambda_{crit}$ सटीकता के एक दृश्य हानि के बिना। चलो अलग करने की कोशिश करते हैं

$\lambda_{crit}$ इष्टतम ब्लॉक आकारों के साथ जो हमने पिछले चरण में निर्धारित किया था।


क) $N=128K$	ख) $N=1M$
अंजीर। 11. फ़ंक्शन गणना समय की निर्भरता। $f_n$ पेड़ के नोड के लिए इसकी त्रिज्या के लिए महत्वपूर्ण दूरी के अनुपात से ( $\lambda_{crit}$ ) विभिन्न GPU ट्री वॉक कार्यान्वयन के लिए

यह देखा जा सकता है कि छोटे के लिए

$\lambda_{crit}$ फ़ंक्शन गणना के सभी तरीके

$f_n$ नजदीकी मानों पर जाएं, केडी-पेड़ के निर्माण और GPU के लिए डेटा तैयार करने के समय से निर्धारित होता है। इसके अलावा, एक पेड़ बनाने का समय कुल समय में महत्वपूर्ण योगदान देता है

$\lambda_{crit}\leq 4$ , तो यह समय उपेक्षित हो सकता है। यह ध्यान रखना दिलचस्प है कि कब

$N=128K$ प्रदर्शन फिर से पहुंचने पर सुधार होता है

$\lambda_{crit}=1024$ सबसे अधिक संभावना यह है कि इस तथ्य के कारण कि GPU के सभी धागे एक ही समय में एक ही पेड़ के चक्कर में घूमते हैं, जो L1 कैश के उपयोग में सुधार करता है और 'शाखा विचलन' को पूरी तरह से समाप्त कर देता है । यह भी देखा जाता है कि स्टैक के बिना कार्यान्वयन के लिए सबसे अच्छा प्रदर्शन (क्यूडा + हीप + नेस्टेड + टेक्स + स्टैकलेस), यह लगभग स्टैक के साथ संस्करण को बेहतर बनाता है।

$1.4-1.5$ बार। अन्य कार्यान्वयन कई बार धीमे होते हैं। परिणाम को मजबूत करने के लिए, हम एक नए आर्किटेक्चर के साथ GPU पर समय की गणना करेंगे।

GeForce GTX 1080 Ti (एकल सटीकता) पर परिणाम लॉन्च करें

7. ( )

$f_n$ GPU

$N=1M$

$\lambda_{crit}=4$

/	8	16	32	64	128	256	512	1024
opencl+dense	47.8	23.4	11.6	11.59	12.8	12.8	12.8	12.8
cuda+dense	49.0	23.8	11.9	11.8	11.7	11.7	11.7	11.7
opencl+heap+cycle	7.48	8.26	7.73	7.36	7.33	7.27	7.25	7.26
opencl+heap+nested	1.33	1.20	1.41	2.46	2.53	2.49	2.44	2.47
cuda+heap+nested	1.23	1.10	1.31	2.28	2.29	2.28	2.23	2.25
cuda+heap+nested+tex	0.88	0.68	0.654	1.6	1.6	1.6	1.6	1.6
cuda+heap+nested+tex+stackless	0.71	0.47	0.45	0.43	0.43	0.42	0.41	0.395

12.

$f_n$ (

$\lambda_{crit}$ ) GPU

गणना के लिए GeForce GTX 1080 Ti का उपयोग करते समय, पेड़ के बीच का अंतर एक स्टैक के साथ चलता है और एक स्टैक के बिना दो बार तक होता है, बशर्ते कि हम उस समय को नजरअंदाज करते हैं जब पेड़ बनाने में समय लगता है। यह तथ्य हमें यह सुनिश्चित करने के लिए प्रेरित करता है कि भविष्य में GPU के लिए पोर्टिंग और एक पेड़ का निर्माण किया जाए। 5-7 की तालिकाओं की तुलना से पता चलता है कि अलग-अलग संख्या में निकायों के लिए कोई एकल इष्टतम ब्लॉक आकार नहीं है, और इससे भी अधिक विभिन्न GPU आर्किटेक्चर के लिए, इसके अलावा, गणना समय अंतर कई बार तक पहुंच सकता है, भले ही आप सीमा मानों को ध्यान में न रखें। इस प्रकार, लंबी गणना से पहले, यह प्रत्येक कॉन्फ़िगरेशन के लिए इष्टतम ब्लॉक आकार निर्धारित करने के लिए समझ में आता है।

मुख्य चीज जो हमने हासिल की है, वह एक लाख से अधिक निकायों के जोड़ीदार अंतःक्रिया करने की क्षमता है (

$2^{20}$ ) एक पर एक उचित समय के लिए, नवीनतम GPU नहीं। नए GPU (टेस्ला V100) पर, जाहिरा तौर पर, उचित समय में लगभग दो मिलियन इंटरेक्टिंग बॉडी को प्रोसेस करना संभव होगा, क्योंकि उनका प्रदर्शन टेस्ला K80 के आधे से लगभग चार गुना अधिक है। यद्यपि यह संख्या बौनी आकाशगंगाओं जैसे कि छोटे मैगेलैनिक क्लाउड में भी सितारों की संख्या के साथ अतुलनीय है , लेकिन, मेरी राय में, प्रभावशाली है।

निष्कर्ष

विभेदक समीकरणों को हल करने के लिए एम्बेडेड तरीकों का उपयोग अपेक्षाकृत कम समय की लागत पर उच्च स्तर की सटीकता के साथ समान समस्याओं को हल करना संभव बनाता है, और युग्मक आकर्षण बलों की गणना के कार्य को अंजाम देने के लिए एल्गोरिदम के आवेदन से हमें अंतःक्रियात्मक निकायों के कॉलोसल संस्करणों को संसाधित करने की अनुमति मिलती है। इस प्रकार, "तीन-शरीर कार्य" से एलियंस के विपरीत, हम समस्या को हल करने में सक्षम हैं

$N$ निकायों, और शरीर की एक छोटी संख्या के लिए और पूरे स्टार समूहों के लिए, हालांकि सटीकता के कुछ नुकसान की कीमत पर।

दृश्य

उन लोगों के लिए जो अंत तक पढ़ चुके हैं, मैं निकायों के विकास के दृश्य के साथ कुछ और वीडियो दूंगा।

दो आकाशगंगाओं के टकराने का अनुकरण। निकायों की कुल संख्या 60 हजार है।

एक लाख तारों से मिलकर एक आकाशगंगा के विकास को मॉडलिंग करना। केंद्र में कुल 99% द्रव्यमान वाला एक निकाय है। एकल कण लगभग अप्रभेद्य हैं। पहले से ही तरल की एक बूंद में एक लहर की तरह। गति मॉड्यूल के अनुसार रंगीन। कम गति - नीला, मध्यम - हरा, उच्च - लाल। यह देखा जाता है कि केंद्र में गति अधिक होती है, और धीरे-धीरे किनारों तक कम हो जाती है, और भूमध्यरेखीय तल में सबसे कम।

एक लाख-सितारा आकाशगंगा के विकास का एक अधिक 'गतिशील' अनुकरण। केंद्र में कुल 10% द्रव्यमान वाला एक शरीर होता है। केंद्रीय निकाय बाकी को काफी कम प्रभावित करता है। सबसे पहले, "सितारे" अलग उड़ते हैं, फिर कई समूहों में वापस इकट्ठा होते हैं, और अंत में फिर से एक बड़ा क्लस्टर बनाते हैं।

मॉडलिंग की प्रक्रिया में, लगभग 5% "सितारों" ने प्रारंभिक क्षेत्र को पूरी तरह से छोड़ दिया।

10 वें सेकंड में, यह मौजूदा कार्टव्हील आकाशगंगा से बहुत मिलता जुलता है ।

प्रोजेक्ट कोड जीथब पर पाया जा सकता है ।

एन निकायों का कार्य या रसोई छोड़ने के बिना आकाशगंगा को कैसे उड़ाया जाए

घटना का भौतिकी

गणित

परिणाम

सीपीयू कार्यान्वयन

GPU कार्यान्वयन

OpenCL

CUDA

समारोह स्वीकृति fnf_n

GPU पेड़ ट्रैवर्सल

निष्कर्ष

दृश्य

More articles:

समारोह स्वीकृति $f_n$