2009-08-07 27 views
8

मैं एक सरल रैखिक प्रणाली को हल करने के लिए बूस्ट उब्लास के लिए न्यूमेरिक लाइब्रेरी बाइंडिंग का उपयोग कर रहा हूं। निम्नलिखित ठीक काम करता है, सिवाय इसके कि यह अपेक्षाकृत छोटे 'एम' के लिए मैट्रिस ए (एम एक्स एम) को संभालने तक ही सीमित है।सी ++ एक्सएक्स के लिए मेमोरी कुशल समाधान = बी रैखिक बीजगणित प्रणाली

प्रैक्टिस में मेरे पास आयाम एम = 10^6 (10^7 तक) के साथ एक बड़ा मैट्रिक्स है।
क्या एक्स = बी को हल करने के लिए मौजूदा सी ++ दृष्टिकोण है जो मेमोरी का कुशलतापूर्वक उपयोग करता है।

#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; 
} 
+2

स्पष्ट बाहर इशारा करते हुए एक मैट्रिक्स है कि आकार सरणी 4x10^12 4x10^14 बाइट्स, या एक ही मैट्रिक्स अकेले के लिए 4 से 400 टेराबाइट्स है:

संबंधित सवाल देखें। (जब तक, जैसा कि नीचे बताया गया है, इसका स्पैस) – cyberconte

उत्तर

13

लघु जवाब: बूस्ट के LAPACK बाइंडिंग का उपयोग न करें, ये घने मैट्रिक्स के लिए डिजाइन किए गए थे, विरल नहीं मेट्रिसेस, UMFPACK बजाय का उपयोग करें।

लंबे उत्तर: UMFPACK ए = = बी को हल करने के लिए सबसे अच्छी पुस्तकालयों में से एक है जब ए बड़ा और विचित्र है।

नीचे नमूना कोड (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 पैदा करता है। मैट्रिक्स पहले तीन गुना रूप में बनाया गया है। वेक्टर टीआई, टीजे, और टीएक्स पूरी तरह से वर्णन करते हैं। ए ट्रिपलेट फॉर्म बनाना आसान है लेकिन कुशल स्पैर मैट्रिक्स विधियों को संपीड़ित स्पैस कॉलम प्रारूप की आवश्यकता होती है। रूपांतरण umfpack_di_triplet_to_col के साथ किया जाता है।

umfpack_di_symbolic के साथ एक प्रतीकात्मक कारक किया जाता है। एक स्पैस A का LU अपघटन umfpack_di_numeric के साथ किया जाता है। निचले और ऊपरी त्रिकोणीय हल umfpack_di_solve के साथ किए जाते हैं।

n के साथ 500,000 के रूप में, मेरी मशीन पर, पूरे कार्यक्रम को चलाने के लिए लगभग एक सेकंड लगता है। वालग्रिंड रिपोर्ट करता है कि 36 9, 3 9, 64 9 बाइट्स (केवल 352 एमबी से थोड़ा अधिक) आवंटित किए गए थे।

नोट करें page ट्रिपलेट (समन्वय) और संपीड़ित प्रारूप में स्पैर मैट्रिस के लिए बूस्ट के समर्थन पर चर्चा करता है। यदि आप चाहें, तो आप इन बूस्ट ऑब्जेक्ट्स को सरल सरणी UMFPACK में इनपुट के रूप में बदलने के लिए दिनचर्या लिख ​​सकते हैं।

+0

+1 :) – ccook

6

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

+2

निष्पक्ष समाधान की ओ (एम^3) समय जटिलता का उल्लेख नहीं है! यहां तक ​​कि चालाक एक नथ वार्ता ओ (एम^2.7ish) है ... यदि ये मैट्रिक्स स्पैस नहीं हैं तो आपको क्लस्टर और प्रथम श्रेणी संख्यात्मक विश्लेषक की आवश्यकता है ...स्पैर मैट्रिक्स विचार के लिए – dmckee

+1

+1। मुझे विभिन्न स्पैस मैट्रिक्स लाइब्रेरी ftp://ftp.numerical.rl.ac.uk/pub/reports/ghsRAL200505.pdf की तुलना करने के बारे में पार्डिसो पेपर में न्यूमेरस लाइब्रेरी और थ्रू तुलना मिली, इसका उपयोग अन्य मान्यता प्राप्त स्पैर मैट्रिक्स पुस्तकालयों को खोजने के लिए किया जा सकता है। –

3
सी ++ कार्यान्वयन के बारे में नहीं

यकीन है, लेकिन वहाँ कई चीजें अगर स्मृति मैट्रिक्स के प्रकार पर निर्भर एक मुद्दा आप साथ काम कर रहे है आप कर सकते हैं: अपने मैट्रिक्स विरल या बैंडेड, आप है

  1. हैं एक स्पैस या बैंडविड्थ सॉल्वर का उपयोग कर सकते हैं। ये बैंड के बाहर शून्य तत्वों को स्टोर नहीं करते हैं।
  2. आप एक वेवफ़्रंट सॉल्वर का उपयोग कर सकते हैं, जो डिस्क पर मैट्रिक्स संग्रहीत करता है और केवल अपघटन के लिए मैट्रिक्स वेवफ़्रंट लाता है।
  3. आप मैट्रिक्स को पूरी तरह से हल करने और पुनरावृत्तियों के तरीकों का उपयोग करने से बच सकते हैं।
  4. आप समाधान के मोंटे कार्लो विधियों का प्रयास कर सकते हैं।
+0

@ डफिमो: धन्यवाद। मैंने सी ++ में पुनरावर्तक दृष्टिकोण कार्यान्वयन को देखा है लेकिन उन्हें अभी भी इसे मैट्रिक्स में संग्रहीत करने की आवश्यकता है। http://freenet-homepage.de/guwi17/ublas/examples/ यदि मैं गलत हूं, तो क्या आप सी ++ में पुनरावृत्ति के लिए कोई भी कुशल कुशल कार्यान्वयन जानते हैं? – neversaint

+0

सही, मूर्खतापूर्ण। मुझे याद रखना चाहिए था। मैं समांतर एल्गोरिदम की जांच करता हूं, क्योंकि एन प्रोसेसर को काम को विभाजित करने और परिणाम प्राप्त करने के लिए इसे एक साथ बुनाई करने की समस्या जर्मनी को अस्थायी रूप से डिस्क पर ले जाने की समस्या के लिए जर्मन है। – duffymo

6

मुझे लगता है कि आपका मैट्रिक्स घना है। यदि यह स्पैस है, तो आप पहले से ही DeusAduro और duffymo द्वारा वर्णित कई विशेष एल्गोरिदम पा सकते हैं।

यदि आपके पास अपने निपटारे में एक (पर्याप्त पर्याप्त) क्लस्टर नहीं है, तो आप आउट-ऑफ-कोर एल्गोरिदम देखना चाहते हैं। ScaLAPACK के पास prototype package के हिस्से के रूप में कुछ आउट-ऑफ-कोर सॉल्वर हैं, अधिक जानकारी के लिए प्रलेखन here और Google देखें। "आउट-ऑफ-कोर LU/(मैट्रिक्स) सॉल्वर/पैकेज" के लिए वेब पर खोज करने से आपको आगे के एल्गोरिदम और टूल की संपत्ति मिल जाएगी। मैं उन पर एक विशेषज्ञ नहीं हूँ।

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

कोडिंग शुरू करने से पहले, शायद आप जल्दी से जांचना चाहते हैं कि आपकी समस्या का समाधान करने में कितना समय लगेगा। एक ठेठ सॉल्वर ओ (3 * एन^3) फ्लॉप (एन मैट्रिक्स का आयाम) लेता है। यदि एन = 100000, तो आप 3000000 Gflops को देख रहे हैं। यह मानते हुए कि आपका इन-मेमोरी सॉल्वर प्रति ग्राफ़ 10 Gflops/s करता है, आप एक कोर पर 3 1/2 दिन देख रहे हैं। जैसे-जैसे एल्गोरिदम अच्छी तरह से स्केल करते हैं, कोर की संख्या में वृद्धि को रैखिक रूप से करीब समय कम करना चाहिए। उस पर सबसे ऊपर आई/ओ आता है।

+0

चेतावनी: उपर्युक्त ओ (3 * एन^3) मानता है कि आप जटिल संख्याओं का उपयोग करते हैं। वास्तविक संख्याओं के लिए, सब कुछ 6 से विभाजित करें, यानी ओ (0.5 * एन^3) के आसपास कहीं। स्कूल गर्व के लिए – stephan

3

जैक डोंगाररा और हैमत लताफ द्वारा संकलित list of freely available software for the solution of linear algebra problems पर एक नज़र डालें।

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

1

स्वीकृत उत्तर के अनुसार यूएमएफपीएक्स है। लेकिन अगर आप बूस्ट का उपयोग कर रहे हैं तो आप अभी भी बूस्ट में कॉम्पैक्ट मैट्रिस का उपयोग कर सकते हैं और सिस्टम को हल करने के लिए यूएमएफपीएक्स का उपयोग कर सकते हैं। वहाँ एक बाध्यकारी जो इसे आसान बना देता है:

दो के बारे में

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

इसके साल दिनांकित है, लेकिन इसकी सिर्फ एक (कुछ अन्य लोगों के साथ) बंधन। UMFPACK and BOOST's uBLAS Sparse Matrix