superlu_z.cpp

Go to the documentation of this file.

/* (w) Andreas Ahland */
/* $Id: superlu_z.cpp,v 1.4.2.12 2009/04/01 12:17:48 garloff Exp $ */

#include "cscmatrix.h"
#include "cplx.h"

#ifdef HAVE_SUPERLU3_INCLUDES
# include <superlu/slu_zdefs.h>
# include <superlu/supermatrix.h>
# include <superlu/slu_util.h>
#else
# include "superlu/util.h"
# include "superlu/zsp_defs.h"
#endif

#if TBCI_NUMLIB_VERSION >= 20600
# define COLPERM_T_DECLARED
# include "solver/superlu.h"
#endif

#ifdef HAVE_UNISTD_H
#include "unistd.h"
#endif

NAMESPACE_TBCI


int lu_solve(CSCMatrix< cplx<double> >& M, 
             Vector< cplx<double> >& x, 
             const Vector< cplx<double> >& rhs, 
             colperm_t permc_spec,
             bool verbose, bool symm)
{
    int nrhs = 1;
    int m    = M.rows();
    int n    = M.columns();
    int nnz  = M.size();
    doublecomplex *a= (doublecomplex*) M.dataPointer();
    unsigned int *asub= M.rowIndexPointer(); 
    unsigned int *xa  = M.columnPointer();

    int      *perm_r; /* row permutations from partial pivoting */
    int      *perm_c; /* column permutation vector */
    SuperMatrix A;
    SuperMatrix B;
    SuperMatrix L;      /* factor L */
    SuperMatrix U;      /* factor U */
    //char     trans[1];
    int  info;

#ifdef HAVE_SUPERLU3_INCLUDES
    superlu_options_t slu_opt;
    SuperLUStat_t slu_stat;
    set_default_options(&slu_opt);
    slu_opt.ColPerm = permc_spec;
    slu_opt.SymmetricMode = (yes_no_t)symm;
    if (symm)
        slu_opt.DiagPivotThresh = 0.1;    
    slu_opt.PrintStat = (yes_no_t)verbose;
    StatInit(&slu_stat);
#elif defined(HAVE_UNISTD_H)
    int savedstdout = 0;
#endif
 
    // Set up matrix structure
    x = rhs;
    //*trans = 'N';
    zCreate_CompCol_Matrix(&A, m, n, nnz, a, (int*)asub, (int*)xa, SLU_NC, SLU_Z, SLU_GE);
    zCreate_Dense_Matrix(&B, m, nrhs, (doublecomplex*) x.get_fortran_vector(),
                         m, SLU_DN, SLU_Z, SLU_GE);
    if ( !(perm_r = intMalloc(m)) ) 
        ABORT("Malloc fails for perm_r[].");
    if ( !(perm_c = intMalloc(n)) ) 
        ABORT("Malloc fails for perm_c[].");

#if defined(HAVE_UNISTD_H) && !defined(HAVE_SUPERLU3_INCLUDES)
    // HACK: SuperLU is way too talkative. Shut it down ...
    if (!verbose) {
        savedstdout = dup(STDOUT_FILENO); close(STDOUT_FILENO);
    }
#endif
#ifndef HAVE_SUPERLU3_INCLUDES    
    // Perform permutation
    get_perm_c(permc_spec, &A, perm_c);
#endif

    // Solve
#ifdef HAVE_SUPERLU3_INCLUDES
    zgssv(&slu_opt, &A, perm_c, perm_r, &L, &U, &B, &slu_stat, &info);
    if (verbose)
            StatPrint(&slu_stat);
#else
    zgssv(&A, perm_c, perm_r, &L, &U, &B, &info);
#endif
#if defined(HAVE_UNISTD_H) && !defined(HAVE_SUPERLU3_INCLUDES)
    // Restore stdout FD
    if (!verbose) {
        dup(savedstdout); close(savedstdout);
    }
#endif
    
    if (info != 0) {
        fprintf(stderr, "SuperLU: zgssv() error returns INFO= %d\n", info);
        if (info > n) { 
                /* We ran out of memory */
                mem_usage_t   mem_usage;
#ifdef HAVE_SUPERLU3_INCLUDES
                zQuerySpace(&L, &U, &mem_usage);
#else
                int panel_size = sp_ienv(1);
                zQuerySpace(&L, &U, panel_size, &mem_usage);
#endif
                fprintf(stderr, "SuperLU: Memory allocation of %i bytes failed!\n",
                        info - n);
                fprintf(stderr, "  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 
                fprintf(stderr, "SuperLU: Singular factor U(%i,%i) found!\n", info, info);
    }

    SUPERLU_FREE(perm_r);
    SUPERLU_FREE(perm_c);
    Destroy_SuperNode_Matrix(&L);
    Destroy_CompCol_Matrix(&U);
    // Huh, memleak :-O
    Destroy_SuperMatrix_Store(&B);
    Destroy_SuperMatrix_Store(&A);
#ifdef HAVE_SUPERLU3_INCLUDES
    StatFree(&slu_stat);
#endif

    return info;
}

#if TBCI_NUMLIB_VERSION < 20600

WEAK(
TVector< cplx<double> > lu_solve(CSCMatrix< cplx<double> >& M, 
                                 const Vector< cplx<double> >& rhs, 
                                 colperm_t permc_spec,
                                 bool verbose, bool symm)
)
{
    TVector< cplx<double> > x(rhs.size());
    const int info = lu_solve(M, (Vector< cplx<double> >&)x,
                              rhs, permc_spec, verbose, symm);
    return x;
}
#endif
NAMESPACE_END

---