C ++ حل فعال لذاكرة للفأس = B نظام الجبر الخطي

Question

أنا أستخدم ارتباطات مكتبة رقمية لتعزيز UBLAs لحل نظام خطي بسيط. ما يلي يعمل بشكل جيد، إلا أنه يقتصر على التعامل مع المصفوفات A (MXM) صغير نسبيا "م".

في الممارسة العملية، لدي مصفوفة أكبر بكثير مع البعد M = 10 ^ 6 (حتى 10 ^ 7).
هل هناك نهج C ++ الموجود لحل الفأس = B يستخدم ذاكرة بكفاءة.

#include<boost/numeric/ublas/matrix.hpp>
#include<boost/numeric/ublas/io.hpp>
#include<boost/numeric/bindings/traits/ublas_matrix.hpp>
#include<boost/numeric/bindings/lapack/gesv.hpp>
#include <boost/numeric/bindings/traits/ublas_vector2.hpp>

// compileable with this command


//g++ -I/home/foolb/.boost/include/boost-1_38 -I/home/foolb/.boostnumbind/include/boost-numeric-bindings solve_Axb_byhand.cc -o solve_Axb_byhand -llapack


namespace ublas = boost::numeric::ublas;
namespace lapack= boost::numeric::bindings::lapack;


int main()
{
    ublas::matrix<float,ublas::column_major> A(3,3);
    ublas::vector<float> b(3);


    for(unsigned i=0;i < A.size1();i++)
        for(unsigned j =0;j < A.size2();j++)
        {
            std::cout << "enter element "<<i << j << std::endl;
            std::cin >> A(i,j);
        }

    std::cout << A << std::endl;

    b(0) = 21; b(1) = 1; b(2) = 17;

    lapack::gesv(A,b);

    std::cout << b << std::endl;


    return 0;
}

Answer 1

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

اجابة طويلة: UMFPACK هي واحدة من أفضل المكتبات لحل الفأس = B عندما يكون كبيرا ومقاذاا.

• http://www.cise.ufl.edu/research/sparse/umfpack/
• http://www.cise.ufl.edu/research/sparse/umfpack/umfpack/doc/quickstart.pdf.

أدناه هو نموذج التعليمات البرمجية (بناء على umfpack_simple.c) ينشئ بسيطة A و b وحل Ax = b.

#include <stdlib.h>
#include <stdio.h>
#include "umfpack.h"

int    *Ap; 
int    *Ai;
double *Ax; 
double *b; 
double *x; 

/* Generates a sparse matrix problem: 
   A is n x n tridiagonal matrix
   A(i,i-1) = -1;
   A(i,i) = 3; 
   A(i,i+1) = -1; 
*/
void generate_sparse_matrix_problem(int n){
  int i;  /* row index */ 
  int nz; /* nonzero index */
  int nnz = 2 + 3*(n-2) + 2; /* number of nonzeros*/
  int *Ti; /* row indices */ 
  int *Tj; /* col indices */ 
  double *Tx; /* values */ 

  /* Allocate memory for triplet form */
  Ti = malloc(sizeof(int)*nnz);
  Tj = malloc(sizeof(int)*nnz);
  Tx = malloc(sizeof(double)*nnz);

  /* Allocate memory for compressed sparse column form */
  Ap = malloc(sizeof(int)*(n+1));
  Ai = malloc(sizeof(int)*nnz);
  Ax = malloc(sizeof(double)*nnz);

  /* Allocate memory for rhs and solution vector */
  x = malloc(sizeof(double)*n);
  b = malloc(sizeof(double)*n);

  /* Construct the matrix A*/
  nz = 0;
  for (i = 0; i < n; i++){
    if (i > 0){
      Ti[nz] = i;
      Tj[nz] = i-1;
      Tx[nz] = -1;
      nz++;
    }

    Ti[nz] = i;
    Tj[nz] = i;
    Tx[nz] = 3;
    nz++;

    if (i < n-1){
      Ti[nz] = i;
      Tj[nz] = i+1;
      Tx[nz] = -1;
      nz++;
    }
    b[i] = 0;
  }
  b[0] = 21; b[1] = 1; b[2] = 17;
  /* Convert Triplet to Compressed Sparse Column format */
  (void) umfpack_di_triplet_to_col(n,n,nnz,Ti,Tj,Tx,Ap,Ai,Ax,NULL);

  /* free triplet format */ 
  free(Ti); free(Tj); free(Tx);
}


int main (void)
{
    double *null = (double *) NULL ;
    int i, n;
    void *Symbolic, *Numeric ;
    n = 500000;
    generate_sparse_matrix_problem(n);
    (void) umfpack_di_symbolic (n, n, Ap, Ai, Ax, &Symbolic, null, null);
    (void) umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, null, null);
    umfpack_di_free_symbolic (&Symbolic);
    (void) umfpack_di_solve (UMFPACK_A, Ap, Ai, Ax, x, b, Numeric, null, null);
    umfpack_di_free_numeric (&Numeric);
    for (i = 0 ; i < 10 ; i++) printf ("x [%d] = %g\n", i, x [i]);
    free(b); free(x); free(Ax); free(Ai); free(Ap);
    return (0);
}

الوظيفة generate_sparse_matrix_problem يخلق المصفوفة A والجانب الأيمن b. وبعد تم بناء المصفوفة لأول مرة في نموذج ثلاثي. تصف الناقلات TI، TJ، و TX بالكامل A. النموذج الثلاثي السهل إنشاء طرق مصفوفة Sparse الفعالة تتطلب تنسيق عمود متقطع مضغوط. يتم إجراء التحويل مع umfpack_di_triplet_to_col.

يتم إجراء عاملات رمزية مع umfpack_di_symbolic. وبعد التحلل lu متفرق من A يتم تنفيذها مع umfpack_di_numeric. وبعد يتم تنفيذ الحلول الثلاثية الدنيا والعليا umfpack_di_solve.

مع n ك 500،000، على جهازي، يستغرق البرنامج بأكمله حوالي ثانية لتشغيله. تقارير Valgrind أن 369،239،649 بايت (أكثر قليلا من 352 ميغابايت) تم تخصيصها.

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

Answer 2

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

Answer 3

أفترض أن مصفوفةك كثيفة. إذا كان متنازع، فيمكنك العثور على العديد من الخوارزميات المتخصصة كما ذكرت بالفعل Deusaduro. و duffymo..

إذا لم يكن لديك كتلة (كبيرة بما فيه الكفاية) تحت تصرفكم، فأنت تريد أن تنظر إلى الخوارزميات خارج النواة. Scalapack. لديه عدد قليل من الحلول خارج النواة كجزء من حزمة النموذج الأولي, ، انظر الوثائق هنا و غوغل لمزيد من التفاصيل. سيعطيك البحث عن الويب من أجل "Out-Core LU / Matrix) لحزم / حزم" الروابط إلى ثروة من الخوارزميات والأدوات الأخرى. أنا لست خبيرا في تلك.

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

قبل أن تبدأ الترميز، ربما تريد التحقق بسرعة من المدة التي ستستغرقها لحل مشكلتك. يحمل Solver نموذجي حول O (3 * n ^ 3) يتخبط (n هو بعد مصفوفة). إذا كان n = 100000، فأنت في النظر إلى 3000000 Gblops. على افتراض أن Solver الخاص بك في الذاكرة يقوم ب 10 gflops / s لكل كور، فأنت تنظر إلى 3 1/2 أيام في جوهر واحد. نظرا لأن الخوارزميات على نطاق واسع، فإن زيادة عدد النوى يجب أن تقلل من الوقت بالقرب من الخطي. علاوة على ذلك يأتي I / O.

Answer 4

لست متأكدا من تطبيقات C ++، ولكن هناك العديد من الأشياء التي يمكنك القيام بها إذا كانت الذاكرة مشكلة اعتمادا على نوع المصفوفة التي تتعامل معها مع:

1. إذا كانت مصفوفةك متناثرة أو نطفة، فيمكنك استخدام حلال متحرك متناثر أو عرض النطاق الترددي. هذه لا تخزن عناصر صفر خارج الفرقة.
2. يمكنك استخدام حلال موجي موجي، الذي يخزن المصفوفة على القرص ويجلب فقط الموجة الموجية المصفوفة للتحلل.
3. يمكنك تجنب حل المصفوفة تماما واستخدام الأساليب التكرارية.
4. يمكنك تجربة Monte Carlo طرق الحل.

Answer 5

إلقاء نظرة على قائمة البرامج المتاحة بحرية لحل مشاكل الجبر الخطي, ، جمعها جاك دونجارا وحاتم لتليم.

أعتقد أنه بالنسبة لحجم المشكلة التي تبحث فيها، ربما تحتاج إلى خوارزمية تكرارية. إذا كنت لا ترغب في تخزين المصفوفة A في تنسيق متقطع، فيمكنك استخدام تطبيق مجاني للمصفوفات. لا تحتاج الخوارزميات التكرارية عادة إلى الوصول إلى إدخالات فردية من المصفوفة أ، فهي تحتاج فقط إلى حساب منتجات مصفوفة ناقلات AV (وأحيانا a ^ t v، المنتج من المصفوفة المقطوعة مع ناقل). لذلك إذا كانت المكتبة مصممة بشكل جيد، فينبغي أن تكون كافية إذا نقلتها فئة تعرف كيفية القيام بمنتجات المصفوفة المتجهات.

Answer 6

كما تشير الإجابة المقبولة إلى وجود UMFACK. ولكن إذا كنت تستخدم دفعة، فلا يزال بإمكانك استخدام المصفوفات المدمجة في زيادة واستخدام Umfpack لحل النظام. هناك ملزمة مما يجعلها سهلة:

http://mathema.tician.de/software/boost-nummeric-bindings.

في حوالي عامين مؤرخة، ولكنها مجرد ملزمة (جنبا إلى جنب مع عدد قليل من الآخرين).

انظر السؤال الوظيفي:مصفوفة Umfpack و Boost's Ublas Sparse