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; } |