C ++ Speicher effiziente Lösung für Ax = b Lineare Algebra System

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

Frage

Ich bin mit Numeric Library Bindings für Boost-UBlas ein einfaches lineares System zu lösen. Die folgende funktioniert gut, außer es auf der Handhabung begrenzte Matrizen A (m x m) für relativ kleine 'm'.

In der Praxis habe ich eine viel größere Matrix mit der Dimension m = 10 ^ 6 (bis zu 10 ^ 7).
Gibt es C ++ Ansatz bestehende Ax = b für die Lösung, die effizient Speicher verwendet. 

#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;
}
War es hilfreich?

Lösung

Kurze Antwort: Sie Boost-die nicht LAPACK Bindungen verwenden, diese für dichte Matrizen entwickelt wurden, nicht dünn besetzte Matrizen, die Verwendung UMFPACK statt.

Lange Antwort:. UMFPACK ist eine der besten Bibliotheken zur Lösung von Ax = b, wenn A ist groß und spärlich

Im Folgenden finden Sie Beispielcode (basierend auf umfpack_simple.c), die eine einfache A und b erzeugt und löst 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);
}

Die Funktion generate_sparse_matrix_problem schafft die Matrix A und die rechten Seite b. Die Matrix wird zuerst in Triplett Form konstruiert. Das Vektoren, Ti, Tj und Tx vollständig A. Triplet Form beschreibt, ist einfach zu erstellen, aber effiziente Sparse Matrix Methoden erfordern Spalte mit geringer Dichte-Format komprimiert. Umwandlung mit umfpack_di_triplet_to_col ausgeführt.

Ein symbolischer Faktorisierung mit umfpack_di_symbolic ausgeführt. Eine spärliche LU-Zerlegung von A mit umfpack_di_numeric durchgeführt. Die untere und die obere Dreieck löst mit umfpack_di_solve durchgeführt werden.

Mit n als 500.000, auf meinem Rechner, nimmt das gesamte Programm um eine zweite zu laufen. Valgrind berichtet, dass 369.239.649 Bytes (nur etwas mehr als 352 MB) zugeteilt wurden.

Beachten Sie diese bespricht Boost-Unterstützung für sparse-Matrizen in Triplet (Koordinaten) und Druckformat. Wenn Sie möchten, können Sie Routinen schreiben, um diese Boost-Objekte zu konvertieren zu dem einfachen Arrays erfordert UMFPACK als Eingabe.

Andere Tipps

Angenommen, Ihre große Matrizen sind spärlich, die ich hoffe, dass sie in dieser Größe sind, haben einen Blick auf die PARDISO Projekt, das eine spärliche lineare Solver ist, ist es das, was Sie brauchen, wenn Sie Matrizen so groß zu handhaben wollen, wie Sie gesagt haben. Ermöglicht effiziente Speicherung von nur Nicht-Null-Werten, und ist viel schneller als das gleiche System von dichten Matrizen zu lösen.

Ich gehe davon aus, dass Ihre Matrix dicht ist. Wenn es spärlich ist, können Sie zahlreiche spezialisierte Algorithmen finden, wie bereits von DeusAduro und duffymo .

Wenn Sie nicht über einen (groß genug) Cluster zur Verfügung haben, möchten Sie auf Out-of-Core-Algorithmen suchen. Scalapack hat ein paar Out-of-Core Solver im Rahmen seines Prototyp-Paket finden Sie in der Dokumentation hier und Google für weitere Details. Sie suchen nach einem „Out-of-Core-LU / (Matrix) Solver / packages“ erhalten Sie Links zu einer Fülle von weiteren Algorithmen und Werkzeuge. Ich bin kein Experte für diese.

Für dieses Problem, würden die meisten Menschen einen Cluster verwenden, jedoch. Das Paket, das Sie fast jeder Cluster finden auf wird ist Scalapack, wieder. Darüber hinaus gibt es in der Regel zahlreiche andere Pakete auf dem typischen Cluster, so können Sie wählen, was passt zu Ihrem Problem (Beispiele hier und hier ).

Vor der Codierung beginnen, möchten Sie wahrscheinlich schnell überprüfen, wie lange es dauern wird, um Ihr Problem zu lösen. Ein typischer Solver dauert etwa O (3 * N ^ 3) Flops (N ist Dimension der Matrix). Wenn N 100000 =, suchen Sie daher bei 3.000.000 GFlops. Unter der Annahme, dass Ihr in-Memory-Löser hat 10 GFlops / s pro Kern, suchen Sie nach 3 1/2 Tagen auf einem einzelnen Kern. Da die Algorithmen gut skalieren, sollte die Anzahl der Kerne linear ansteigend schließt die Zeit reduzieren. Hinzu kommt, dass kommt die I / O.

Nicht sicher C ++ Implementierungen, aber es gibt mehrere Dinge, die Sie tun können, wenn der Speicher ist ein Problem von der Art der Matrix je mit Ihnen zu tun hat:

  1. Wenn Sie Ihre Matrix spärlich oder umreift ist, können Sie einen spärlichen oder Bandbreite Solver verwenden. Diese speichern keine Null-Elemente außerhalb des Bandes.
  2. Sie können einen Wellenfront-Solver verwenden, welche die Matrix auf der Festplatte speichert und bringen nur in der Matrix Wellenfront zur Zersetzung.
  3. Sie vermeiden können die Matrix zusammen und verwenden iterative Methoden zu lösen.
  4. Sie können versuchen, Monte-Carlo-Methoden der Lösung aus.

Haben Sie einen Blick auf die Liste der frei verfügbaren Software für die Lösung der linearen Algebra Probleme , von Jack Dongarra und Hatem Ltaief zusammengestellt.

Ich denke, dass für die Problemgröße an Sie suchen, werden Sie wahrscheinlich einen iterativen Algorithmus benötigen. Wenn Sie die Matrix A in einem spärlichen Format nicht speichern möchten, können Sie eine matrixfreie Implementierung verwenden. Iterative Algorithmen der Regel braucht nicht einzelne Einträge der Matrix A zugreifen zu können, benötigen sie nur Matrix-Vektor-Produkte Av berechnen (und manchmal A ^ T v, das Produkt aus der transponierten Matrix mit dem Vektor). Also, wenn die Bibliothek gut gestaltet, sollte es genug sein, wenn Sie ihm eine Klasse übergeben, die wissen, wie Matrix-Vektor-Produkte zu tun.

Wie die akzeptierte Antwort gibt schon sagt ist UMFPACK. Aber wenn Sie BOOST verwenden, können Sie immer noch die kompakten Matrizen in BOOST verwenden und UMFPACK das System zu lösen. Es gibt eine Bindung, die es macht einfach:

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

Die etwa zwei Jahre alt, aber es ist nur eine Bindung (zusammen mit ein paar anderen).

siehe zugehörige Frage: UMFPACK und uBLAS Sparse Matrix des BOOST

Lizenziert unter: CC-BY-SA mit Zuschreibung
Nicht verbunden mit StackOverflow
scroll top