SuperLU

De Wiki de Calcul Québec
Aller à : Navigation, rechercher
Autres langues :anglais 100% • ‎français 100%

Description

SuperLU (Supernodal LU) est une bibliothèque mathématique pour la résolution directe de grands systèmes d'équations linéaires représentés par des matrices creuses. Tel que son nom l'indique la résolution procède par décomposition LU. Cette partie du code peut également traiter des matrices non carrées.

La bibliothèque est écrite en C, mais peut également être utilisée en Fortran. Il existe deux variantes parallèles à SuperLU : SuperLU_MT pour les systèmes à mémoire partagée et SuperLU_DIST pour les systèmes à mémoire distribuée. Cette bibliothèque utilise la bibliothèque MPI pour la parallélisation.

Un guide pour les utilisateurs est disponible en version PDF sur le site de Lawrence Berkeley National Laboratory. La bibliothèque requiert la présence d'une bibliothèque BLAS.

Exemple

Le fichier suivant contient un exemple en langage C de résolution d'équations linéaires complexes avec la bibliothèque SuperLU :

Fichier : ex_superlu.c
/*
 * -- SuperLU routine (version 3.0) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * October 15, 2003
 *
 */
#include "slu_zdefs.h"
 
main(int argc, char *argv[])
{
    SuperMatrix A;
    NCformat *Astore;
    doublecomplex   *a;
    int      *asub, *xa;
    int      *perm_c; /* column permutation vector */
    int      *perm_r; /* row permutations from partial pivoting */
    SuperMatrix L;      /* factor L */
    SCformat *Lstore;
    SuperMatrix U;      /* factor U */
    NCformat *Ustore;
    SuperMatrix B;
    int      nrhs, ldx, info, m, n, nnz;
    doublecomplex   *xact, *rhs;
    mem_usage_t   mem_usage;
    superlu_options_t options;
    SuperLUStat_t stat;
 
#if ( DEBUGlevel>=1 )
    CHECK_MALLOC("Enter main()");
#endif
 
    /* Set the default input options:
	options.Fact = DOFACT;
        options.Equil = YES;
    	options.ColPerm = COLAMD;
	options.DiagPivotThresh = 1.0;
    	options.Trans = NOTRANS;
    	options.IterRefine = NOREFINE;
    	options.SymmetricMode = NO;
    	options.PivotGrowth = NO;
    	options.ConditionNumber = NO;
    	options.PrintStat = YES;
     */
    set_default_options(&options);
 
    /* Read the matrix in Harwell-Boeing format. */
    zreadhb(&m, &n, &nnz, &a, &asub, &xa);
 
    zCreate_CompCol_Matrix(&A, m, n, nnz, a, asub, xa, SLU_NC, SLU_Z, SLU_GE);
    Astore = A.Store;
    printf("Dimension %dx%d; # nonzeros %d\n", A.nrow, A.ncol, Astore->nnz);
 
    nrhs   = 1;
    if ( !(rhs = doublecomplexMalloc(m * nrhs)) ) ABORT("Malloc fails for rhs[].");
    zCreate_Dense_Matrix(&B, m, nrhs, rhs, m, SLU_DN, SLU_Z, SLU_GE);
    xact = doublecomplexMalloc(n * nrhs);
    ldx = n;
    zGenXtrue(n, nrhs, xact, ldx);
    zFillRHS(options.Trans, nrhs, xact, ldx, &A, &B);
 
    if ( !(perm_c = intMalloc(n)) ) ABORT("Malloc fails for perm_c[].");
    if ( !(perm_r = intMalloc(m)) ) ABORT("Malloc fails for perm_r[].");
 
    /* Initialize the statistics variables. */
    StatInit(&stat);
 
    zgssv(&options, &A, perm_c, perm_r, &L, &U, &B, &stat, &info);
 
    if ( info == 0 ) {
 
	/* This is how you could access the solution matrix. */
        doublecomplex *sol = (doublecomplex*) ((DNformat*) B.Store)->nzval; 
 
	 /* Compute the infinity norm of the error. */
	zinf_norm_error(nrhs, &B, xact);
 
	Lstore = (SCformat *) L.Store;
	Ustore = (NCformat *) U.Store;
    	printf("No of nonzeros in factor L = %d\n", Lstore->nnz);
    	printf("No of nonzeros in factor U = %d\n", Ustore->nnz);
    	printf("No of nonzeros in L+U = %d\n", Lstore->nnz + Ustore->nnz - n);
 
	zQuerySpace(&L, &U, &mem_usage);
	printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
	       mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
	       mem_usage.expansions);
 
    } else {
	printf("zgssv() error returns INFO= %d\n", info);
	if ( info <= n ) { /* factorization completes */
	    zQuerySpace(&L, &U, &mem_usage);
	    printf("L\\U MB %.3f\ttotal MB needed %.3f\texpansions %d\n",
		   mem_usage.for_lu/1e6, mem_usage.total_needed/1e6,
		   mem_usage.expansions);
	}
    }
 
    if ( options.PrintStat ) StatPrint(&stat);
    StatFree(&stat);
 
    SUPERLU_FREE (rhs);
    SUPERLU_FREE (xact);
    SUPERLU_FREE (perm_r);
    SUPERLU_FREE (perm_c);
    Destroy_CompCol_Matrix(&A);
    Destroy_SuperMatrix_Store(&B);
    Destroy_SuperNode_Matrix(&L);
    Destroy_CompCol_Matrix(&U);
 
#if ( DEBUGlevel>=1 )
    CHECK_MALLOC("Exit main()");
#endif
}


Il est possible de compiler l'exemple avec une ligne de commande ressemblant à ceci, si la bibliothèque SuperLU se trouve dans vos modules :


[nom@serveur $] 
icc ex_superlu.c -lsuperlu -lmkl_lapack95_lp64 -lmkl_intel_lp64 -lmkl_intel_thread -lmkl_core -liomp5 -lpthread -lm


Outils personnels
Espaces de noms

Variantes
Actions
Navigation
Ressources de Calcul Québec
Outils
Partager