Ax = B 선형 대수 시스템에 대한 C ++ 메모리 효율적인 솔루션

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

문제

나는 단순한 선형 시스템을 해결하기 위해 부스트 UBLA를 위해 숫자 라이브러리 바인딩을 사용하고 있습니다. 다음은 비교적 작은 '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 A가 크고 희소 할 때 AX = 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. 희소 한 lu 분해 A 수행됩니다 umfpack_di_numeric. 하부 및 상부 삼각형 해결은 함께 수행됩니다 umfpack_di_solve.

와 함께 n 내 컴퓨터에서 500,000으로 전체 프로그램을 실행하는 데 약 1 초가 걸립니다. Valgrind는 369,239,649 바이트 (352MB 이상)가 할당되었다고보고했다.

이것에 주목하십시오 페이지 트리플렛 (좌표) 및 압축 형식의 희소 행렬에 대한 Boost의 지원에 대해 논의합니다. 원하는 경우 루틴을 작성하여 이러한 부스트 객체를 간단한 배열로 변환 할 수 있습니다. UMFPACK 입력이 필요합니다.

다른 팁

당신의 거대한 행렬이 드물다고 가정하면, 나는 그들이 그 크기에 있기를 바랍니다. Pardiso 희소 한 선형 솔버 인 프로젝트, 이것은 당신이 말한 것처럼 매트릭스를 처리하려면 필요한 것입니다. 0이 아닌 값 만 효율적으로 저장할 수 있으며 동일한 조밀 한 매트릭스 시스템을 해결하는 것보다 훨씬 빠릅니다.

나는 당신의 매트릭스가 조밀하다고 가정합니다. 드문 경우, 이미 언급 한 수많은 특수 알고리즘을 찾을 수 있습니다. Deusaduro 그리고 더 피니 모.

처분 할 때 (충분히 큰) 클러스터가없는 경우, 코어 외 알고리즘을보고 싶습니다. Scalapack 그 일부로 코어 외곽 솔버가 몇 개 있습니다 프로토 타입 패키지, 문서를 참조하십시오 여기 그리고 Google 자세한 사항은. "코어 외부 LU / (매트릭스) 솔버 / 패키지"에 대한 웹을 검색하면 풍부한 추가 알고리즘 및 도구에 대한 링크가 제공됩니다. 나는 그것에 대한 전문가가 아닙니다.

그러나이 문제의 경우 대부분의 사람들은 클러스터를 사용합니다. 거의 모든 클러스터에서 찾을 수있는 패키지는 Scalapack입니다. 또한 일반적인 클러스터에는 일반적으로 많은 다른 패키지가 있으므로 문제에 맞는 것을 선택하고 선택할 수 있습니다 (예제. 여기 그리고 여기).

코딩을 시작하기 전에 문제를 해결하는 데 걸리는 시간을 빨리 확인하고 싶을 것입니다. 전형적인 솔버는 약 O (3*n^3) 플롭 (N is dimension of matrix)을 취합니다. n = 100000이라면 3000000 gflops를보고 있습니다. 메모리 인 솔버가 코어 당 10 GFLOPS/S를 수행한다고 가정하면 단일 코어에서 3 1/2 일을보고 있습니다. 알고리즘이 잘 균형을 유지함에 따라 코어 수를 늘리면 선형에 가까운 시간을 줄여야합니다. 그 외에도 I/O가 온다.

C ++ 구현에 대해서는 확실하지 않지만 처리중인 매트릭스 유형에 따라 메모리가 문제가되는 경우 몇 가지 할 수있는 일이 있습니다.

  1. 매트릭스가 드문 또는 밴드 인 경우 희소 또는 대역폭 솔버를 사용할 수 있습니다. 이들은 밴드 밖에서 제로 요소를 저장하지 않습니다.
  2. 디스크에 매트릭스를 저장하고 분해를 위해 매트릭스 파면에만 가져 오는 웨이브 프론트 솔버를 사용할 수 있습니다.
  3. 매트릭스를 모두 해결하지 않고 반복적 인 방법을 사용할 수 있습니다.
  4. 몬테 카를로 솔루션 방법을 사용해 볼 수 있습니다.

살펴보십시오 선형 대수 문제 솔루션을위한 자유롭게 사용할 수있는 소프트웨어 목록, Jack Dongarra와 Hatem ltaief가 편집했습니다.

보고있는 문제 크기의 경우 반복 알고리즘이 필요할 것입니다. 매트릭스 A를 드문 형식으로 저장하지 않으려면 행렬없는 구현을 사용할 수 있습니다. 반복 알고리즘은 일반적으로 매트릭스 A의 개별 항목에 액세스 할 필요가 없으며, 매트릭스 벡터 제품 AV (때로는 벡터와의 전치 된 매트릭스의 생성물 인 A^t V) 만 계산하면됩니다. 따라서 라이브러리가 잘 설계된 경우 매트릭스 벡터 제품을 수행하는 방법을 알고있는 클래스를 통과하면 충분해야합니다.

허용 된 답변에서 알 수 있듯이 Umfpack이 있습니다. 그러나 부스트를 사용하는 경우 여전히 Compact Matrices를 부스트에서 사용하고 Umfpack을 사용하여 시스템을 해결할 수 있습니다. 쉽게 만드는 바인딩이 있습니다.

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

약 2 년이 지났지 만 그냥 구속력이 있습니다 (다른 사람들과 함께).

관련 질문 참조 :Umfpack 및 Boost의 Ublas 스파 스 매트릭스

라이센스 : CC-BY-SA ~와 함께 속성
제휴하지 않습니다 StackOverflow
scroll top