C ++ Memoria solución eficiente para Ax = b Sistema de Álgebra Lineal

https://stackoverflow.com/questions/1242190

12-09-2019
|

Pregunta

Estoy utilizando numérico Biblioteca enlaces para UBlas Boost para resolver un sistema lineal simple. La siguiente funciona bien, excepto que se limita a la manipulación de matrices A (m m x) para relativamente pequeño 'm'.

En la práctica tengo una matriz mucho más grande con dimensión m = 10 ^ 6 (hasta 10 ^ 7).
¿Hay enfoque C ++ existentes para la solución de Ax = b que utiliza memoria de manera 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;
}

Solución

Respuesta corta: No utilice fijaciones LAPACK de Boost, estos fueron diseñados para matrices densas, no matrices dispersas, el uso UMFPACK su lugar.

Respuesta larga:. UMFPACK es una de las mejores librerías para la solución de Ax = b cuando A es grande y escaso

A continuación se muestra el código de ejemplo (basado en umfpack_simple.c) que genera una A simple y b y resuelve 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);
}

El generate_sparse_matrix_problem función crea el A matriz y la b lado derecho. La matriz se construye primero en forma de triplete. los vectores de Ti, Tj, y Tx describe completamente A. forma Triplet es fácil crear pero métodos eficientes de matriz dispersa requieren formato de columna dispersa comprimido. Conversión se realiza con umfpack_di_triplet_to_col.

A factorización simbólica se realiza con umfpack_di_symbolic. Una escasa Descomposición LU de A se realiza con umfpack_di_numeric. Los resuelve triangulares inferiores y superiores se realizan con umfpack_di_solve.

Con n 500.000, en mi máquina, todo el programa tarda aproximadamente un segundo para correr. Valgrind informa que se asignaron 369,239,649 bytes (sólo un poco más de 352 MB).

Tenga en cuenta esta página href="http://www.boost.org/doc/libs/1_39_0/libs/numeric/ublas/doc/matrix_sparse.htm" discute El apoyo de impulso para matrices dispersas en trío (Coordinate) y un formato comprimido. Si lo desea, puede escribir rutinas para convertir estos objetos a impulsar a las matrices simples UMFPACK requiere como entrada.

Otros consejos

Si se asume sus enormes matrices son escasos, lo que espero que sean en ese tamaño, echar un vistazo a la Pardiso proyecto que es un solucionador lineal dispersa, esto es lo que necesita si desea gestionar las matrices tan grandes como usted ha dicho. Permite el almacenamiento eficiente de los únicos valores que no son cero, y es mucho más rápido que la solución del mismo sistema de matrices densas.

Asumo que su matriz es densa. Si es escasa, se pueden encontrar numerosos algoritmos especializados como ya se ha mencionado por DeusAduro y duffymo .

Si no se dispone de un clúster (bastante grande) a su disposición, que desee ver en los algoritmos de fuera del núcleo. ScaLAPACK tiene algunos solucionadores de fuera del núcleo como parte de su paquete prototipo , consulte la documentación de aquí y Google para más detalles. Buscando en la web "fuera de núcleo LU / (matriz) solucionadores / paquetes" le dará enlaces a una gran cantidad de nuevos algoritmos y herramientas. No soy un experto en ellos.

Para este problema, la mayoría de la gente utiliza un clúster, sin embargo. El paquete encontrará en casi cualquier grupo es ScaLAPACK, otra vez. Además, por lo general hay muchos otros paquetes en el clúster típico, por lo que pueden escoger y elegir lo que se adapte a su problema (ejemplos aquí y aquí ).

Antes de empezar a programar, es probable que desee comprobar rápidamente el tiempo que se necesita para resolver su problema. Un solucionador típica toma alrededor de O (3 * N ^ 3) flop (N es la dimensión de la matriz). Si n = 100.000, por lo tanto, usted está buscando en Gflops 3000000. Suponiendo que el solucionador en memoria hace 10 GFLOPS / s por núcleo, que busca a los 3 días de media en un solo núcleo. Como los algoritmos de escala bien, aumentando el número de núcleos debería reducir el tiempo cerca linealmente. Además de eso viene la E / S.

No está seguro acerca de las implementaciones de C ++, pero hay varias cosas que puede hacer si la memoria no es un problema dependiendo del tipo de matriz que está tratando con:

Si su matriz es escasa o bandas, se puede utilizar un solucionador de escasa o ancho de banda. Estos no almacenan cero elementos fuera de la banda.
Puede utilizar un solucionador de frente de onda, que almacena la matriz en el disco y sólo lleva en el frente de onda de la matriz para la descomposición.
Usted puede evitar la solución de la matriz por completo y utilizar métodos iterativos.
Puede probar métodos de Monte Carlo de solución.

Tener un vistazo a la lista de software de libre disposición para la solución de álgebra lineal problemas , compilado por Jack Dongarra y Hatem Ltaief.

Creo que para el tamaño del problema que estamos viendo, es probable que tenga un algoritmo iterativo. Si no desea almacenar la matriz A en un formato escaso, se puede utilizar una aplicación sin matriz. algoritmos iterativos típicamente no necesitan acceder a las entradas individuales de la matriz A, que sólo necesitan para calcular la matriz-vector productos Av (y algunas veces A ^ T v, el producto de la matriz transpuesta con el vector). Así que si la biblioteca está bien diseñado, que debería ser suficiente si se le pasa una clase que sabe cómo hacer productos matriz-vector.

A medida que la respuesta aceptada sugiere que hay UMFPACK. Pero si está utilizando BOOST puede seguir utilizando las matrices compactas en impulsar y utilizar UMFPACK para resolver el sistema. Hay una unión que hace que sea fácil:

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

Sus cerca de dos años anticuadas, pero es sólo una unión (junto con algunos otros).

véase pregunta relacionada: UMFPACK y uBLAS matrices dispersas de BOOST

Licenciado bajo: CC-BY-SA con atribución

No afiliado a StackOverflow