问题描述
我使用的数字图书馆绑定升压的uBLAS解决一个简单的线性系统。
下面的工作正常，但它仅限于处理矩阵A（米×米）的比较
小M。
在实践中我有维数m = 10 ^ 6（最多10 ^ 7）大得多的矩阵。
是否有现有的C ++算法求解能够有效地使用内存Ax = b的。
＃包括LT＆;升压/数字/ uBLAS库/ matrix.hpp＆GT;
＃包括LT＆;升压/数字/ uBLAS库/ io.hpp＆GT;
＃包括LT＆;升压/数字/绑定/特征/ ublas_matrix.hpp＆GT;
＃包括LT＆;升压/数字/绑定/ LAPACK / gesv.hpp＆GT;
＃包括LT＆;升压/数字/绑定/特征/ ublas_vector2.hpp＆GT;// compileable使用此命令
// G ++ -I /家庭/ foolb / .boost /有/升压1_38 -I /家庭/ foolb / .boostnumbind /有/升压数字绑定solve_Axb_byhand.cc -o solve_Axb_byhand -llapack
命名空间的uBLAS =提高::数字:: uBLAS库;
命名空间LAPACK =提高::数字::绑定:: LAPACK;
诠释的main（）
{
:: uBLAS库矩阵LT;浮动，uBLAS库:: column_major＆GT; A（3,3）;
:: uBLAS库矢量＆lt;浮球GT; B（3）;
对（无符号I = 0; I＆LT; A.size1（）;我++）
对（无符号J = 0; J＆LT; A.size2（）; J ++）
{
性病::法院LT＆;＆LT; 输入元素＆LT;＆LT; I＆LT;＆LT; J＆LT;＆LT;的std :: ENDL;
给std :: cin＆GT;＆GT; A（I，J）;
} 性病::法院LT＆;＆LT; A＆LT;＆LT;的std :: ENDL; B（0）= 21; B（1）= 1; B（2）= 17; LAPACK :: gesv（A，B）; 性病::法院LT＆;＆LT; B＆LT;＆LT;的std :: ENDL;
返回0;
}
简短的回答：不使用Boost的 LAPACK 绑定，这些都是专为稠密矩阵，
不稀疏矩阵，使用 UMFPACK 代替。
龙回答： UMFPACK 是解决Ax = b的时候，A是大型稀疏最好的图书馆之一
下面是样本code（根据 umfpack_simple.c ）生成一个简单的 A 和 b
解决了 Ax = b的。
的#include＆LT;＆stdlib.h中GT;
＃包括LT＆;＆stdio.h中GT;
＃包括umfpack.hINT *鸭;
INT *艾;
双* AX;
双* B;
双* X;/ *生成一个稀疏矩阵的问题：
A是N×N的三角矩阵
A（I，I-1）= -1;
A（I，I）= 3;
A（I，I + 1）= -1;
* /
无效generate_sparse_matrix_problem（INT N）{
INT I; / *行索引* /
INT NZ; / *非零指数* /
诠释NNZ = 2 + 3 *（N-2）+ 2; / *非零数* /
INT *钛; / *行*指数/
INT * TJ; / * COL指数* /
双*的Tx; / * *值/ / *分配内存的三重形式* /
TI =的malloc（sizeof的（INT）* NNZ）;
TJ =的malloc（sizeof的（INT）* NNZ）;
TX =的malloc（sizeof的（双）* NNZ）; / *分配内存为COM pressed稀疏列形式* /
AP = malloc的（的sizeof（int）的*（N + 1））;
艾=的malloc（sizeof的（INT）* NNZ）;
AX =的malloc（sizeof的（双）* NNZ）; / *分配内存RHS及解向量* /
X =的malloc（sizeof的（双）* N）;
B =的malloc（sizeof的（双）* N）; / *构建矩阵A * /
NZ = 0;
对于（i = 0; I＆LT; N;我++）{
如果（I 0）{
钛[NZ] =我;
TJ [NZ] = I-1;
TX [NZ] = -1;
NZ ++;
} 钛[NZ] =我;
TJ [NZ] =我;
TX [NZ] = 3;
NZ ++; 如果（ⅰ＆下;正 - 1）{
钛[NZ] =我;
TJ [NZ] = I + 1;
TX [NZ] = -1;
NZ ++;
}
B〔Ⅰ〕= 0;
}
B [0] = 21; B〔1] = 1; B〔2] = 17;
/ *三重转换对COM pressed稀疏列格式* /
（无效）umfpack_di_triplet_to_col（N，N，NNZ，钛，TJ，TX，鸭，艾，斧，NULL）; / *免费三重格式* /
免费（TI）;免费（TJ）;免费（TX）;
}
INT主要（无效）
{
双* NULL =（双*）NULL;
INT I，N;
无效*符号，数字*;
N = 500000;
generate_sparse_matrix_problem（N）;
（无效）umfpack_di_symbolic（N，N，鸭，艾，斧头，和放大器;象征，NULL，NULL）;
（无效）umfpack_di_numeric（AP，艾，斧，象征性的，和放大器;数字，NULL，NULL）;
umfpack_di_free_symbolic（安培;象征）;
（无效）umfpack_di_solve（UMFPACK_A，鸭，艾，斧，X，B，数字，NULL，NULL）;
umfpack_di_free_numeric（安培;数字）;
对于（I = 0; I＆小于10;我++）printf的（×[％d个] =％G \\ N，I，X [I]）;
免费（B）;免费（X）;免费（AX）;免费（AI）;免费（AP）;
返回（0）;
}
功能 generate_sparse_matrix_problem 创建Matrix A 和
右侧 B 。该矩阵位于三重的形式首先构建。该
向量将钛，TJ和TX全面地描述A.三胞胎的形式很容易创建，但
高效稀疏矩阵方法需要的COM pressed稀疏列格式。转变
与 umfpack_di_triplet_to_col 执行。
一个符号分解与 umfpack_di_symbolic 执行。稀疏 A的LU分解与 umfpack_di_numeric 执行。
下部和上部三角求解的umfpack_di_solve 与执行。
使用 N 50万，我的机器上，整个程序大约需要一秒钟来运行。
Valgrind的报告说，369239649字节（刚刚超过352有点MB）被分配。
请注意这个讨论Boost的对三联稀疏矩阵支持（坐标）
并且通讯pressed格式。如果你喜欢，你可以写程序来提升这些对象转换
对简单数组 UMFPACK 需要输入。
I am using Numeric Library Bindings for Boost UBlas to solve a simple linear system.The following works fine, except it is limited to handling matrices A(m x m) for relativelysmall 'm'.
In practice I have a much larger matrix with dimension m= 10^6 (up to 10^7).
Is there existing C++ approach for solving Ax=b that uses memory efficiently.
#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;
}
Short answer: Don't use Boost's LAPACK bindings, these were designed for dense matrices,not sparse matrices, use UMFPACK instead.
Long answer: UMFPACK is one of the best libraries for solving Ax=b when A is large and sparse.
- http://www.cise.ufl.edu/research/sparse/umfpack/
- http://www.cise.ufl.edu/research/sparse/umfpack/UMFPACK/Doc/QuickStart.pdf
Below is sample code (based on umfpack_simple.c) that generates a simple A and band solves 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);
}
The function generate_sparse_matrix_problem creates the matrix A and theright-hand side b. The matrix is first constructed in triplet form. Thevectors Ti, Tj, and Tx fully describe A. Triplet form is easy to create butefficient sparse matrix methods require Compressed Sparse Column format. Conversionis performed with umfpack_di_triplet_to_col.
A symbolic factorization is performed with umfpack_di_symbolic. A sparseLU decomposition of A is performed with umfpack_di_numeric.The lower and upper triangular solves are performed with umfpack_di_solve.
With n as 500,000, on my machine, the entire program takes about a second to run.Valgrind reports that 369,239,649 bytes (just a little over 352 MB) were allocated.
Note this page discusses Boost's support for sparse matrices in Triplet (Coordinate)and Compressed format. If you like, you can write routines to convert these boost objectsto the simple arrays UMFPACK requires as input.
这篇关于对于Ax = b的线性代数系统C ++内存有效解的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！