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

StackOverflow https://stackoverflow.com/questions/1242190

题

我正在使用 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;
}
有帮助吗？

解决方案

简短的回答：不用提升的 LAPACK 绑定，这些都是设计用于密集的矩阵， 不稀疏的矩阵，用 UMFPACK 代替。

只要回答： UMFPACK 是一个最好的图书馆为解决斧=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.稀疏 鲁分解 A 执行与 umfpack_di_numeric.下而上三角解决执行与 umfpack_di_solve.

与 n 作为500,000,我的机器上，整个程序需要有关的第二个运行。才报告说，369,239,649字节(只是略超过352MB)进行分配。

注意这个 页面 讨论提高的支持对于稀疏的矩阵中的三重(协调) 并压制格式。如果你喜欢，你可以编写程序转换这些提高的对象 到简单的阵列 UMFPACK 需要作为输入。

其他提示

假设你的巨大的矩阵是稀疏的，我希望他们在这个大小，看看在 PARDISO 项目，这是一个稀疏线性解算器，这就是，如果你想为你说要处理的矩阵一样大，你需要什么。只允许非零值的有效存储，并且比解决稠密矩阵的同一系统快得多。

我假设你的矩阵是稠密的。如果它很稀疏，您可以找到许多专门的算法，如已经提到的 神之阿杜罗 和 达菲莫.

如果您没有（足够大的）集群可供使用，您需要考虑核外算法。 斯卡拉帕克 有一些核外求解器作为其一部分 原型包, ，参见文档 这里 和 谷歌 更多细节。在网络上搜索“out-of-core LU /（矩阵）求解器/包”将为您提供大量其他算法和工具的链接。我不是这些方面的专家。

然而，对于这个问题，大多数人会使用集群。您几乎可以在任何集群上找到的软件包是 ScaLAPACK。此外，典型集群上通常还有许多其他软件包，因此您可以挑选适合您问题的软件包（示例 这里 和 这里).

在开始编码之前，您可能想快速检查解决问题需要多长时间。典型的求解器大约需要 O(3*N^3) 次触发器（N 是矩阵的维度）。如果 N = 100000，那么您将看到 3000000 Gflops。假设您的内存求解器每个核心的处理速度为 10 Gflops/s，则单个核心的运行时间为 3 1/2 天。由于算法可扩展性良好，增加内核数量应该可以接近线性地减少时间。最重要的是 I/O。

不知道C++实现，但有几件事你可以做的，如果记忆是一个问题根据不同类型的矩阵你处理：

  1. 如果你的矩阵是疏或联，可以使用稀疏或带宽求解。这些不储存的零件以外的频带。
  2. 你可以使用波求解，其中存储的矩阵在磁盘和只带来的矩阵的波前用于分解。
  3. 你可以避免的解决矩阵共和使用的迭代方法。
  4. 你可以试着蒙特卡洛的方法的解决方案。

有一个看看列表线性代数问题时，由Jack唐加拉和哈特姆Ltaief编译。

我认为，问题的规模你看，你可能需要一个迭代算法。如果你不想矩阵A存储在稀疏的格式，你可以使用矩阵自由的实现。迭代算法通常并不需要访问矩阵A的各个条目，他们只需要计算的矩阵矢量乘积AV（有时甲^ T V，用该载体转置矩阵的乘积）。因此，如果库是精心设计的，它应该，如果你通过它知道如何做矩阵向量积一类是不够的。

作为接受的答案提示有UMFPACK。但是，如果你正在使用BOOST你仍然可以使用小型矩阵在BOOST和使用UMFPACK解决系统。有一个结合这使得很容易：

http://mathema.tician.de/software/boost-numeric-bindings

其两年左右过时，但它只是一个（连同其他几个人）结合。

看到相关的问题： UMFPACK和BOOST的uBLAS库稀疏矩阵

许可以下： CC-BY-SA 和 归因
不隶属于 StackOverflow
scroll top