gecrs-superlu.cc (flens/examples/gecrs-superlu.cc)

       1
       2
       3
       4
       5
       6
       7
       8
       9
      10
      11
      12
      13
      14
      15
      16
      17
      18
      19
      20
      21
      22
      23
      24
      25
      26
      27
      28
      29
      30
      31
      32
      33
      34
      35
      36
      37
      38
      39
      40
      41
      42
      43
      44
      45
      46
      47
      48
      49
      50
      51
      52
      53
      54
      55
      56
      57
      58
      59
      60
      61
      62
      63
      64
      65
      66
      67
      68
      69
      70
      71
      72
      73
      74
      75
      76
      77
      78
      79
      80
      81
      82
      83
      84
      85
      86
      87
      88
      89
      90
      91
      92
      93
      94
      95
      96
      97
      98
      99
     100
     101
     102
     103
     104
     105
     106
     107
     108
     109
     110
     111
     112
     113
     114
     115
     116
     117
     118
     119
     120
     121
     122
     123
     124
     125
     126
     127
     128
     129
     130
     131
     132
     133
     134
     135
     136
     137
     138
     139
     140
     141
     142
     143
     144
     145
     146
     147
     148
     149
     150
     151
     152
     153
     154
     155
     156
     157
     158
     159
     160
     161
     162
     163
     164
     165
     166
     167
     168
     169
     170
     171
     172
     173
     174
     175
     176

#include <iostream>
#include <flens/flens.cxx>

using namespace flens;
using namespace std;

///
/// Include the SuperLu stuff for double precision.
///
#include <slu_ddefs.h>

///
/// Again we hack a quick and dirty wrapper for SuperLU.  This time we call the
/// SuperLU expert driver routine *dgssvx*: As our sparse matrix comes in
/// compressed row storage format but SuperLU expects the compressed col
/// storage format we actually habe to solve $A^T X = B$ instead of $AX = B$.
/// This operation is not available when using $dgssvx'.  The wrapper is based
/// on `dlinsolx.c` from the SuperLU examples directory.
///

template <typename MA, typename PR, typename PC, typename VB>
int
dgssvx(GeCRSMatrix<MA>  &A,
       DenseVector<PC>  &pc,
       DenseVector<PR>  &pr,
       DenseVector<VB>  &b)
{
///
/// We keep it as simple as possible and assume $A$ is square.
///
    ASSERT(A.numRows()==A.numCols());
    ASSERT(b.length()==A.numRows());
    const int n = A.numRows();

///
/// Some *temporary* vectors for intermediate results.
///
    typedef typename DenseVector<PC>::NoView  TempIntVec;
    typedef typename DenseVector<VB>::NoView  TempDoubleVec;
    TempDoubleVec   x(n), r(n), c(n), ferr(1), berr(1);
    TempIntVec      etree(n);

    superlu_options_t   options;
    SuperLUStat_t       stat;
    SuperMatrix         _A, _L, _U, _B, _X;

///
/// Set the transpose flag so that we solve $A^T X = B$.  Check the SuperLU
/// documentation for further details.
///
    set_default_options(&options);
    options.Equil           = YES;
    options.Trans           = TRANS;
    options.DiagPivotThresh = double(1);
    // Add more functionalities than the defaults.
    options.PivotGrowth = YES;        // Compute reciprocal pivot growth
    options.ConditionNumber = YES;    // Compute reciprocal condition number
    options.IterRefine = SLU_DOUBLE;  // Perform double-precision refinement

///
/// Because `A` has *CRS* format first pass `A.engine().cols().data()`and then
/// `A.engine().rows().data()`.  After creation `A_` holds a reference of $A^T$.
///
    dCreate_CompCol_Matrix(&_A, n, n, A.engine().numNonZeros(),
                           A.engine().values().data(),
                           A.engine().cols().data(),
                           A.engine().rows().data(),
                           SLU_NC, SLU_D, SLU_GE);
    dCreate_Dense_Matrix(&_B, n, 1, b.data(), b.length(),
                         SLU_DN, SLU_D, SLU_GE);
    dCreate_Dense_Matrix(&_X,
                         x.length(), 1, x.data(), x.length(),
                         SLU_DN, SLU_D, SLU_GE);

    StatInit(&stat);
    int             info;
    char            equed;
    int             lwork = 0; // let the solver allocate its own memory
    void            *work = 0;
    double          rpg, rcond;
    mem_usage_t     mem_usage;

///
/// Solve the system and compute the condition number and error bounds using
/// dgssvx.  Except for the computed solution this trashy SuperLU wrapper
/// ignores most of these computations.  However, in a real application some
/// of this stuff could be really useful (e.g. `rcond`!).
///
    dgssvx(&options, &_A, pc.data(), pr.data(), etree.data(), &equed,
           r.data(), c.data(), &_L, &_U, work, lwork, &_B, &_X,
           &rpg, &rcond, ferr.data(), berr.data(), &mem_usage, &stat, &info);

///
/// Overwrite $b$ with the solution $x$. (Quick and dirty, remember?)
///
    b = x;

    Destroy_SuperMatrix_Store(&_A);
    Destroy_SuperMatrix_Store(&_B);
    Destroy_SuperNode_Matrix(&_L);
    Destroy_CompCol_Matrix(&_U);
    StatFree(&stat);
    return info;
}

int
main()
{
///
/// SuperLU requires an index base of zero.  Here we set the default index base
/// of the storage scheme to zero via `IndexBaseZero`.  Note that the template
/// *Coord* contains the *CoordRowColCmp* class as third parameter.  That's
/// because we want to convert the storage to *compressed row storage* later.
///
    typedef int                                              IndexType;
    typedef IndexBaseZero<IndexType>                         IndexBase;
    typedef CoordStorage<double, CoordRowColCmp, IndexBase>  Coord;

    const IndexType m = 5;
    const IndexType n = 5;
    GeCoordMatrix<Coord>  A_(m, n);

///
/// Again we setup the matrix from the SuperLU user guide.  So here the values.
///
    const double s = 19,
                 u = 21,
                 p = 16,
                 e =  5,
                 r = 18,
                 l = 12;
    A_(0,0) += s;
    A_(1,1) += u;
    A_(2,2) += p;
    A_(3,3) += e;
    A_(4,4) += r;
    A_(1,0) += l; A_(2,1) += l; A_(4,0) += l; A_(4,1) += l;
    A_(0,2) += u; A_(0,3) += u; A_(3,4) += u;

///
/// Convert to *compressed row storage*.
///
    GeCRSMatrix<CRS<double, IndexBase> >  A = A_;

///
/// Just for curiosity: Compare the two formats.
///
    std::cout << "A_ = " << A_ << endl;
    std::cout << "A = " << A << endl;

///
/// Setup the right-hand side $b$.  We set all values to $1$.
///
    DenseVector<Array<double, IndexBase> >   b(m);
    b = 1;

///
/// Call the SuperLU solver for solving $Ax = b$.  Note that on exit $b$ is
/// overwritten with the solution $x$.
///
    DenseVector<Array<IndexType, IndexBase> >  pr(m), pc(n);
    int info = dgssvx(A, pr, pc, b);

///
/// Check the info code
///
    if (info==0) {
        cout << "x = " << b << endl;
    } else {
        cout << "SuperLU dgssv:  info = " << info << endl;
    }

    return info;
}