C ++ Mémoire Solution efficace pour Ax = b linéaire Algebra System

https://stackoverflow.com/questions/1242190

12-09-2019
|

Question

J'utilise la bibliothèque numérique bindings pour uBLAS Boost pour résoudre un système linéaire simple. Ci-dessous fonctionne très bien, mais elles sont limitées à la manipulation des matrices A (m x m) de relativement petit 'm'.

Dans la pratique, j'ai une matrice beaucoup plus grande à la dimension m = 10 ^ 6 (jusqu'à 10 ^ 7).
Y at-il approche actuelle de C ++ pour la résolution Ax = b qui utilise efficacement la mémoire.

#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;
}

La solution

Réponse courte: Ne pas utiliser les fixations de LAPACK de Boost, ceux-ci ont été conçus pour des matrices denses, pas des matrices creuses, l'utilisation UMFPACK à la place.

Réponse longue:. UMFPACK est l'une des meilleures bibliothèques pour résoudre Ax = b lorsque A est grande et rare

est un exemple de code ci-dessous (sur la base de umfpack_simple.c) qui génère simple A et b et permet de résoudre 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 fonction crée la matrice A et la droite b. La matrice est d'abord construit sous forme de triplet. le vecteurs de Ti, Tj et Tx décrire complètement A. forme Triplet est facile à créer, mais méthodes de matrices creuses efficaces nécessitent comprimé Sparse format de colonne. Conversion est effectuée avec umfpack_di_triplet_to_col.

A factorisation symbolique est réalisée avec umfpack_di_symbolic. A clairsemée LU décomposition de A est réalisée avec umfpack_di_numeric. Les résout triangulaires inférieurs et supérieurs sont réalisés avec umfpack_di_solve.

Avec n 500.000, sur ma machine, l'ensemble du programme prend environ une seconde pour exécuter. Valgrind rapporte que 369,239,649 octets (un peu plus de 352 Mo) ont été attribués.

Notez cette page discute le soutien de Boost pour les matrices creuses dans triplet (coordonnées) et le format compressé. Si vous le souhaitez, vous pouvez écrire des routines pour convertir ces objets boost aux tableaux simples UMFPACK nécessite en entrée.

Autres conseils

Si l'on suppose vos matrices énormes sont rares, que j'espère qu'ils sont à cette taille, jetez un oeil à la PARDISO projet qui est un solveur linéaire clairsemée, c'est ce que vous aurez besoin si vous voulez gérer les matrices aussi grandes que vous avez dit. Permet un stockage efficace des seules valeurs non nulles, et est beaucoup plus rapide que la résolution du même système de matrices denses.

Je suppose que votre matrice est dense. Si elle est rare, vous pouvez trouver de nombreux algorithmes spécialisés comme déjà mentionné par DeusAduro et duffymo .

Si vous ne disposez pas d'un cluster (assez grand) à votre disposition, vous voulez regarder des algorithmes hors-of-core. ScaLAPACK a quelques solveurs hors du noyau dans le cadre de son paquet prototype , consultez la documentation ici et Google pour plus de détails. Recherche sur le Web pour « out-of-core LU / (matrice) solveurs / packages » vous donnera des liens vers une multitude d'algorithmes et d'outils supplémentaires. Je ne suis pas un expert sur ceux-ci.

Pour ce problème, la plupart des gens utiliseraient un cluster, cependant. Le package que vous trouverez sur presque tous les cluster est ScaLAPACK, encore une fois. De plus, il y a généralement de nombreux autres paquets sur le cluster typique, vous pouvez donc choisir ce qui convient à votre problème (exemples ici et ).

Avant de commencer à coder, vous voulez probablement vérifier rapidement combien de temps il faudra pour résoudre votre problème. Un résolveur typique prend environ O (3 * N ^ 3) bascules (N est la dimension de la matrice). Si N = 100000, vous cherchez donc à 3000000 Gflops. En supposant que votre solveur en mémoire ne 10 Gflops / s par cœur, vous êtes à la recherche à 3 1/2 jours sur un seul noyau. Comme les algorithmes échelle bien, ce qui augmente le nombre de cœurs devrait réduire le temps près de façon linéaire. En plus de cela vient d'E / S.

Je ne sais pas sur les implémentations C +, mais il y a plusieurs choses que vous pouvez faire si la mémoire est un problème en fonction du type de matrice que vous avez affaire à:

Si votre matrice est clairsemée ou bagués, vous pouvez utiliser un solveur ou la bande passante. Ceux-ci ne stockent pas des éléments nuls en dehors du groupe.
Vous pouvez utiliser un solveur de front d'onde, qui stocke la matrice sur le disque et apporte seulement dans la matrice pour la décomposition du front d'onde.
Vous pouvez éviter de résoudre la matrice tout à fait et utiliser des méthodes itératives.
Vous pouvez essayer des méthodes de Monte Carlo de solution.

Jetez un oeil à la liste des logiciels disponibles librement pour la solution des problèmes d'algèbre linéaire , compilé par Jack Dongarra et Hatem Ltaief.

Je pense que pour la taille du problème que vous regardez, vous avez probablement besoin d'un algorithme itératif. Si vous ne voulez pas stocker la matrice A dans un format creux, vous pouvez utiliser une implémentation sans matrice. algorithmes itératifs ne généralement pas besoin d'accéder à des entrées individuelles de la matrice A, ils ne doivent calculer les produits matrice-vecteur Av (et parfois ^ T v, le produit de la matrice transposée avec le vecteur). Donc, si la bibliothèque est bien conçu, il devrait être suffisant si vous passez une classe qui sait comment faire des produits matrice-vecteur.

Comme la réponse acceptée suggère qu'il ya UMFPACK. Mais si vous utilisez BOOST vous pouvez toujours utiliser les matrices compactes en BOOST et utiliser UMFPACK pour résoudre le système. Il y a une liaison qui le rend facile:

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

Ses environ deux ans en date, mais juste un liant (avec quelques autres).

voir la question connexe: UMFPACK et uBLAS Sparse Matrix BOOST

Licencié sous: CC-BY-SA avec attribution

Non affilié à StackOverflow