Memória C ++ solução eficiente para Ax = b álgebra linear sistema

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

Pergunta

Eu estou usando Numeric Biblioteca Ligações para UBlas impulso para resolver um sistema linear simples. O seguinte funciona muito bem, excepto que se limita a um manuseamento matrizes A (m x m) para relativamente pequena 'm'.

Na prática, têm uma matriz muito maior dimensão com m = 10 ^ 6 (até 10 ^ 7).
Existe existente abordagem C ++ para resolver Ax = b que a memória usos de forma eficiente. 

#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;
}
Foi útil?

Solução

Resposta curta: Não utilize ligações LAPACK de impulso, estes foram projetados para matrizes densas, não matrizes esparsas, o uso UMFPACK vez.

Resposta longa:. UMFPACK é uma das melhores bibliotecas para resolver Ax = b quando A é grande e escassa

Abaixo está o código da amostra (baseado em umfpack_simple.c) que gera um A simples e b e resolve 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);
}

A função generate_sparse_matrix_problem cria a A matriz eo Do lado direito b lado. A matriz é primeiro construído em forma de tripleto. o vetores Ti, Tj, e Tx descrever completamente a forma A. Triplet é fácil de criar, mas eficientes métodos matriciais esparsos exigem formato de coluna Compressed esparsas. Conversão é realizada com umfpack_di_triplet_to_col.

A fatoração simbólica é realizada com umfpack_di_symbolic. A escassa LU decomposição de A é realizada com umfpack_di_numeric. Os resolve inferior e superior triangulares são realizadas com umfpack_di_solve.

Com n como 500.000, na minha máquina, todo o programa leva cerca de um segundo para executar. Valgrind relata que 369,239,649 bytes (um pouco mais de 352 MB) foram alocados.

Nota discute esta página suporte de impulso para matrizes esparsas em Triplet (coordenadas) e formato comprimido. Se quiser, você pode escrever rotinas para converter esses impulsionar objetos às matrizes simples UMFPACK exige como entrada.

Outras dicas

Assumindo que o seu enormes matrizes são esparsas, que eu espero que eles estão nesse tamanho, ter um olhar para o Pardiso projeto que é uma escassa linear solver, é isso que você vai precisar se você quiser lidar com matrizes tão grande como você disse. Permite o armazenamento eficiente de apenas valores diferentes de zero, e é muito mais rápido do que resolver o mesmo sistema de matrizes densas.

Eu supor que seu matriz é densa. Se é escassa, você pode encontrar vários algoritmos especializados como já mencionado por DeusAduro e duffymo .

Se você não tiver um cluster (grande o suficiente) à sua disposição, você quer olhar para algoritmos de out-of-core. ScaLAPACK tem alguns solucionadores de out-of-core como parte de sua protótipo pacote , consulte a documentação aqui e Google para mais detalhes. Pesquisando na web para "out-of-core CN / (matriz) solucionadores / packages" lhe dará links para uma riqueza de novos algoritmos e ferramentas. Eu não sou um especialista sobre aqueles.

Para este problema, a maioria das pessoas iria usar um cluster, no entanto. O pacote que você vai encontrar em quase todo o cluster é ScaLAPACK, novamente. Além disso, geralmente há vários outros pacotes no cluster típico, para que possa escolher o que se adapte ao seu problema (exemplos aqui e aqui ).

Antes de começar a codificar, você provavelmente vai querer verificar rapidamente quanto tempo vai demorar para resolver o seu problema. Um solucionador típico tem cerca de O (3 * N ^ 3) da aleta (N é a dimensão da matriz). Se n = 100000, que são, portanto, olhando 3000000 Gflops. Assumindo que o seu in-memory solver faz 10 Gflops / s por núcleo, você está olhando para 3 1/2 dias em um único núcleo. Como os algoritmos escala bem, aumentando o número de núcleos deve reduzir o tempo perto de linearmente. No topo do que vem a I / O.

Não tenho certeza implementações C ++, mas há várias coisas que você pode fazer se a memória não é um problema, dependendo do tipo de matriz que você está lidando com:

  1. Se a sua matriz é escassa ou em faixas, você pode usar uma escassa ou largura de banda solver. Estes não armazenar zero elementos fora da banda.
  2. Você pode usar um solucionador de frente de onda, que armazena a matriz no disco e só traz na frente de onda de matriz para a decomposição.
  3. Você pode evitar resolver a matriz completamente e usar métodos iterativos.
  4. Você pode tentar métodos de Monte Carlo de solução.

Tenha uma olhada na lista de software disponíveis gratuitamente para a solução da álgebra linear problemas , compilado por Jack Dongarra e Hatem Ltaief.

Eu acho que para o tamanho do problema que você está olhando, você provavelmente precisará de um algoritmo iterativo. Se você não quiser armazenar a matriz A em um formato escasso, você pode usar uma implementação livre-matriz. algoritmos iterativos normalmente não precisam de acesso entradas individuais da matriz A, eles só precisam de produtos de computação matriz-vetor Av (e às vezes um ^ T V, os produtos da matriz transposta com o vector). Então, se a biblioteca está bem concebido, que deve ser suficiente se você passá-lo uma classe que sabe fazer produtos matriz-vetor.

Como a resposta aceita sugere que há UMFPACK. Mas se você estiver usando o impulso que você ainda pode usar as matrizes compactas em BOOST e usar UMFPACK para resolver o sistema. Há uma ligação que torna mais fácil:

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

Seu aproximadamente dois anos datado, mas é apenas uma ligação (juntamente com alguns outros).

ver questão relacionada: UMFPACK e Matrix uBLAS Sparse do IMPULSO

Licenciado em: CC-BY-SA com atribuição
Não afiliado a StackOverflow
scroll top