9

我正在使用 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;
}
4

6 回答 6

14

简短的回答:不要使用 Boost 的LAPACK绑定,这些是为密集矩阵设计的,而不是稀疏矩阵,请UMFPACK改用。

长答案:UMFPACK当 A 大且稀疏时,它是解决 Ax=b 的最佳库之一。

umfpack_simple.c下面是生成简单的示例代码(基于)Ab 解决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需要作为输入的简单数组。

于 2009-08-14T19:30:28.837 回答
6

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

于 2009-08-07T00:09:17.333 回答
6

我假设你的矩阵是密集的。如果它是稀疏的,您可以找到DeusAduroduffymo已经提到的许多专用算法。

如果您没有(足够大的)集群可供您使用,则需要查看核外算法。ScaLAPACK有一些核外求解器作为其原型包的一部分,请参阅此处的文档和Google了解更多详细信息。在网上搜索“核外 LU /(矩阵)求解器 / 包”将为您提供指向更多算法和工具的链接。我不是这些方面的专家。

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

在开始编码之前,您可能想快速检查解决问题需要多长时间。一个典型的求解器大约需要 O(3*N^3) 次翻转(N 是矩阵的维数)。如果 N = 100000,那么您将看到 3000000 Gflops。假设您的内存求解器每个内核执行 10 Gflops/s,那么您在单个内核上查看 3 1/2 天。随着算法的良好扩展,增加内核数量应该会减少接近线性的时间。最重要的是 I/O。

于 2009-08-11T19:07:42.473 回答
3

不确定 C++ 实现,但如果内存是一个问题,取决于您正在处理的矩阵类型,您可以做几件事:

  1. 如果您的矩阵是稀疏或带状的,您可以使用稀疏或带宽求解器。这些不存储带外的零元素。
  2. 您可以使用波前求解器,它将矩阵存储在磁盘上,并且只引入矩阵波前进行分解。
  3. 您可以完全避免求解矩阵并使用迭代方法。
  4. 您可以尝试蒙特卡洛的解决方法。
于 2009-08-07T00:10:54.023 回答
3

查看由 Jack Dongarra 和 Hatem Ltaief 编写的用于解决线性代数问题的免费软件列表。

我认为对于您正在查看的问题大小,您可能需要一个迭代算法。如果您不想以稀疏格式存储矩阵 A,则可以使用无矩阵实现。迭代算法通常不需要访问矩阵 A 的各个条目,它们只需要计算矩阵向量乘积 Av(有时是 A^T v,转置矩阵与向量的乘积)。因此,如果该库设计良好,那么如果您向它传递一个知道如何进行矩阵向量乘积的类就足够了。

于 2009-08-12T09:43:13.793 回答
1

正如公认的答案所暗示的那样,有 UMFPACK。但是如果您使用 BOOST,您仍然可以使用 BOOST 中的紧凑矩阵并使用 UMFPACK 来求解系统。有一个绑定使它变得容易:

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

它大约有两年的历史,但它只是一个绑定(连同其他一些)。

请参阅相关问题: UMFPACK 和 BOOST 的 uBLAS 稀疏矩阵

于 2010-10-22T04:03:37.890 回答