المداعبة

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

في كثير من الأحيان العمل مع ما يسمى المصفوفات المتفرقة - تلك الأصفار التي لها أوامر بحجم أكبر من القيم الأخرى. هذا ، على سبيل المثال ، أمر لا مفر منه إذا كنت تتعامل مع معادلات تفاضلية جزئية أو بعض العمليات الأخرى التي ترتبط فيها العناصر الناشئة في العلاقات الخطية المحددة فقط بـ "الجيران". فيما يلي مثال محتمل لمصفوفة متناثرة لمعادلة بواسون أحادية البعد المعروفة في الفيزياء الكلاسيكية $$$- \ phi ^ {''} = f$ على شبكة موحدة (نعم ، على الرغم من عدم وجود الكثير من الأصفار فيها ، ولكن عند طحن الشبكة ستكون صحية!):

$$

$$\ تبدأ {pmatrix} 2 & -1 & 0 & 0 & 0 \\ -1 & 2 & -1 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & -1 & 2 & -1 \\ 0 & 0 & 0 & -1 & 2 \ end {pmatrix}$$

المصفوفات المقابلة لها - تلك التي لا تولي اهتماما لعدد الأصفار وتأخذ في الاعتبار جميع المكونات دون استثناء - تسمى الكثيفة.

غالبًا ما يتم تقديم المصفوفات المتناثرة ، لأسباب مختلفة ، بتنسيق عمود مضغوط - عمود متفرق مضغوط (CSC). يستخدم هذا الوصف صفيفين صحيحين ونقطة واحدة عائمة. دع المصفوفة لديها كل شيء $$$nnz_A$ غير الصفر و $$$N$ الأعمدة. يتم سرد عناصر المصفوفة في أعمدة من اليسار إلى اليمين ، دون استثناء. مجموعة الأولى $$$i_A$ أطوال $$$nnz_A$ يحتوي على أرقام الصفوف لمكونات المصفوفة غير الصفرية. المجموعة الثانية $$$j_A$ أطوال $$$N + 1 دولا$ يحتوي على ... لا ، ليس أرقام الأعمدة ، لأنه بعد ذلك ستكتب المصفوفة بتنسيق إحداثي (إحداثي) ، أو ثلاثي (ثلاثي). يحتوي المصفوفة الثانية على الأرقام التسلسلية لمكونات المصفوفة هذه التي تبدأ بها الأعمدة ، بما في ذلك عمود دمية إضافي في النهاية. وأخيرا ، المجموعة الثالثة $$$v_A$ أطوال $$$nnz_A$ يحتوي بالفعل على المكونات نفسها ، بالترتيب المقابل للصفيف الأول. هنا ، على سبيل المثال ، في ظل افتراض أن ترقيم الصفوف والأعمدة من المصفوفات يبدأ من الصفر ، لمصفوفة معينة

$$

$$A = \ تبدأ {pmatrix} 0 & 1 & 0 & 4 & -1 \\ 7 & 0 & 2.1 & 0 & 3 \\ 0 & 0 & 0 & 10 & 0 \ end {pmatrix}$$

سوف تكون صفائف $$$i_A = \ {1 ، 0 ، 1 ، 0 ، 2 ، 0 ، 1 \}$ ، $$$j_A = \ {0 ، 1 ، 2 ، 3 ، 5 ، 7 \}$ ، $$$v_A = \ {7 ، 1 ، 2.1 ، 4 ، 10 ، -1 ، 3 \}$ .

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

دعنا نقول لمصفوفة مربعة $$$دولا$ التحلل من النموذج متاح لنا $$$A = LU$ اين $$$L$ و $$$يو$ - على التوالي ، المصفوفات الثلاثية السفلية والمثلثية العليا ، على التوالي. الأول يعني أنه يحتوي على أصفار واحد أعلى المائل ، والثاني يعني أنه أسفل المائل. كيف بالضبط حصلنا على هذا التحلل - نحن لسنا مهتمين الآن. فيما يلي مثال بسيط للتحلل:

$$

$$\ تبدأ {pmatrix} 1 & -1 & -1 \\ 2 & - 1 & -0.5 \\ 4 & -2 & -1.5 \ end {pmatrix} = \ تبدأ {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ start {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5 \ end {pmatrix}$$

كيف في هذه الحالة لحل نظام المعادلات $$ على سبيل المثال ، مع الجانب الأيمن $$$f = \ تبدأ {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix}$ ؟ المرحلة الأولى هي خطوة مباشرة (حل للأمام = استبدال للأمام). نحن نشير $$$y: = Ux$ والعمل مع النظام $$$Ly = f$ . بسبب $$$L$ - انخفاض الثلاثي ، على التوالي في دورة نجد جميع المكونات $$$ذ$ أعلى إلى أسفل:

$$

$$\ تبدأ {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ تبدأ {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ تبدأ {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ يتضمن ص = \ تبدأ {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix}$$

الفكرة المركزية هي أنه عند العثور عليه $$$i$ مكونات المتجه $$$ذ$ يتم ضربه بواسطة عمود بنفس رقم المصفوفة $$$L$ الذي يطرح بعد ذلك من الجانب الأيمن. يبدو أن المصفوفة نفسها تنهار من اليسار إلى اليمين ، وتناقص في الحجم مع وجود المزيد والمزيد من مكونات المتجهات $$$ذ$ . وتسمى هذه العملية القضاء العمود.

المرحلة الثانية هي خطوة إلى الوراء (حل الوراء = استبدال الوراء). وجود ناقلات وجدت $$$ذ$ نحن نقرر $$$Ux = y$ . نحن هنا نتصفح المكونات من الأسفل إلى الأعلى ، لكن الفكرة تظل كما هي: $$$i$ يتم ضرب العمود ال بواسطة المكون الموجود للتو $$$x_i$ وينتقل إلى اليمين ، وتنهار المصفوفة من اليمين إلى اليسار. يتم توضيح العملية برمتها للمصفوفة المذكورة في المثال في الصورة أدناه:

$$

$$\ small \ تبدأ {pmatrix} 1 & 0 & 0 \\ 2 & 1 & 0 \\ 4 & 2 & 1 \ end {pmatrix} \ start {pmatrix} y_1 \\ y_2 \\ y_3 \ end {pmatrix} = \ تبدأ {pmatrix} 4 \\ 2 \\ 3 \ end {pmatrix} \ ضمني \ start {pmatrix} 1 & 0 \\ 2 & 1 \ end {pmatrix} \ تبدأ {pmatrix} y_2 \\ y_3 \ end {pmatrix } = \ تبدأ {pmatrix} -6 \\ -13 \ end {pmatrix} \ ضمني \ تبدأ {pmatrix} 1 \ end {pmatrix} \ start {pmatrix} y_3 \ end {pmatrix} = \ تبدأ {pmatrix} -1 \ end {pmatrix}$$

$$

$$\ small \ تبدأ {pmatrix} 1 & -1 & -1 \\ 0 & 1 & 1.5 \\ 0 & 0 & -0.5 \ end {pmatrix} \ start {pmatrix} x_1 \\ x_2 \\ x_3 \ end { pmatrix} = \ تبدأ {pmatrix} 4 \\ -6 \\ -1 \ end {pmatrix} \ Rightarrow \ start {pmatrix} 1 & -1 \\ 0 & 1 \ end {pmatrix} \ start {pmatrix} x_1 \ \ x_2 \ end {pmatrix} = \ start {pmatrix} 6 \\ -9 \ end {pmatrix} \ implies \ start {pmatrix} 1 \ end {pmatrix} \ start {pmatrix} x_1 \ end {pmatrix} = \ start {pmatrix} -3 \ end {pmatrix}$$

وقرارنا سيكون $$$x = \ تبدأ {pmatrix} -3 \\ - 9 \\ 2 \ end {pmatrix}$ .

إذا كانت المصفوفة كثيفة ، أي ، فهي ممثلة بالكامل في شكل نوع من الصفيف ، أحادي البعد أو ثنائي الأبعاد ، ويتم الوصول إلى عنصر محدد فيه بمرور الوقت $$$O (1)$ وبالتالي ، فإن إجراء حل مشابه مع التحلل الحالي ليس صعبًا ، كما أنه سهل الترميز ، لذلك لن نضيع الوقت فيه. ماذا لو كانت المصفوفة متناثرة؟ هنا ، من حيث المبدأ ، هو أيضا ليست صعبة. هنا هو رمز C ++ للحركة إلى الأمام التي الحل $$$x$ هو مكتوب للمكان على الجانب الأيمن ، دون التحقق من صحة بيانات الإدخال (صفائف CSC تتوافق مع المصفوفة الثلاثي السفلي):

الخوارزمية 1:

void forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { x[j] /= vA[jA[j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j] ; } } 

جميع المناقشات الإضافية سوف تتعلق فقط بالسكتة الدماغية الأمامية لحل النظام الثلاثي السفلي $$$Lx = f$ .

التعادل

ولكن ماذا لو كان الجانب الأيمن ، على سبيل المثال متجه إلى يمين علامة المساواة $$$Lx = f$ هل هو نفسه يحتوي على عدد كبير من الأصفار؟ ثم من المنطقي تخطي الحسابات المرتبطة بمراكز الصفر. التغييرات في الكود في هذه الحالة ضئيلة:

الخوارزمية 2:

 void fake_sparse_forward_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, std::vector<double>& x) { size_t j, p, n = x.size(); for (j = 0; j < n; j++) //    { if(x[j]) { x[j] /= vA[jA [j]]; for (p = jA[j]+1; p < jA[j+1]; p++) x[iA[p]] -= vA[p] * x[j]; } } } 

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

يبدو كل شيء رائعًا ، وبالطبع سيعمل ، لكن هنا ، بعد المداعبة الطويلة ، أصبحت المشكلة واضحة أخيرًا - الأداء المنخفض بشكل غير متماثل لهذا المحلل في حالة الأنظمة الكبيرة. وكل ذلك بسبب حقيقة وجود حلقة. ما هي المشكلة؟ حتى لو كانت الحالة داخل الحالة "نادرًا" نادرة للغاية ، فلا يوجد مفر من الدورة نفسها ، وهذا يخلق تعقيد الخوارزمية $$$O (N)$ اين $$$N$ هو حجم المصفوفة. بغض النظر عن كيفية تحسين الحلقات من قبل المترجمين الحديثين ، فإن هذا التعقيد سيجعل نفسه يشعر بأحجام كبيرة $$$N$ . الهجومية خاصة عندما تكون ناقلات كاملة $$$و$ يتكون كلياً من أصفار ، لأنه ، كما قلنا ، لا تفعل شيئًا في هذه الحالة! ليست عملية حسابية واحدة! ماذا بحق الجحيم $$$O (N)$ العمل؟

حسنا ، دعنا نقول. ومع ذلك ، لماذا لا يمكنك تحمل الخمول البعيد ، لأن الحسابات الحقيقية بأرقام حقيقية ، مثل تلك التي تندرج تحت if سوف لا تزال قليلة؟ ولكن الحقيقة هي أن هذا الإجراء التدفق المباشر مع الجانب الأيمن غريب نفسه غالبا ما تستخدم في الدورات الخارجية وتحت التحلل Cholesky $$$A = LL ^ T$ و LU المظهر الأيسر التحلل . نعم ، نعم ، أحد هذه التوسعات ، دون القدرة على القيام بكل هذه التحركات المباشرة والعكسية في حل الأنظمة الخطية ، يفقد كل معنى عملي.

نظرية إذا كانت المصفوفة محددة بشكل إيجابي متماثل (SPD) ، فيمكن تمثيلها كـ $$$A = LL ^ T$ الطريقة الوحيدة أين $$$L$ - انخفاض مصفوفة الثلاثي مع عناصر إيجابية على قطري.

بالنسبة لمصفوفات SPD المتفرق للغاية ، يتم استخدام تحليل Cholesky ذي المظهر العالي. يمثل بشكل تخطيطي التحلل في شكل كتلة المصفوفة

$$

$$\ تبدأ {pmatrix} \ mathbf {L} _ {11} & \ mathbf {0} \\ \ mathbf {l} _ {12} ^ T & l_ {22} \ end {pmatrix} \ start {pmatrix} \ mathbf {L} _ {11} ^ T & \ mathbf {l} _ {12} \\ \ mathbf {0} & l_ {22} \ end {pmatrix} = \ تبدأ {pmatrix} \ mathbf {A} _ { 11} & \ mathbf {a} _ {12} \\ \ mathbf {a} _ {12} ^ T & a_ {22} \ end {pmatrix}،$$

يمكن تقسيم عملية التعمير بأكملها منطقياً إلى ثلاث خطوات فقط.

الخوارزمية 3:

1. الطريقة العليا لتحلل Cholesky ذات البعد الأصغر $$$\ mathbf {L} _ {11} \ mathbf {L} _ {11} ^ T = \ mathbf {A} _ {11}$ (العودية!)
2. نظام الثلاثي السفلي مع الجانب الأيمن متفرق $$$\ mathbf {L} _ {11} \ mathbf {l} _ {12} = \ mathbf {a} _ {12}$ ( ها هو! )
3. حساب $$$l_ {22} = \ sqrt {a_ {22} - \ mathbf {l} _ {12} ^ T \ mathbf {l} _ {12}}$ (عملية تافهة)

في الممارسة العملية ، يتم تنفيذ هذا بطريقة أن يتم تنفيذ الخطوتين 3 و 2 في دورة واحدة كبيرة ، وبهذا الترتيب. وبالتالي ، يتم تشغيل المصفوفة قطريا من أعلى إلى أسفل ، وزيادة المصفوفة $$$L$ سطرا سطرا في كل تكرار للحلقة.

إذا كان في هذا aglorithm التعقيد إلى الأمام في الخطوة 2 سيكون $$$O (N)$ اين $$$N$ هو حجم المصفوفة الثلاثي السفلي $$$\ mathbf {L} _ {11}$ على التكرار التعسفي لدورة كبيرة ، فإن تعقيد التحلل بأكمله سيكون على الأقل $$$O (N ^ 2)$ ! أوه ، أنا لا أحب ذلك!

حالة من الفن

تعتمد العديد من الخوارزميات بطريقة ما على محاكاة الإجراءات البشرية في حل المشكلات. إذا أعطيت أي شخص نظامًا خطيًا أدنى مثلث ذو جانب يمين ، بحيث يكون 1-2 فقط غير صفري ، فإنه بالطبع يدير متجه الجانب الأيمن مع عينيه من أعلى إلى أسفل (تلك الدورة اللعينة من التعقيد $$$O (N)$ ! ) للعثور على هذه nonzeros. عندها سيستخدمها فقط ، دون إضاعة الوقت على المكونات الصفرية ، لأن هذا الأخير لن يؤثر على الحل: ليس من المنطقي تقسيم الأصفار إلى مكونات قطرية للمصفوفة ، وكذلك نقل العمود مضروبًا في الصفر إلى اليمين. هذه هي الخوارزمية أعلاه 2. لا توجد معجزات. ولكن ماذا لو تم إعطاء شخص على الفور أعداد المكونات غير الصفرية من بعض المصادر الأخرى؟ على سبيل المثال ، إذا كان الجانب الأيمن عبارة عن عمود في بعض المصفوفات الأخرى ، كما هو الحال مع تحلل Cholesky ، عندها يمكننا الوصول الفوري إلى مكوناته غير الصفرية ، إذا طلب ذلك بالتتابع:

 //     j,     : for (size_t p = jA[j]; p < jA[j+1]; p++) //    vA[p]     iA[p]   j 

تعقيد هذا الوصول هو $$$O (nnz_j)$ اين $$$nnz_j$ هو عدد المكونات غير الصفرية في العمود j. الحمد لله على تنسيق CSC! كما ترون ، لا يستخدم فقط لحفظ الذاكرة.

في هذه الحالة ، يمكننا أن نلاحظ جوهر ما يحدث عند الحل بطريقة الدورة التدريبية المباشرة $$$Lx = f$ لمصفوفات متفرق $$$L$ والجانب الأيمن $$$و$ . احبس أنفاسك: نأخذ مكونًا غير صفري $$$f_i$ على الجانب الأيمن ، نجد المتغير المقابل $$$x_i$ تقسيم على $$$L_ {ii}$ ومن ثم ، بضرب عمود ith بأكمله في هذا المتغير الذي تم العثور عليه ، فإننا نقدم عناصر غير صفرية إضافية على الجانب الأيمن بطرح هذا العمود من الجانب الأيمن! تم وصف هذه العملية جيدًا بلغة ... الرسوم البيانية. وعلاوة على ذلك ، الموجهة وغير دوري.

بالنسبة للمصفوفة الثلاثية السفلية التي لا تحتوي على أصفار على المائل ، نحدد رسم بياني متصل. نحن نفترض أن ترقيم الصف والعمود يبدأ من الصفر.
التعريف رسم بياني للاتصال لمصفوفة ثلاثية أقل $$$L$ الحجم $$$N$ ، والتي لا يوجد بها أصفار على القطر ، تسمى مجموعة العقد $$$V = \ {0، ...، N-1 \}$ والحواف الموجهة $$$E = \ {(j، i) | L_ {ij} \ ne 0 \}$ .

هنا ، على سبيل المثال ، هو الشكل الذي يبدو عليه الرسم البياني للاتصال لمصفوفة ثلاثية أقل.

حيث تتوافق الأرقام ببساطة مع الرقم الترتيبي للعنصر المائل ، وتشير النقاط إلى مكونات غير صفرية أسفل المائل:



التعريف إمكانية الوصول الرسم البياني المنحى $$$G$ على مؤشرات متعددة $$$W \ subset V$ استدعاء الكثير من المؤشرات $$$R_G (W) \ subset V$ هذا في أي مؤشر $$$z \ in R_G (W)$ يمكنك الحصول عليها من خلال تمرير الرسم البياني من بعض الفهرس $$$w (z) \ in W$ .

مثال: للحصول على رسم بياني من صورة $$$R_G (\ {0 ، 4 \}) = \ {0 ، 4 ، 5 ، 6 \}$ . من الواضح أنه يحمل دائمًا $$$W \ subset R_G (W)$ .

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

هيا $$$nnz (x) = \ {j | x_j \ ne 0 \}$ يدل على مجموعة من المؤشرات المقابلة لمواقع غير صفرية في المتجه x.

نظرية هلبرت (لا ، ليس باسمها يتم تسمية المسافات) $$$nnz (x)$ اين $$$x$ هناك متجه حل لنظام الثلاثي أقل متناثر $$$Lx = f$ مع الجانب الأيمن متفرق $$$و$ يتزامن مع $$$R_G (nnz (f))$ .

الإضافة من تلقاء نفسها: في النظرية ، لا نأخذ في الاعتبار الاحتمال المحتمل أن الأرقام غير الصفرية على الجانب الأيمن ، عند تدمير الأعمدة ، يمكن تخفيضها إلى صفر ، على سبيل المثال ، 3 - 3 = 0. يسمى هذا التأثير بالإلغاء العددي. إن أخذ مثل هذه الأصفار التي تحدث تلقائيًا في الاعتبار هو مضيعة للوقت ، ويُنظر إليها على أنها جميع الأرقام الأخرى في مواقع غير صفرية.

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

1. ننفذ الرسم البياني باستخدام طريقة "البحث أولاً عن العمق" ، بدءًا من الفهرس بالتتابع $$$i \ in nnz (f)$ كل مكون غير صفري من الجانب الأيمن. الكتابة وجدت العقد إلى مجموعة $$$R_G (nnz (f))$ أثناء القيام بذلك بالترتيب الذي عدنا فيه إلى الرسم البياني! بالقياس إلى جيش الغزاة: نحن نحتل البلد بدون قسوة ، لكن عندما بدأنا في العودة ، نحن الغاضبون ندمر كل شيء في طريقه.
   تجدر الإشارة إلى أنه ليس من الضروري على الإطلاق أن قائمة المؤشرات $$$nnz (f)$ تم الفرز حسب الصعود عندما تم تغذيته لإدخال خوارزمية "البحث في العمق الأول". يمكنك البدء في أي ترتيب على المجموعة $$$nnz (f)$ . تسلسل مختلف ينتمي إلى المجموعة $$$nnz (f)$ لا تؤثر الفهارس على الحل النهائي ، كما سنرى في المثال.
   
   هذه الخطوة لا تتطلب أي معرفة مجموعة حقيقية. $$$v_A$ ! مرحبًا بكم في عالم التحليل الرمزي عند العمل مع محللي متفرق مباشر!
   
2. ننتقل إلى الحل نفسه ، حيث لدينا تحت تصرفنا صفيف $$$R_G (nnz (f))$ من الخطوة السابقة. نحن ندمر الأعمدة بالترتيب العكسي لسجل هذه المجموعة . فويلا!

مثال

فكر في مثال يوضح فيه كلا الخطوتين بالتفصيل. افترض أن لدينا مصفوفة 12 × 12 من النموذج التالي:

الرسم البياني اتصال المقابلة لديه النموذج:

دع غير الصفر في الجزء الأيمن يكون في المواضع 3 و 5 فقط ، أي $$$nnz (f) = \ {3، 5 \}$ . دعنا نذهب من خلال الرسم البياني ، بدءا من هذه المؤشرات في ترتيب مكتوب. ثم ستبدو طريقة "البحث العميق" على النحو التالي. أولاً سنقوم بزيارة المؤشرات بالترتيب $$$3 \ إلى 8 \ إلى 11 دولا$ لا ننسى وضع علامات على الفهارس كما تمت زيارتها. في الشكل أدناه مظللة باللون الأزرق:

عند العودة إلى الوراء ، ندخل الأرقام القياسية في مجموعتنا من مكونات مكونات الحلول غير الصفرية $$$nnz (x): = \ {\ color {blue} {11} ، \ color {blue} 8 ، \ color {blue} 3 \}$ . بعد ذلك ، حاول الجري $$$5 \ to8 \ إلى ...$ ، لكننا صادفنا علامة 8 ، لذلك لا نلمس هذا الطريق ونذهب إلى الفرع $$$5 \ to9 \ to ...$ . ستكون نتيجة هذا المدى $$$5 \ to9 \ to10$ . لا يمكننا زيارة Node 11 ، لأنه مُسمى بالفعل. لذا ، عد واستكمل الصفيف $$$nnz (x)$ الصيد الجديد ملحوظ باللون الأخضر: $$$nnz (x): = \ {\ color {blue} {11} ، \ color {blue} 8 ، \ color {blue} 3 ، \ color {green} {10} ، \ color {green} 9 ، \ color {أخضر} 5 \}$ . وهنا هي الصورة:

العقدتان الأخضرتان الطنانتان 8 و 11 هما العقدتان اللتان أردنا زيارتهما خلال الجولة الثانية ، لكننا لم نستطع ذلك بسبب ذلك زار بالفعل خلال الأول.

لذلك مجموعة $$$nnz (x) = R_G (nnz (f))$ شكلت. نمر إلى الخطوة الثانية. خطوة بسيطة: نحن نمر من خلال مجموعة $$$nnz (x)$ بالترتيب العكسي (من اليمين إلى اليسار) ، والعثور على المكونات المقابلة لمتجه الحل $$$x$ القسمة على المكونات القطرية للمصفوفة $$$L$ ونقل الأعمدة إلى اليمين. مكونات أخرى $$$x$ كما كانوا أصفار ، ظلوا كذلك. بشكل تخطيطي ، يظهر هذا أدناه ، حيث يشير السهم السفلي إلى الترتيب الذي يتم به تدمير الأعمدة:

ملاحظة: بترتيب تدمير الأعمدة ، سيتم مواجهة الرقم 3 بعد الأرقام 5 و 9 و 10! لا يتم تدمير الأعمدة بترتيب تصاعدي ، مما سيكون خطأً بالنسبة للأعمدة الكثيفة ، أي المصفوفات غير المكتملة. لكن بالنسبة للمصفوفات المتفرقة ، هذا أمر شائع. كما يتضح من بنية المصفوفة غير الصفرية في هذا المثال ، فإن استخدام الأعمدة 5 و 9 و 10 إلى العمود 3 لن يشوه المكونات في الاستجابة $$$x_5 دولا$ ، $$$x_9 دولا$ ، $$$x_ {10}$ لا على الإطلاق ، لديهم ج $$$x_3 دولا$ فقط "لا التقاطعات". أخذت طريقتنا هذا في الاعتبار! وبالمثل ، فإن استخدام العمود 8 بعد العمود 9 لن يفسد المكون $$$x_9 دولا$ . خوارزمية رائعة ، أليس كذلك؟

إذا ذهبنا حول الرسم البياني أولاً من الفهرس 5 ، ثم من الفهرس 3 ، فستكون صفحتنا $$$nnz (x): = \ {\ color {blue} {11} ، \ color {blue} 8 ، \ color {blue} {10} ، \ color {blue} {9} ، \ color {blue} 5 ، \ اللون {الأخضر} 3 \}$ . إن تدمير الأعمدة بالترتيب العكسي لهذه المجموعة سيعطي نفس الحل تمامًا كما في الحالة الأولى!

يتم زيادة تعقيد جميع عملياتنا بعدد العمليات الحسابية الفعلية وعدد المكونات غير الصفرية على الجانب الأيمن ، ولكن ليس بحجم المصفوفة! لقد حققنا هدفنا.

النقد

ومع ذلك! سوف يلاحظ القارئ الناقد أن الافتراض ذاته في البداية بأن لدينا "إمكانية الوصول المباشر إلى المكونات غير الصفرية من الجانب الأيمن من الفهرس" يعني بالفعل أننا في وقت سابق ركضنا على الجانب الأيمن من أعلى إلى أسفل للعثور على هذه المؤشرات ذاتها وتنظيمها في الصفيف $$$nnz (f)$ ، وهذا هو ، لقد قضى بالفعل $$$O (N)$ العمل! علاوة على ذلك ، يتطلب تشغيل الرسم البياني نفسه أن نخصص أولاً ذاكرة الحد الأقصى للطول الممكن (في مكان آخر تحتاج إلى كتابة الفهارس التي تم العثور عليها من خلال البحث المتعمق!) حتى لا تعاني من مزيد من إعادة التخصيص مع نمو الصفيف $$$nnz (x)$ ، وهذا يتطلب أيضا $$$O (N)$ العمليات! لماذا ، إذن ، يقولون ، كل هذا العناء؟

ولكن في الواقع ، بالنسبة إلى حل لمرة واحدة لنظام مثلثي منخفض متفرق مع جانب يميني متفرق ، تم تعريفه أصلاً على أنه ناقل كثيف ، لا يوجد أي نقطة في إضاعة وقت مطور البرامج في جميع الخوارزميات المذكورة أعلاه. قد تكون أبطأ من طريقة الجبين المقدمة من الخوارزمية 2 أعلاه. ولكن ، كما ذكرنا سابقًا ، لا غنى عن هذا الجهاز لعامل Cholesky ، لذلك يجب ألا ترمي الطماطم في وجهي. في الواقع ، قبل تشغيل الخوارزمية 3 ، يتم تخصيص كل الذاكرة اللازمة ذات الطول الأقصى على الفور وتتطلب $$$O (N)$ الوقت في دورة عمود أخرى $$$دولا$ تتم الكتابة فوق كافة البيانات فقط في صفيف ذي طول ثابت $$$N$ وفقط في تلك المواضع التي تكون فيها إعادة الكتابة ذات صلة ، وذلك بفضل الوصول المباشر إلى العناصر غير الصفرية. ولهذا السبب بالتحديد تنشأ الكفاءة!

تنفيذ C ++

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

 /* j -   iA, jA -    ,    CSC top -     nnz(x);     0 result -   N,    nnz(x)   0  top-1  marked -    N;      */ void DepthFirstSearch(size_t j, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t& top, std::vector<size_t>& result, std::vector<int>& marked) { size_t p; marked[j] = 1; //   j   for (p = jA[j]; p < jA[j+1]; p++) //       j { if (!marked[iA[p]]) //  iA[p]   { DepthFirstSearch(iA[p], iA, jA, top, result, marked); //      iA[p] } } result[top++] = j; //  j   nnz(x) } 

إذا قمت بتمرير المتغير top يساوي الصفر في المكالمة الأولى إلى DepthFirstSearch ، وبعد الانتهاء من الإجراء DepthFirstSearch بأكمله ، فإن المتغير top يساوي عدد الفهارس الموجودة في صفيف result . الواجب المنزلي: اكتب وظيفة أخرى تأخذ مجموعة من مؤشرات المكونات غير الصفرية على الجانب الأيمن DepthFirstSearch بالتتابع إلى وظيفة DepthFirstSearch . بدون هذا ، الخوارزمية غير كاملة. يرجى ملاحظة: مجموعة من الأرقام الحقيقية $$$vA$ نحن لا ننقلها إلى الوظيفة على الإطلاق ، لأن ليست هناك حاجة في عملية التحليل الرمزي.

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

 /* j -   N -   iA, jA -    ,    CSC top -     nnz(x) result -   N,     nnz(x)   top  ret-1 ,  ret -    marked -    N work_stack -     N */ size_t DepthFirstSearch(size_t j, size_t N, const std::vector<size_t>& iA, const std::vector<size_t>& jA, size_t top, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, head = N - 1; int done; result[N - 1] = j; //    while (head < N) { j = result[head]; //  j     if (!marked[j]) { marked[j] = 1; //   j   work_stack[head] = jA[j]; } done = 1; //    j      for (p = work_stack[head] + 1; p < jA[j+1]; p++) //    j { //   c   iA[p] if (marked[iA[p]]) continue; //  iA[p]  ,   work_stack[head] = p; //       j result[--head] = iA[p]; //       iA[p] done = 0; //   j    break; } if (done) //      j  { head++; //  j    result[top++] = j; //  j    } } return (top); } 

وهذا هو ما يبدو عليه بالفعل مولد بنية ناقلات الحل غير الصفري $$$nnz (x)$ :

 /* iA, jA -   ,    CSC iF -        f nnzf -      f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N */ size_t NnzX(const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<size_t>& iF, size_t nnzf, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack) { size_t p, N, top; N = jA.size() - 1; top = 0; for (p = 0; p < nnzf; p++) if (!marked[iF[p]]) //        top = DepthFirstSearch(iF[p], N, iA, jA, vA, top, result, marked, work_stack); for (p = 0; p < top; p++) marked[result[p]] = 0; //   return (top); } 

أخيرًا ، بدمج كل شيء في واحد ، نكتب الحامل الثلاثي السفلي لحالة الجانب الأيمن المتناثر:

 /* iA, jA, vA -  ,    CSC iF -        f nnzf -      f vF -       f result -   N,    nnz(x)   0  ret-1 ,  ret -    marked -    N,    work_stack -     N x -    N,    ;     */ size_t lower_triangular_solve (const std::vector<size_t>& iA, const std::vector<size_t>& jA, const std::vector<double>& vA, const std::vector<size_t>& iF, size_t nnzf, const std::vector<double>& vF, std::vector<size_t>& result, std::vector<int>& marked, std::vector<size_t>& work_stack, std::vector<double>& x) { size_t top, p, j; ptrdiff_t px; top = NnzX(iA, jA, iF, nnzf, result, marked, work_stack); for (p = 0; p < nnzf; p++) x[iF[p]] = vF[p]; //    for (px = top; px > -1; px--) //     { j = result[px]; // x(j)   x [j] /= vA[jA[j]]; //   x(j) for (p = jA[j]+1; p < jA[j+1]; p++) { x[iA[p]] -= vA[p]*x[j]; //  j-  } } return (top) ; } 

نرى أن دورتنا لا تعمل إلا من خلال مؤشرات المجموعة $$$nnz (x)$ ، وليس المجموعة بأكملها $$$0 ، 1 ، ... ، N - 1 دولا$ . انتهى

هناك تطبيق لا يستخدم مجموعة من العلامات markedلحفظ الذاكرة. بدلاً من ذلك ، يتم استخدام مجموعة إضافية من المؤشرات.$$$V_1$ لا تتقاطع مع مجموعة $$$V = \{0, ..., N - 1\}$التي يوجد فيها مناظرة رأس برأس من خلال عملية جبرية بسيطة كإجراء لتعليم العقدة. ومع ذلك ، في عصرنا من الذاكرة الرخيصة ، وحفظه على مجموعة واحدة من الطول$$$N$ يبدو زائدة تماما.

في الختام

تنقسم عملية حل نظام متفرق من المعادلات الخطية بالطريقة المباشرة ، كقاعدة عامة ، إلى ثلاث مراحل:

1. تحليل الشخصية
2. عامل العددية بناء على بيانات تحليل sivol
3. حل النظم الثلاثية التي تم الحصول عليها مع الجانب الأيمن كثيف

الخطوة الثانية ، التخصيم العددي ، هي الجزء الأكثر استهلاكا للموارد وتلتهم معظم (> 90 ٪) من الوقت المقدر. الغرض من الخطوة الأولى هو تقليل التكلفة العالية للخطوة الثانية. تم تقديم مثال للتحليل الرمزي في هذا المنشور. ومع ذلك ، فهذه هي الخطوة الأولى التي تتطلب أطول وقت للتطوير وأقصى تكاليف العقلية من جانب المطور. يتطلب التحليل الرمزي الجيد معرفة نظرية الرسوم والأشجار وامتلاك "غريزة حسابية". والخطوة الثانية أسهل في التنفيذ.
حسنًا ، الخطوة الثالثة ، سواء في التنفيذ أو في الوقت المقدر ، هي في معظم الحالات هي الأكثر بسيطة.

مقدمة جيدة للطرق المباشرة للأنظمة المتناثرة هي في تيم ديفيس ، أستاذ في قسم علوم الكمبيوتر بجامعة تكساس إيه آند إم بعنوان "الطرق المباشرة للأنظمة الخطية المتناثرة ". تم توضيح الخوارزمية الموصوفة أعلاه هناك. لست على دراية لديفيز بنفسي ، رغم أنني كنت أعمل في الجامعة نفسها (وإن كان ذلك في كلية مختلفة). إذا لم أكن مخطئًا ، فقد شارك ديفيس بنفسه في تطوير المذيبات (!). المستخدمة في Matlab وعلاوة على ذلك، كان متورطا حتى بالنسبة الصورة مولدات الكهرباء في جوجل ستريت فيو (شريان الحياة) بالمناسبة، في نفس الكلية أيا المدرجة بخلاف نفسه بيارن ستروستروب - .. خالق C ++

ربما في يوم من الأيام سأكتب شيئًا آخر عن المصفوفات المتفرقة أو الطرق العددية عام. جيد للجميع!