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 A và b 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.
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