Ax = b线性代数系统的C++内存有效解决方案

Question

我正在使用用于Boost UBlas的数​​值库绑定来解决简单的线性系统.以下工作正常,但它仅限于处理相对较小的'm'的矩阵A(mxm).

在实践中,我有一个更大的矩阵,维数m = 10 ^ 6(最多10 ^ 7).
是否存在用于解决有效使用内存的Ax = b的现有C++方法.

#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;
}

Answer 1

简短回答:不要使用Boost的LAPACK绑定,这些是针对密集矩阵设计的,而不是稀疏矩阵,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.矩阵首先以三重形式构建.矢量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).

请注意,本页讨论了Boost对Triplet(坐标)和压缩格式的稀疏矩阵的支持.如果您愿意,可以编写例程将这些boost对象转换为简单数组UMFPACK所需的输入.

Answer 2

假设你的巨大矩阵是稀疏的,我希望它们是那么大,看看PARDISO项目是一个稀疏线性求解器,如果你想处理像你说的那么大的矩阵,这就是你需要的.允许仅有效存储非零值,并且比解决相同的密集矩阵系统快得多.

Answer 3

我假设你的矩阵是密集的.如果它很稀疏,你可以找到许多专门的算法,如DeusAduro和duffymo已经提到的那样.

如果您没有可用的(足够大的)群集,则需要查看核心外算法.ScaLAPACK有一些核外解算器作为其原型包的一部分,请参阅此处的文档和Google以获取更多详细信息.在网上搜索"核外LU /(矩阵)求解器/包"将为您提供大量其他算法和工具的链接.我不是那些专家.

但是,对于这个问题,大多数人会使用群集.您将在几乎任何群集上找到的包再次是ScaLAPACK.此外,典型群集上通常还有许多其他软件包,因此您可以选择适合您问题的软件包(此处和此处的示例).

在开始编码之前,您可能希望快速检查解决问题所需的时间.典型的求解器需要大约O(3*N ^ 3)个触发器(N是矩阵的维数).如果N = 100000,那么您将看到3000000 Gflops.假设您的内存中求解器每个核心的速度为10 Gflops/s,那么您在单个核心上看到的是3 1/2天.由于算法可以很好地扩展,因此增加内核数量可以减少接近线性的时间.最重要的是I/O.