HermesCommon  3.0
superlu_solver.h
Go to the documentation of this file.
1 // This file is part of HermesCommon
2 //
3 // Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
5 //
6 // Hermes2D is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published
8 // by the Free Software Foundation; either version 2 of the License,
9 // or (at your option) any later version.
10 //
11 // Hermes2D is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Hermes2D; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 #ifndef __HERMES_COMMON_SUPERLU_SOLVER_H_
23 #define __HERMES_COMMON_SUPERLU_SOLVER_H_
24
25 #include "config.h"
26 #ifdef WITH_SUPERLU
27 typedef int int_t;
29 #include "algebra/cs_matrix.h"
30
31 #include <supermatrix.h>
32 #include <slu_util.h>
33 namespace Hermes
34 {
35  namespace Solvers
36  {
37  template <typename Scalar> class SuperLUSolver;
38 #ifdef SLU_MT
39  template <typename Scalar>
40  class SuperLu
41  {
42  public:
43  void gsequ(SuperMatrix *A, double *r, double *c, double *rowcnd, double *colcnd, double *amax, int *info);
44  void laqgs(SuperMatrix *A, float *r, float *c, float rowcnd, float colcnd, float amax, char *equed);
45  int gstrf(superlu_options_t *options, int m, int n, double anorm, LUstruct_t *LUstruct, gridinfo_t *grid, SuperLUStat_t *stat, int *info);
46  float pivotGrowth(int ncols, SuperMatrix *A, int *perm_c, SuperMatrix *L, SuperMatrix *U);
47  float langs(char *norm, SuperMatrix *A);
48  void gscon(char *norm, SuperMatrix *L, SuperMatrix *U, float anorm, float *rcond, SuperLUStat_t *stat, int *info);
49  void gstrs(trans_t trans, SuperMatrix *L, SuperMatrix *U, int *perm_c, int *perm_r, SuperMatrix *B, SuperLUStat_t *stat, int *info);
50  double lamch_(char *cmach);
51  int querySpace(SuperMatrix *, SuperMatrix *, slu_memusage_t *);
52  }
53 #else //SLU_MT
54
55  typedef superlu_options_t slu_options_t;
56  typedef SuperLUStat_t slu_stat_t;
57  typedef struct
58  {
59  float for_lu;
60  float total_needed;
61  } slu_memusage_t;
62 #define SLU_DESTROY_L Destroy_SuperNode_Matrix
63 #define SLU_DESTROY_U Destroy_CompCol_Matrix
64 #define SLU_INIT_STAT(stat_ptr) StatInit(stat_ptr)
65 #define SLU_PRINT_STAT(stat_ptr) StatPrint(stat_ptr)
66
67 #define SLU_DTYPE SLU_D
68
69 #define SLU_PRINT_CSC_MATRIX zPrint_CompCol_Matrix
70 #define Scalar_MALLOC doublecomplexMalloc
71
73  template<typename Scalar> struct SuperLuType;
74
76  template<>
77  struct SuperLuType < double >
78  {
80  typedef double Scalar;
81  };
82
84  template<>
85  struct SuperLuType < std::complex<double> >
86  {
88  typedef struct { double r, i; } Scalar;
89  };
90 #endif //SLU_MT
91  }
92 }
93
94 namespace Hermes
95 {
96  namespace Solvers
97  {
100  template <typename Scalar>
101  class HERMES_API SuperLUSolver : public DirectSolver < Scalar >
102  {
103  public:
107  SuperLUSolver(CSCMatrix<Scalar> *m, SimpleVector<Scalar> *rhs);
108  virtual ~SuperLUSolver();
109
110  virtual void solve();
111  virtual int get_matrix_size();
112  void free();
113
114  protected:
116  CSCMatrix<Scalar> *m;
117
119  SimpleVector<Scalar> *rhs;
120
122  bool has_A, has_B;
124  bool inited;
126  bool A_changed;
127  // internally during factorization or externally by
128  // the user.
129
132  bool check_status(unsigned int info);
133
136  int *local_Ai, *local_Ap;
137  typename SuperLuType<Scalar>::Scalar *local_Ax, *local_rhs;
138
139  bool setup_factorization();
140  void free_factorization_data();
141  void free_matrix();
142  void free_rhs();
143
145  SuperMatrix A, B;
147  SuperMatrix L, U;
149  double *R, *C;
151  int *perm_r;
153  int *perm_c;
155  int *etree;
157  slu_options_t options;
158
159  private:
160 #ifndef SLU_MT
161  void create_csc_matrix(SuperMatrix *A, int m, int n, int nnz, typename SuperLuType<Scalar>::Scalar *nzval, int *rowind, int *colptr,
162  Stype_t stype, Dtype_t dtype, Mtype_t mtype);
163  void solver_driver(superlu_options_t *options, SuperMatrix *A, int *perm_c, int *perm_r, int *etree, char *equed, double *R,
164  double *C, SuperMatrix *L, SuperMatrix *U, void *work, int lwork, SuperMatrix *B, SuperMatrix *X, double *recip_pivot_growth,
165  double *rcond, double *ferr, double *berr, slu_memusage_t *mem_usage, SuperLUStat_t *stat, int *info);
166  void create_dense_matrix(SuperMatrix *X, int m, int n, typename SuperLuType<Scalar>::Scalar *x, int ldx, Stype_t stype, Dtype_t dtype, Mtype_t mtype);
167 #endif //SLU_MT
168
169 #ifndef SLU_MT
170  char equed[1];
172 #else
173  equed_t equed;
176  SuperMatrix AC;
177 #endif //SLU_MT
178  template<typename T> friend LinearMatrixSolver<T>* create_linear_solver(Matrix<T>* matrix, Vector<T>* rhs, bool use_direct_solver = false);
179  };
180  }
181 }
182 #endif //SUPER_LU
183 #endif
General namespace for the Hermes library.
Linear matrix solver functionality.
Definition: compat.h:85
Basic cs (Compressed sparse) matrix classes and operations.
HERMES_API LinearMatrixSolver< Scalar > * create_linear_solver(Matrix< Scalar > *matrix, Vector< Scalar > *rhs, bool use_direct_solver=false)
Function returning a solver according to the users's choice.