b2superlu_solver.H

//------------------------------------------------------------------------
// b2superlu_solver.H --
//
//
// written by Mathias Doreille
//
// Copyright (c) 2004-2012 SMR Engineering & Development SA
//                         2502 Bienne, Switzerland
//
// All Rights Reserved.  Proprietary source code.  The contents of
// this file may not be disclosed to third parties, copied or
// duplicated in any form, in whole or in part, without the prior
// written permission of SMR.
//------------------------------------------------------------------------
 
#ifndef __B2SUPERLU_SOLVER_H__
#define __B2SUPERLU_SOLVER_H__
 
extern "C" {
#include <SuperLU/dsp_defs.h>
}
 
#include "b2sparse_solver.H"
 
namespace b2000 { namespace b2linalg {
class Generic_SuperLU {
public:
    void Init_SuperMatrix_Store(SuperMatrix* A) { A->Store = 0; }
 
    void Create_CompCol_Matrix(
          SuperMatrix* A, int m, int n, int nnz, double* nzval, int* rowind, int* colptr,
          Stype_t stype, Mtype_t mtype) {
        dCreate_CompCol_Matrix(A, m, n, nnz, nzval, rowind, colptr, stype, SLU_D, mtype);
    }
 
    void Destroy_CompCol_Matrix(SuperMatrix* A) {
        if (A->Store) {
            SUPERLU_FREE(((NCformat*)A->Store)->rowind);
            SUPERLU_FREE(((NCformat*)A->Store)->colptr);
            SUPERLU_FREE(((NCformat*)A->Store)->nzval);
        }
        SUPERLU_FREE(A->Store);
    }
 
    void Create_Dense_Matrix(
          SuperMatrix* X, int m, int n, double* x, int ldx, Stype_t stype, Mtype_t mtype) {
        dCreate_Dense_Matrix(X, m, n, x, ldx, stype, SLU_D, mtype);
    }
 
    void Destroy_Dense_Matrix(SuperMatrix* A) {
        if (A->Store) { SUPERLU_FREE(((DNformat*)A->Store)->nzval); }
        SUPERLU_FREE(A->Store);
    }
 
    void Destroy_SuperNode_Matrix(SuperMatrix* A) {
        if (A->Store) {
            SUPERLU_FREE(((SCformat*)A->Store)->rowind);
            SUPERLU_FREE(((SCformat*)A->Store)->rowind_colptr);
            SUPERLU_FREE(((SCformat*)A->Store)->nzval);
            SUPERLU_FREE(((SCformat*)A->Store)->nzval_colptr);
            SUPERLU_FREE(((SCformat*)A->Store)->col_to_sup);
            SUPERLU_FREE(((SCformat*)A->Store)->sup_to_col);
        }
        SUPERLU_FREE(A->Store);
    }
 
    void gssvx(
          superlu_options_t* options, SuperMatrix* A, int* perm_c, int* perm_r, int* etree,
          char* equed, double* R, double* C, SuperMatrix* L, SuperMatrix* U, void* work, int lwork,
          SuperMatrix* B, SuperMatrix* X, double* recip_pivot_growth, double* rcond, double* ferr,
          double* berr, mem_usage_t* mem_usage, SuperLUStat_t* stat, int* info) {
        dgssvx(
              options, A, perm_c, perm_r, etree, equed, R, C, L, U, work, lwork, B, X,
              recip_pivot_growth, rcond, ferr, berr, mem_usage, stat, info);
    }
};
 
template <typename T>
class SuperLU_LDLt_sparse_direct_solver : public LDLt_sparse_solver<T>, public Generic_SuperLU {
public:
    SuperLU_LDLt_sparse_direct_solver()
        : size(0),
          nnz(0),
          colind(0),
          rowind(),
          value(0),
          colptr1(0),
          rowind1(0),
          nzval1(0),
          perm_c(0),
          perm_r(0),
          etree(0),
          equed(0),
          R(0),
          C(0) {
        Init_SuperMatrix_Store(&L);
        Init_SuperMatrix_Store(&U);
        Init_SuperMatrix_Store(&A);
    }
 
    ~SuperLU_LDLt_sparse_direct_solver() {
        delete[] perm_c;
        delete[] perm_r;
        delete[] etree;
        delete[] R;
        delete[] C;
        Destroy_SuperNode_Matrix(&L);
        Destroy_CompCol_Matrix(&U);
        Destroy_CompCol_Matrix(&A);
    }
 
    void init(
          size_t size_, size_t nnz_, const size_t* colind_, const size_t* rowind_, const T* value_,
          int connectivity, const Dictionary& dictionary) {
        size = size_;
        nnz = nnz_;
        colind = colind_;
        rowind = rowind_;
        value = value_;
 
        delete[] perm_c;
        perm_c = new int[size];
        delete[] perm_r;
        perm_r = new int[size];
        delete[] etree;
        etree = new int[size];
        equed = 'N';
        delete[] R;
        R = new T[size];
        delete[] C;
        C = new T[size];
 
        set_default_options(&options);
        options.SymmetricMode = YES;
        options.Equil = NO;
        // options.IterRefine = EXTRA;
        // options.DiagPivotThresh = 0.0;
        options.ConditionNumber = YES;
    }
 
    void update_value() { options.Fact = SamePattern; }
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx,
          char left_or_right = ' ') {
        if (left_or_right == 'L' || left_or_right == 'R') { UnimplementedError() << THROW; }
        SuperLUStat_t stat;
        StatInit(&stat);
        int info;
        mem_usage_t mem_usage;
        T rpg;
        T rcond;
        std::vector<T> ferr(s);
        std::vector<T> berr(s);
 
        SuperMatrix B, X;
        Create_Dense_Matrix(&B, s, nrhs, const_cast<double*>(b), ldb, SLU_DN, SLU_GE);
        Create_Dense_Matrix(&X, s, nrhs, x, ldx, SLU_DN, SLU_GE);
 
        if (options.Fact != FACTORED) {
            Destroy_SuperNode_Matrix(&L);
            Destroy_CompCol_Matrix(&U);
            Destroy_CompCol_Matrix(&A);
            colptr1 = (int*)SUPERLU_MALLOC((size + 2) * sizeof(int));
            colptr1[0] = colptr1[1] = 0;
            colptr1 += 2;
            for (size_t j = 0; j != size; ++j) { colptr1[j] = colind[j + 1] - colind[j]; }
            for (size_t j = 0; j != size; ++j) {
                const size_t i_end = colind[j + 1];
                for (size_t i = colind[j] + 1; i < i_end; ++i) { ++colptr1[rowind[i]]; }
            }
            --colptr1;
            for (size_t j = 0; j != size; ++j) { colptr1[j + 1] += colptr1[j]; }
            int nnz1 = colptr1[size];
            int* rowind1 = (int*)SUPERLU_MALLOC(nnz1 * sizeof(int));
            T* nzval1 = (T*)SUPERLU_MALLOC(nnz1 * sizeof(T));
            for (size_t j = 0; j != size; ++j) {
                size_t i = colind[j];
                const size_t i_end = colind[j + 1];
                int ii = colptr1[j];
                colptr1[j] += i_end - i;
                std::copy(rowind + i, rowind + i_end, rowind1 + ii);
                std::copy(value + i, value + i_end, nzval1 + ii);
                for (++i; i < i_end; ++i) {
                    int ii = colptr1[rowind[i]]++;
                    rowind1[ii] = j;
                    nzval1[ii] = value[i];
                }
            }
            --colptr1;
 
            /*
            for (size_t j = 0; j != size; ++j) {
                for (int i = colptr1[j]; i != colptr1[j+1]; ++i)
                    std::cout << "(" << rowind1[i] << "," << nzval1[i] << ") ";
                std::cout << std::endl;
            }
            std::cout << std::endl;
            */
 
            Create_CompCol_Matrix(&A, size, size, nnz1, nzval1, rowind1, colptr1, SLU_NC, SLU_GE);
        }
 
        gssvx(&options, &A, perm_c, perm_r, etree, &equed, R, C, &L, &U, 0, 0, &B, &X, &rpg, &rcond,
              &ferr[0], &berr[0], &mem_usage, &stat, &info);
 
        std::cout << "factorisation, rcond = " << rcond << std::endl;
 
        Destroy_SuperMatrix_Store(&B);
        Destroy_SuperMatrix_Store(&X);
        StatFree(&stat);
 
        if (info > 0 && info <= int(s)) {
            Exception() << "The matrix entrie L(" << info << "," << info
                        << ") of the factorised matrix is exactly zero. The factorization has been "
                           "completed, but the factor U is exactly singular, so the solution and "
                           "error bounds could not be computed."
                        << THROW;
        }
        if (info == int(s) + 1) {
            Exception() << "The factorised matrix L is nonsingular, but this condition number is "
                           "less than machine precision, meaning that the matrix is singular to "
                           "working precision."
                        << THROW;
        }
        if (info) { Exception() << THROW; }
 
        options.Fact = FACTORED;
    }
 
    void resolve(size_t s, size_t nrhs, T* bx, size_t ldbx, char left_or_right = ' ') {
        size_t ss = s * nrhs;
        T* x = new T[ss];
        resolve(s, nrhs, bx, ldbx, x, s, left_or_right);
        if (ldbx == s) {
            std::copy(x, x + ss, bx);
        } else {
            T* xx = x;
            T* xx_end = xx + ss;
            while (xx < xx_end) {
                std::copy(xx, xx + s, bx);
                xx += s;
                bx += ldbx;
            }
        }
        delete[] x;
    }
 
private:
    size_t size;
    size_t nnz;
    const size_t* colind;
    const size_t* rowind;
    const T* value;
 
    int* colptr1;
    int* rowind1;
    T* nzval1;
 
    SuperMatrix A, L, U;
    int* perm_c;
    int* perm_r;
    int* etree;
    char equed;
    T* R;
    T* C;
    superlu_options_t options;
};
 
template <typename T>
class SuperLU_LDLt_extension_sparse_direct_solver : public LDLt_extension_sparse_solver<T>,
                                                    public Generic_SuperLU {
public:
    SuperLU_LDLt_extension_sparse_direct_solver()
        : size(0),
          nnz(0),
          colind(0),
          rowind(),
          value(0),
          size_ext(0),
          colptr1(0),
          rowind1(0),
          nzval1(0),
          perm_c(0),
          perm_r(0),
          etree(0),
          equed(0),
          R(0),
          C(0) {
        Init_SuperMatrix_Store(&L);
        Init_SuperMatrix_Store(&U);
        Init_SuperMatrix_Store(&A);
    }
 
    ~SuperLU_LDLt_extension_sparse_direct_solver() {
        delete[] perm_c;
        delete[] perm_r;
        delete[] etree;
        delete[] R;
        delete[] C;
        Destroy_SuperNode_Matrix(&L);
        Destroy_CompCol_Matrix(&U);
        Destroy_CompCol_Matrix(&A);
    }
 
    void init(
          size_t size_, size_t nnz_, const size_t* colind_, const size_t* rowind_, const T* value_,
          size_t size_ext_, int connectivity, const Dictionary& dictionary) {
        size = size_;
        nnz = nnz_;
        colind = colind_;
        rowind = rowind_;
        value = value_;
        size_ext = size_ext_;
 
        delete[] perm_c;
        perm_c = new int[size + size_ext];
        delete[] perm_r;
        perm_r = new int[size + size_ext];
        delete[] etree;
        etree = new int[size + size_ext];
        equed = 'N';
        delete[] R;
        R = new T[size + size_ext];
        delete[] C;
        C = new T[size + size_ext];
 
        set_default_options(&options);
        options.SymmetricMode = YES;
        options.Equil = YES;
        options.IterRefine = EXTRA;
        // options.DiagPivotThresh = 0.0;
        options.ConditionNumber = YES;
    }
 
    void update_value() { options.Fact = SamePattern; }
 
    void resolve(
          size_t s, size_t nrhs, const T* b, size_t ldb, T* x, size_t ldx, const T* ma = 0,
          const T* mb = 0, const T* mc = 0, char left_or_right = ' ') {
        if (left_or_right == 'L' || left_or_right == 'R') { UnimplementedError() << THROW; }
        SuperLUStat_t stat;
        StatInit(&stat);
        int info;
        mem_usage_t mem_usage;
        T rpg;
        T rcond;
        std::vector<T> ferr(s);
        std::vector<T> berr(s);
 
        SuperMatrix B, X;
        Create_Dense_Matrix(&B, s, nrhs, const_cast<double*>(b), ldb, SLU_DN, SLU_GE);
        Create_Dense_Matrix(&X, s, nrhs, x, ldx, SLU_DN, SLU_GE);
 
        if (ma != 0 && mb != 0 && mc != 0 && options.Fact == FACTORED) {
            options.Fact = SamePattern;
        }
 
        if (options.Fact != FACTORED) {
            Destroy_SuperNode_Matrix(&L);
            Destroy_CompCol_Matrix(&U);
            Destroy_CompCol_Matrix(&A);
            colptr1 = (int*)SUPERLU_MALLOC((size + size_ext + 1) * sizeof(int));
            colptr1[0] = colptr1[1] = 0;
            colptr1 += 2;
            for (size_t j = 0; j != size; ++j) {
                colptr1[j] = colind[j + 1] - colind[j] + size_ext;
            }
            for (size_t j = 0; j != size; ++j) {
                const size_t i_end = colind[j + 1];
                for (size_t i = colind[j] + 1; i < i_end; ++i) { ++colptr1[rowind[i]]; }
            }
            --colptr1;
            for (size_t j = 0; j != size; ++j) { colptr1[j + 1] += colptr1[j]; }
            int nnz1 = colptr1[size] + size + 1;
            int* rowind1 = (int*)SUPERLU_MALLOC(nnz1 * sizeof(int));
            T* nzval1 = (T*)SUPERLU_MALLOC(nnz1 * sizeof(T));
            for (size_t j = 0; j != size; ++j) {
                size_t i = colind[j];
                const size_t i_end = colind[j + 1];
                int ii = colptr1[j];
                colptr1[j] += i_end - i;
                std::copy(rowind + i, rowind + i_end, rowind1 + ii);
                std::copy(value + i, value + i_end, nzval1 + ii);
                for (++i; i < i_end; ++i) {
                    int ii = colptr1[rowind[i]]++;
                    rowind1[ii] = j;
                    nzval1[ii] = value[i];
                }
            }
 
            // Add the extension matrix mb, ma and mc
            if (ma == 0 || mb == 0 || mc == 0) { Exception() << THROW; }
            for (size_t j = 0; j != size; ++j) {
                int& ii = colptr1[j];
                for (size_t i = 0; i != size_ext; ++i, ++ii) {
                    rowind1[ii] = size + i;
                    nzval1[ii] = mb[j + i * size];
                }
            }
 
            --colptr1;
 
            {
                for (size_t i = 0; i != size_ext; ++i) {
                    size_t ii = colptr1[size + i];
                    colptr1[size + i + 1] = ii + size + size_ext;
                    for (size_t j = 0; j != size; ++j, ++ii) {
                        rowind1[ii] = j;
                        nzval1[ii] = ma[j + i * size];
                    }
                    for (size_t j = 0; j != size_ext; ++j, ++ii) {
                        rowind1[ii] = j + size;
                        nzval1[ii] = mc[j + i * size_ext];
                    }
                }
            }
 
            for (size_t j = 0; j != size + size_ext; ++j) {
                for (int i = colptr1[j]; i != colptr1[j + 1]; ++i) {
                    std::cout << "(" << rowind1[i] << "," << nzval1[i] << ") ";
                }
                std::cout << std::endl;
            }
            std::cout << std::endl;
 
            Create_CompCol_Matrix(
                  &A, size + size_ext, size + size_ext, nnz1, nzval1, rowind1, colptr1, SLU_NC,
                  SLU_GE);
        }
 
        gssvx(&options, &A, perm_c, perm_r, etree, &equed, R, C, &L, &U, 0, 0, &B, &X, &rpg, &rcond,
              &ferr[0], &berr[0], &mem_usage, &stat, &info);
 
        std::cout << "factorisation, rcond = " << rcond << std::endl;
 
        Destroy_SuperMatrix_Store(&B);
        Destroy_SuperMatrix_Store(&X);
        StatFree(&stat);
 
        if (info > 0 && info <= int(s)) {
            Exception() << "The matrix entrie L(" << info << "," << info
                        << ") of the factorised matrix is exactly zero. The factorization has been "
                           "completed, but the factor U is exactly singular, so the solution and "
                           "error bounds could not be computed."
                        << THROW;
        }
        if (info == int(s) + 1) {
            Exception() << "The factorised matrix L is nonsingular, but this condition number is "
                           "less than machine precision, meaning that the matrix is singular to "
                           "working precision."
                        << THROW;
        }
        if (info) { Exception() << THROW; }
 
        options.Fact = FACTORED;
    }
 
    void resolve(
          size_t s, size_t nrhs, T* bx, size_t ldbx, const T* ma = 0, const T* mb = 0,
          const T* mc = 0, char left_or_right = ' ') {
        size_t ss = s * nrhs;
        T* x = new T[ss];
        resolve(s, nrhs, bx, ldbx, x, s, ma, mb, mc, left_or_right);
        if (ldbx == s) {
            std::copy(x, x + ss, bx);
        } else {
            T* xx = x;
            T* xx_end = xx + ss;
            while (xx < xx_end) {
                std::copy(xx, xx + s, bx);
                xx += s;
                bx += ldbx;
            }
        }
        delete[] x;
    }
 
private:
    size_t size;
    size_t nnz;
    const size_t* colind;
    const size_t* rowind;
    const T* value;
 
    size_t size_ext;
 
    int* colptr1;
    int* rowind1;
    T* nzval1;
 
    SuperMatrix A, L, U;
    int* perm_c;
    int* perm_r;
    int* etree;
    char equed;
    T* R;
    T* C;
    superlu_options_t options;
};
 
}}  // namespace b2000::b2linalg
 
#endif