我正在使用 Boost UBlas 的数值库绑定(bind)来求解一个简单的线性系统。 以下工作正常,除了它仅限于处理矩阵 A(m x m) 相对 小“m”。
在实践中,我有一个更大的矩阵,维度 m= 10^6(最多 10^7)。
是否存在有效使用内存的现有 C++ 方法来解决 Ax=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;
}
最佳答案
简短回答:不要使用 Boost 的 LAPACK
绑定(bind),这些绑定(bind)是为密集矩阵设计的,
不是稀疏矩阵,请改用 UMFPACK
。
长答案:UMFPACK
是解决 A 大且稀疏时 Ax=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
。该矩阵首先以三元组形式构建。这
vector Ti、Tj 和 Tx 完全描述了 A。三元组形式很容易创建,但
高效的稀疏矩阵方法需要压缩稀疏列格式。转换
使用 umfpack_di_triplet_to_col
执行。
使用 umfpack_di_symbolic
执行符号分解。稀疏的
A
的 LU 分解是使用 umfpack_di_numeric
执行的。
使用 umfpack_di_solve
执行下三角和上三角求解。
n
为 500,000,在我的机器上,整个程序运行大约需要一秒钟。
Valgrind 报告分配了 369,239,649 字节(略高于 352 MB)。
注意这个page讨论 Boost 对 Triplet 中稀疏矩阵的支持(坐标)
和压缩格式。如果你愿意,你可以编写例程来转换这些 boost 对象
到简单数组 UMFPACK
需要作为输入。
关于Ax=b 线性代数系统的 C++ 内存高效解决方案,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/1242190/