superlu_d.cpp

Go to the documentation of this file.

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

#include "tbci/cscmatrix.h"

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

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

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

NAMESPACE_TBCI

int lu_solve(CSCMatrix<double>& M,
             Vector<double>& x,
             const Vector<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();
    double *a= 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 */
    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;
    dCreate_CompCol_Matrix(&A, m, n, nnz, a, (int*)asub, (int*)xa, SLU_NC, SLU_D, SLU_GE);
    dCreate_Dense_Matrix(&B, m, nrhs, x.get_fortran_vector(), m, SLU_DN, SLU_D, 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
    dgssv(&slu_opt, &A, perm_c, perm_r, &L, &U, &B, &slu_stat, &info);
    if (verbose)
            StatPrint(&slu_stat);
#else
    dgssv(&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: dgssv() error returns INFO= %d\n", info);
        if (info > n) {
                /* We ran out of memory */
                mem_usage_t   mem_usage;
#ifdef HAVE_SUPERLU3_INCLUDES
                dQuerySpace(&L, &U, &mem_usage);
#else
                int panel_size = sp_ienv(1);
                dQuerySpace(&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<double> lu_solve(CSCMatrix<double>& M,
                         const Vector<double>& rhs,
                         colperm_t permc_spec,
                         bool verbose, bool symm)
)
{
    TVector<double> x(rhs.size());
    const int info = lu_solve(M, (Vector<double>&)x, rhs, permc_spec, verbose, symm);
    return x;
}
#endif

NAMESPACE_END