2009-08-07 38 views
8

Tôi đang sử dụng Bindings thư viện số để tăng cường UBlas để giải quyết một hệ thống tuyến tính đơn giản. Các công trình sau đây tốt, ngoại trừ nó được giới hạn để xử lý ma trận A (m x m) cho tương đối nhỏ 'm'.Giải pháp hiệu quả bộ nhớ C++ cho Ax = b Hệ thống đại số tuyến tính

Trong thực tế, tôi có ma trận lớn hơn nhiều với kích thước m = 10^6 (tối đa 10^7).
Có cách tiếp cận C++ hiện có để giải quyết Ax = b sử dụng bộ nhớ hiệu quả không.

#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

Chỉ ra rõ ràng, một ma trận mà mảng kích thước là 4x10^12 đến 4x10^14 byte, hoặc 4-400 terabyte cho một ma trận đơn một mình. (trừ khi, như đã nêu dưới đây, thưa thớt) – cyberconte

Trả lời

13

Câu trả lời ngắn: Không sử dụng các liên kết LAPACK của Boost, chúng được thiết kế cho các ma trận dày đặc, không phải là ma trận thưa thớt, thay vào đó hãy sử dụng UMFPACK.

Câu trả lời dài: UMFPACK là một trong những thư viện tốt nhất để giải quyết Ax = b khi A lớn và thưa thớt.

Dưới đây là mẫu mã (dựa trên umfpack_simple.c) mà tạo ra một đơn giản Ab và giải quyết 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); 
} 

Chức năng generate_sparse_matrix_problem tạo ra ma trận A và phía bên phải b. Ma trận được xây dựng lần đầu ở dạng ba. Các vectơ Ti, Tj và Tx mô tả đầy đủ A. Biểu mẫu Triplet rất dễ tạo nhưng các phương thức ma trận thưa thớt hiệu quả yêu cầu định dạng Cột thưa thớt nén. Chuyển đổi được thực hiện với umfpack_di_triplet_to_col.

Một yếu tố tượng trưng được thực hiện với umfpack_di_symbolic. Phân đoạn thưa thớt LU của A được thực hiện với umfpack_di_numeric. Giải quyết hình tam giác dưới và trên được thực hiện với umfpack_di_solve.

Với n là 500.000, trên máy của tôi, toàn bộ chương trình mất khoảng một giây để chạy. Valgrind báo cáo rằng 369,239,649 byte (chỉ hơn một chút 352 MB) đã được phân bổ.

Lưu ý điều này page thảo luận về hỗ trợ của Boost cho ma trận thưa thớt trong Triplet (Tọa độ) và định dạng Nén. Nếu bạn thích, bạn có thể viết thường trình để chuyển đổi các đối tượng tăng này thành các mảng đơn giản UMFPACK yêu cầu làm đầu vào.

+0

+1 cho niềm tự hào của trường :) – ccook

6

Giả sử ma trận khổng lồ của bạn là thưa thớt, mà tôi hy vọng họ đang ở kích thước đó, có một cái nhìn tại dự án PARDISO mà là một giải tuyến tính thưa thớt, đây là những gì bạn cần nếu bạn muốn xử lý ma trận lớn như bạn đã nói. Cho phép lưu trữ hiệu quả các giá trị khác không, và nhanh hơn nhiều so với việc giải quyết cùng một hệ thống các ma trận dày đặc.

+2

Chưa kể đến độ phức tạp thời gian O (m^3) của giải pháp ngây thơ! Ngay cả một trong những thông minh Knuth nói về là O (m^2.7ish) ... Nếu những ma trận không thưa thớt bạn cần một cụm và một phân tích số lớp học đầu tiên ... – dmckee

+1

+1 cho ý tưởng ma trận thưa thớt. Tôi đã tìm thấy các thư viện số và so sánh có trong bài báo PARDISO về việc so sánh các thư viện ma trận thưa thớt khác nhau ftp://ftp.numerical.rl.ac.uk/pub/reports/ghsRAL200505.pdf Điều này có thể được sử dụng để tìm các thư viện ma trận thưa thớt được công nhận khác. –

3

Không chắc về triển khai C++, nhưng có một số điều bạn có thể làm gì nếu bộ nhớ là một vấn đề tùy thuộc vào loại ma trận bạn đang làm việc với:

  1. Nếu ma trận của bạn là thưa thớt hoặc dải, bạn có thể sử dụng trình giải quyết băng thông hoặc thưa thớt. Chúng không lưu trữ các phần tử zero bên ngoài ban nhạc.
  2. Bạn có thể sử dụng bộ giải mã sóng, lưu trữ ma trận trên đĩa và chỉ đưa vào mặt trận ma trận để phân hủy.
  3. Bạn có thể tránh giải quyết ma trận hoàn toàn và sử dụng các phương pháp lặp.
  4. Bạn có thể thử phương pháp giải pháp Monte Carlo.
+0

@duffymo: cảm ơn. Tôi đã xem xét thực hiện tiếp cận lặp trong C + + nhưng họ vẫn yêu cầu lưu trữ nó trong ma trận. http://freenet-homepage.de/guwi17/ublas/examples/ Nếu tôi sai, Bạn có biết bất kỳ triển khai thực hiện mem nào hiệu quả trong C++ cho lặp lại không? – neversaint

+0

Đúng, ngu xuẩn. Tôi nên nhớ điều đó. Tôi muốn điều tra các thuật toán song song, bởi vì vấn đề phân vùng công việc cho N bộ vi xử lý và đan nó lại với nhau để có được kết quả là germane cho vấn đề di chuyển nó tạm thời ra đĩa. – duffymo

6

Tôi cho rằng ma trận của bạn dày đặc. Nếu nó thưa thớt, bạn có thể tìm thấy nhiều thuật toán chuyên ngành như đã được đề cập bởi DeusAduroduffymo.

Nếu bạn không có cụm (đủ lớn) theo ý của mình, bạn muốn xem các thuật toán ngoài lõi. ScaLAPACK có một vài bộ giải mã ngoài lõi như một phần của prototype package, xem tài liệu hereGoogle để biết thêm chi tiết. Tìm kiếm trên web cho "các gói/gói giải quyết LU/(ma trận ngoài lõi)" sẽ cung cấp cho bạn các liên kết đến vô số các thuật toán và công cụ khác. Tôi không phải là một chuyên gia về chúng.

Đối với vấn đề này, hầu hết mọi người sẽ sử dụng một cụm, tuy nhiên. Các gói phần mềm bạn sẽ tìm thấy trên hầu như bất kỳ cụm là ScaLAPACK, một lần nữa. Ngoài ra, thường có nhiều gói khác trên cụm điển hình, vì vậy bạn có thể chọn và chọn những gì phù hợp với vấn đề của mình (ví dụ: herehere).

Trước khi bắt đầu mã hóa, bạn có thể muốn nhanh chóng kiểm tra xem sẽ mất bao lâu để giải quyết vấn đề của mình. Một người giải quyết điển hình mất khoảng O (3 * N^3) thất bại (N là kích thước của ma trận). Nếu N = 100000, bạn đang xem xét 3000000 Gflops. Giả sử rằng bộ giải mã trong bộ nhớ của bạn có 10 Gflops/s trên mỗi lõi, bạn đang xem 3 1/2 ngày trên một lõi đơn. Khi các thuật toán mở rộng, việc tăng số lượng lõi sẽ giảm thời gian gần với tuyến tính. Trên hết là I/O.

+0

Lưu ý: O (3 * N^3) ở trên giả định rằng bạn sử dụng các số phức. Đối với các số thực, chia mọi thứ cho 6, tức là ở đâu đó xung quanh O (0,5 * N^3). – stephan

3

Hãy xem list of freely available software for the solution of linear algebra problems, do Jack Dongarra và Hatem Ltaief biên soạn.

Tôi nghĩ rằng đối với kích thước vấn đề bạn đang xem, bạn có thể cần một thuật toán lặp lại. Nếu bạn không muốn lưu trữ ma trận A ở định dạng thưa thớt, bạn có thể sử dụng triển khai ma trận miễn phí.Thuật toán lặp lại thường không cần truy cập các mục riêng lẻ của ma trận A, chúng chỉ cần tính toán các sản phẩm vector ma trận Av (và đôi khi A^T v, sản phẩm của ma trận transposed với vec-tơ). Vì vậy, nếu thư viện được thiết kế tốt, nó sẽ là đủ nếu bạn vượt qua nó một lớp học biết làm thế nào để làm các sản phẩm ma trận-vector.

1

Khi câu trả lời được chấp nhận cho thấy có UMFPACK. Nhưng nếu bạn đang sử dụng BOOST, bạn vẫn có thể sử dụng các ma trận nhỏ gọn trong BOOST và sử dụng UMFPACK để giải quyết hệ thống. Có một ràng buộc mà làm cho nó dễ dàng:

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

của nó khoảng hai năm hẹn hò, nhưng nó chỉ là một ràng buộc (cùng với một vài người khác).

xem câu hỏi liên quan: UMFPACK and BOOST's uBLAS Sparse Matrix

Các vấn đề liên quan