简短的回答:不要使用 Boost 的LAPACK
绑定,这些是为密集矩阵设计的,而不是稀疏矩阵,请UMFPACK
改用。
长答案:UMFPACK
当 A 大且稀疏时,它是解决 Ax=b 的最佳库之一。
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。三元组形式很容易创建,但有效的稀疏矩阵方法需要压缩稀疏列格式。使用 执行转换umfpack_di_triplet_to_col
。
使用 执行符号分解umfpack_di_symbolic
。的稀疏 LU 分解A
是用 执行的umfpack_di_numeric
。使用 执行下三角和上三角求解umfpack_di_solve
。
在n
我的机器上,有 500,000 个,整个程序大约需要一秒钟才能运行。Valgrind 报告说分配了 369,239,649 字节(略多于 352 MB)。
请注意,此页面讨论了 Boost 对 Triplet(坐标)和压缩格式的稀疏矩阵的支持。如果您愿意,您可以编写例程将这些 boost 对象转换为UMFPACK
需要作为输入的简单数组。