Risposta breve: Non utilizzare i binding di Boost LAPACK, questi sono stati progettati per matrici dense, matrici non sparse, utilizzare invece UMFPACK.
Risposta lunga: UMFPACK è una delle migliori librerie per risolvere Ax = b quando A è grande e sparse.
seguito è il codice di esempio (basato su 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 funzione generate_sparse_matrix_problem crea la matrice A e il lato destro b. La matrice viene prima costruita in forma di terzine. I vettori Ti, Tj e Tx descrivono completamente A. La forma di triplet è facile da creare ma i metodi di matrice sparse efficienti richiedono il formato di colonna sparsa compressa. La conversione viene eseguita con umfpack_di_triplet_to_col.
Una fattorizzazione simbolica viene eseguita con umfpack_di_symbolic. Una scarsa decomposizione LU di A viene eseguita con umfpack_di_numeric. Le soluzioni triangolari inferiore e superiore vengono eseguite con umfpack_di_solve.
Con n su 500.000, sulla mia macchina, l'intero programma impiega circa un secondo per essere eseguito. Valgrind segnala che sono stati allocati 369.239.649 byte (poco più di 352 MB).
Nota questo page descrive il supporto di Boost per le matrici sparse in Triplet (Coordinate) e il formato compresso. Se lo si desidera, è possibile scrivere routine per convertire questi oggetti boost nei matrici semplici UMFPACK come input.
Indicando l'ovvio, una matrice che è matrice di dimensione 4x10^12 a 4x10^14 byte o 4 a 400 terabyte per solo una singola matrice. (a meno che, come indicato sotto, sia scarso) – cyberconte