Memoria C ++ soluzione efficiente per Ax = b Algebra lineare Sistema

Question

Sto usando numerico Biblioteca Associazioni per uBLAS Boost per risolvere un sistema lineare semplice. Il seguente funziona bene, tranne che viene limitato al trattamento matrici A (m x m) per relativamente piccolo 'm'.

In pratica ho una matrice molto più grande di dimensione m = 10 ^ 6 (fino a 10 ^ 7).
È l'approccio C ++ non esistenti per risolvere Ax = b che utilizza la memoria in modo efficiente.

#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

Risposta breve: non utilizzare attacchi LAPACK di Boost, questi sono stati progettati per le matrici dense, non matrici sparse, uso UMFPACK invece.

Long risposta:. UMFPACK è uno dei migliori librerie per risolvere Ax = b quando A è grande e rada

• http://www.cise.ufl.edu/research/sparse/umfpack /
• http: //www.cise.ufl. edu / ricerca / sparse / umfpack / uMFPACK / doc / QuickStart.pdf

Di seguito è riportato il codice di esempio (sulla base di umfpack_simple.c) che genera un semplice A e b e risolve 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);
}

La generate_sparse_matrix_problem funzione crea l'A matrice e la destro b lato. La matrice viene prima costruito in forma di tripletto. Il vettori di Ti, Tj e Tx descrivere pienamente A. forma Triplet è facile da creare, ma efficaci metodi matrici sparse richiedono formato compresso colonna di tipo sparse. Conversione viene eseguita con umfpack_di_triplet_to_col.

Una fattorizzazione simbolica viene eseguita con umfpack_di_symbolic. A sparse LU decomposizione di A viene eseguita con umfpack_di_numeric. Le risolve triangolari inferiori e superiori sono eseguite con umfpack_di_solve.

Con n come 500.000, sulla mia macchina, l'intero programma dura circa un secondo per l'esecuzione. Valgrind riferisce che 369,239,649 byte (poco più di 352 MB) sono stati assegnati.

Si noti questa pagina href="http://www.boost.org/doc/libs/1_39_0/libs/numeric/ublas/doc/matrix_sparse.htm" rel="noreferrer"> discute Il supporto di Boost per matrici sparse in Triplet (coordinate) e il formato compresso. Se si desidera, è possibile scrivere routine per convertire questi oggetti Boost alle semplici matrici UMFPACK richiede come input.

Answer 2

Assumendo che il enormi matrici sono scarsi, che spero siano in quel formato, date un'occhiata alla Pardiso progetto che è un risolutore lineari sparsi, questo è quello che vi serve se si desidera gestire matrici grandi come lei ha detto. Consente la memorizzazione efficiente dei soli valori diversi da zero, ed è molto più veloce di risolvere lo stesso sistema di matrici dense.

Answer 3

Presumo che la vostra matrice è denso. Se è scarsa, è possibile trovare numerosi algoritmi specializzati come già accennato da DeusAduro e duffymo .

Se non si dispone di una (abbastanza grande) grappolo a vostra disposizione, si consiglia di guardare algoritmi in out-of-core. ScaLAPACK ha un paio di solutori out-of-core come parte della sua pacchetto prototipo , consultare la documentazione qui e Google per maggiori dettagli. La ricerca sul web per "out-of-core LU / (matrice) risolutori / pacchetti" vi darà link ad una vasta gamma di altri algoritmi e strumenti. Io non sono un esperto di quelli.

Per questo problema, la maggior parte delle persone sarebbe utilizzare un cluster, tuttavia. Il pacchetto si trovano su quasi qualsiasi cluster è ScaLAPACK, ancora una volta. Inoltre, di solito ci sono numerosi altri pacchetti del cluster tipico, in modo da poter scegliere quello che si adatta il vostro problema (esempi qui e qui ).

Prima di iniziare la codifica, probabilmente si desidera controllare rapidamente quanto tempo ci vorrà per risolvere il problema. Un tipico risolutore dura circa O (3 * N ^ 3) flop (N è la dimensione della matrice). Se N = 100000, si sta quindi guardando 3000000 Gflops. Supponendo che il risolutore in memoria fa 10 GFLOPS / s per core, siete alla ricerca di 3 1/2 giorni su un singolo core. Come gli algoritmi scala bene, aumentando il numero di nuclei dovrebbe ridurre il tempo vicino al linearmente. In cima a che arriva il I / O.

Answer 4

Non sei sicuro di C ++ implementazioni, ma ci sono diverse cose che puoi fare se la memoria è un problema a seconda del tipo di matrice hai a che fare con:

1. Se la matrice è sparsa o fasciato, è possibile utilizzare un risolutore rada o di larghezza di banda. Questi non memorizzano zero elementi al di fuori della band.
2. È possibile utilizzare un risolutore di fronte d'onda, che memorizza la matrice su disco e solo porta nel fronte d'onda matrice per la decomposizione.
3. È possibile evitare di risolvere la matrice del tutto e utilizzare metodi iterativi.
4. Potete provare metodi Monte Carlo di soluzione.

Answer 5

Date un'occhiata alla elenco di software liberamente disponibile per la soluzione di algebra lineare problemi , compilati da Jack Dongarra e Hatem Ltaief.

Credo che per la dimensione del problema che stai guardando, probabilmente avete bisogno di un algoritmo iterativo. Se non si desidera memorizzare la matrice A in un formato sparse, è possibile utilizzare un'implementazione libera matrice. algoritmi iterativi tipicamente non devono accedere iscrizioni individuali della matrice A, hanno solo bisogno di calcolare matrice-vettore prodotti Av (e talvolta A ^ T v, il prodotto della matrice trasposta con il vettore). Quindi, se la libreria è ben progettato, dovrebbe essere sufficiente se si passa una classe che sa fare prodotti matrice-vettore.

Answer 6

Come suggerisce la risposta accettata c'è UMFPACK. Ma se si utilizza BOOST è comunque possibile utilizzare le matrici compatte in BOOST e utilizzare UMFPACK per risolvere il sistema. C'è un legame che lo rende facile:

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

I suoi circa due anni vecchie, ma il suo solo un legame (insieme a pochi altri).

vedi domanda relativa: UMFPACK e di BOOST uBLAS Sparse Matrix