superlu_stdcplx.cpp

Go to the documentation of this file.

/* (w) Andreas Ahland */
/* $Id: superlu_stdcplx.cpp,v 1.1.2.18 2022/11/03 20:40:08 garloff Exp $ */

#include "tbci/cscmatrix.h"
#include "tbci/std_cplx.h"

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

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

#ifdef HAVE_UNISTD_H
#include <unistd.h>
#endif

NAMESPACE_TBCI

int lu_solve(CSCMatrix< CPLX__ complex<double> >& M,
             Vector< CPLX__ complex<double> >& x,
             const Vector< CPLX__ complex<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) {
        int newout = dup(savedstdout); close(savedstdout);
        if (newout != STDOUT_FILENO)
            fprintf(stderr, "Could not reopen stdout (%i)\n", newout);
    }
#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",
                        mem_usage.for_lu/1e6, mem_usage.total_needed/1e6);
#ifdef HAVE_MEM_USAGE_EXPANSIONS
                fprintf(stderr, "\texpansions %d\n",
                        mem_usage.expansions);
#else
                fprintf(stderr, "\n");
#endif
        } 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__ complex<double> > lu_solve(CSCMatrix< CPLX__ complex<double> >& M,
                                           const Vector< CPLX__ complex<double> >& rhs,
                                           colperm_t permc_spec,
                                           bool verbose, bool symm)
)
{
    TVector< CPLX__ complex<double> > x(rhs.size());
    const int info = lu_solve(M, (Vector< CPLX__ complex<double> >&)x,
                              rhs, permc_spec, verbose, symm);
    return x;
}
#endif

NAMESPACE_END