dune-istl  2.8.0
superlu.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_SUPERLU_HH
4 #define DUNE_ISTL_SUPERLU_HH
5 
6 #if HAVE_SUPERLU
7 
8 #include "superlufunctions.hh"
9 #include "solvers.hh"
10 #include "supermatrix.hh"
11 #include <algorithm>
12 #include <functional>
13 #include "bcrsmatrix.hh"
14 #include "bvector.hh"
15 #include "istlexception.hh"
16 #include <dune/common/fmatrix.hh>
17 #include <dune/common/fvector.hh>
18 #include <dune/common/stdstreams.hh>
19 #include <dune/istl/solvertype.hh>
21 
22 namespace Dune
23 {
24 
35  template<class M, class T, class TM, class TD, class TA>
36  class SeqOverlappingSchwarz;
37 
38  template<class T, bool tag>
39  struct SeqOverlappingSchwarzAssemblerHelper;
40 
41  template<class T>
43  {};
44 
45  template<class T>
47  {};
48 
49  template<class T>
51  {};
52 
53  template<class T>
55  {};
56 
57 #if __has_include("slu_sdefs.h")
58  template<>
59  struct SuperLUDenseMatChooser<float>
60  {
61  static void create(SuperMatrix *mat, int n, int m, float *dat, int n1,
62  Stype_t stype, Dtype_t dtype, Mtype_t mtype)
63  {
64  sCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
65 
66  }
67 
68  static void destroy(SuperMatrix*)
69  {}
70 
71  };
72  template<>
73  struct SuperLUSolveChooser<float>
74  {
75  static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
76  char *equed, float *R, float *C, SuperMatrix *L, SuperMatrix *U,
77  void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
78  float *rpg, float *rcond, float *ferr, float *berr,
79  mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
80  {
81  GlobalLU_t gLU;
82  sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
83  L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
84  &gLU, memusage, stat, info);
85  }
86  };
87 
88  template<>
89  struct QuerySpaceChooser<float>
90  {
91  static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
92  {
93  sQuerySpace(L,U,memusage);
94  }
95  };
96 
97 #endif
98 
99 #if __has_include("slu_ddefs.h")
100 
101  template<>
102  struct SuperLUDenseMatChooser<double>
103  {
104  static void create(SuperMatrix *mat, int n, int m, double *dat, int n1,
105  Stype_t stype, Dtype_t dtype, Mtype_t mtype)
106  {
107  dCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
108 
109  }
110 
111  static void destroy(SuperMatrix * /* mat */)
112  {}
113  };
114  template<>
115  struct SuperLUSolveChooser<double>
116  {
117  static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
118  char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U,
119  void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
120  double *rpg, double *rcond, double *ferr, double *berr,
121  mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
122  {
123  GlobalLU_t gLU;
124  dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
125  L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
126  &gLU, memusage, stat, info);
127  }
128  };
129 
130  template<>
131  struct QuerySpaceChooser<double>
132  {
133  static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
134  {
135  dQuerySpace(L,U,memusage);
136  }
137  };
138 #endif
139 
140 #if __has_include("slu_zdefs.h")
141  template<>
142  struct SuperLUDenseMatChooser<std::complex<double> >
143  {
144  static void create(SuperMatrix *mat, int n, int m, std::complex<double> *dat, int n1,
145  Stype_t stype, Dtype_t dtype, Mtype_t mtype)
146  {
147  zCreate_Dense_Matrix(mat, n, m, reinterpret_cast<doublecomplex*>(dat), n1, stype, dtype, mtype);
148 
149  }
150 
151  static void destroy(SuperMatrix*)
152  {}
153  };
154 
155  template<>
156  struct SuperLUSolveChooser<std::complex<double> >
157  {
158  static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
159  char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U,
160  void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
161  double *rpg, double *rcond, double *ferr, double *berr,
162  mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
163  {
164  GlobalLU_t gLU;
165  zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
166  L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
167  &gLU, memusage, stat, info);
168  }
169  };
170 
171  template<>
172  struct QuerySpaceChooser<std::complex<double> >
173  {
174  static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
175  {
176  zQuerySpace(L,U,memusage);
177  }
178  };
179 #endif
180 
181 #if __has_include("slu_cdefs.h")
182  template<>
183  struct SuperLUDenseMatChooser<std::complex<float> >
184  {
185  static void create(SuperMatrix *mat, int n, int m, std::complex<float> *dat, int n1,
186  Stype_t stype, Dtype_t dtype, Mtype_t mtype)
187  {
188  cCreate_Dense_Matrix(mat, n, m, reinterpret_cast< ::complex*>(dat), n1, stype, dtype, mtype);
189 
190  }
191 
192  static void destroy(SuperMatrix* /* mat */)
193  {}
194  };
195 
196  template<>
197  struct SuperLUSolveChooser<std::complex<float> >
198  {
199  static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
200  char *equed, float *R, float *C, SuperMatrix *L, SuperMatrix *U,
201  void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
202  float *rpg, float *rcond, float *ferr, float *berr,
203  mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
204  {
205  GlobalLU_t gLU;
206  cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
207  L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
208  &gLU, memusage, stat, info);
209  }
210  };
211 
212  template<>
213  struct QuerySpaceChooser<std::complex<float> >
214  {
215  static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
216  {
217  cQuerySpace(L,U,memusage);
218  }
219  };
220 #endif
221 
222  namespace Impl
223  {
224  template<class M>
225  struct SuperLUVectorChooser
226  {};
227 
228  template<typename T, typename A, int n, int m>
229  struct SuperLUVectorChooser<BCRSMatrix<FieldMatrix<T,n,m>,A > >
230  {
232  using domain_type = BlockVector<
233  FieldVector<T,m>,
234  typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,m> > >;
236  using range_type = BlockVector<
237  FieldVector<T,n>,
238  typename std::allocator_traits<A>::template rebind_alloc<FieldVector<T,n> > >;
239  };
240 
241  template<typename T, typename A>
242  struct SuperLUVectorChooser<BCRSMatrix<T,A> >
243  {
245  using domain_type = BlockVector<T, A>;
247  using range_type = BlockVector<T, A>;
248  };
249  }
250 
264  template<typename M>
265  class SuperLU
266  : public InverseOperator<
267  typename Impl::SuperLUVectorChooser<M>::domain_type,
268  typename Impl::SuperLUVectorChooser<M>::range_type >
269  {
270  using T = typename M::field_type;
271  public:
273  using Matrix = M;
274  using matrix_type = M;
280  using domain_type = typename Impl::SuperLUVectorChooser<M>::domain_type;
282  using range_type = typename Impl::SuperLUVectorChooser<M>::range_type;
283 
286  {
287  return SolverCategory::Category::sequential;
288  }
289 
304  explicit SuperLU(const Matrix& mat, bool verbose=false,
305  bool reusevector=true);
306 
307 
318  SuperLU(const Matrix& mat, const ParameterTree& config)
319  : SuperLU(mat, config.get<bool>("verbose", false), config.get<bool>("reuseVector", true))
320  {}
321 
328  SuperLU();
329 
330  ~SuperLU();
331 
336 
340  void apply (domain_type& x, range_type& b, [[maybe_unused]] double reduction, InverseOperatorResult& res)
341  {
342  apply(x,b,res);
343  }
344 
348  void apply(T* x, T* b);
349 
351  void setMatrix(const Matrix& mat);
352 
353  typename SuperLUMatrix::size_type nnz() const
354  {
355  return mat.nonzeroes();
356  }
357 
358  template<class S>
359  void setSubMatrix(const Matrix& mat, const S& rowIndexSet);
360 
361  void setVerbosity(bool v);
362 
367  void free();
368 
369  const char* name() { return "SuperLU"; }
370  private:
371  template<class Mat,class X, class TM, class TD, class T1>
372  friend class SeqOverlappingSchwarz;
374 
375  SuperLUMatrix& getInternalMatrix() { return mat; }
376 
378  void decompose();
379 
381  SuperMatrix L, U, B, X;
382  int *perm_c, *perm_r, *etree;
383  typename GetSuperLUType<T>::float_type *R, *C;
384  T *bstore;
385  superlu_options_t options;
386  char equed;
387  void *work;
388  int lwork;
389  bool first, verbose, reusevector;
390  };
391 
392  template<typename M>
393  SuperLU<M>
395  {
396  if(mat.N()+mat.M()>0)
397  free();
398  }
399 
400  template<typename M>
402  {
403  delete[] perm_c;
404  delete[] perm_r;
405  delete[] etree;
406  delete[] R;
407  delete[] C;
408  if(lwork>=0) {
409  Destroy_SuperNode_Matrix(&L);
410  Destroy_CompCol_Matrix(&U);
411  }
412  lwork=0;
413  if(!first && reusevector) {
414  SUPERLU_FREE(B.Store);
415  SUPERLU_FREE(X.Store);
416  }
417  mat.free();
418  }
419 
420  template<typename M>
421  SuperLU<M>
422  ::SuperLU(const Matrix& mat_, bool verbose_, bool reusevector_)
423  : work(0), lwork(0), first(true), verbose(verbose_),
424  reusevector(reusevector_)
425  {
426  setMatrix(mat_);
427 
428  }
429  template<typename M>
431  : work(0), lwork(0),verbose(false),
432  reusevector(false)
433  {}
434  template<typename M>
436  {
437  verbose=v;
438  }
439 
440  template<typename M>
441  void SuperLU<M>::setMatrix(const Matrix& mat_)
442  {
443  if(mat.N()+mat.M()>0) {
444  free();
445  }
446  lwork=0;
447  work=0;
448  //a=&mat_;
449  mat=mat_;
450  decompose();
451  }
452 
453  template<typename M>
454  template<class S>
456  const S& mrs)
457  {
458  if(mat.N()+mat.M()>0) {
459  free();
460  }
461  lwork=0;
462  work=0;
463  //a=&mat_;
464  mat.setMatrix(mat_,mrs);
465  decompose();
466  }
467 
468  template<typename M>
469  void SuperLU<M>::decompose()
470  {
471 
472  first = true;
473  perm_c = new int[mat.M()];
474  perm_r = new int[mat.N()];
475  etree = new int[mat.M()];
476  R = new typename GetSuperLUType<T>::float_type[mat.N()];
477  C = new typename GetSuperLUType<T>::float_type[mat.M()];
478 
479  set_default_options(&options);
480  // Do the factorization
481  B.ncol=0;
482  B.Stype=SLU_DN;
483  B.Dtype=GetSuperLUType<T>::type;
484  B.Mtype= SLU_GE;
485  DNformat fakeFormat;
486  fakeFormat.lda=mat.N();
487  B.Store=&fakeFormat;
488  X.Stype=SLU_DN;
489  X.Dtype=GetSuperLUType<T>::type;
490  X.Mtype= SLU_GE;
491  X.ncol=0;
492  X.Store=&fakeFormat;
493 
494  typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr=1e10;
495  int info;
496  mem_usage_t memusage;
497  SuperLUStat_t stat;
498 
499  StatInit(&stat);
500  SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
501  &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
502  &berr, &memusage, &stat, &info);
503 
504  if(verbose) {
505  dinfo<<"LU factorization: dgssvx() returns info "<< info<<std::endl;
506 
507  auto nSuperLUCol = static_cast<SuperMatrix&>(mat).ncol;
508 
509  if ( info == 0 || info == nSuperLUCol+1 ) {
510 
511  if ( options.PivotGrowth )
512  dinfo<<"Recip. pivot growth = "<<rpg<<std::endl;
513  if ( options.ConditionNumber )
514  dinfo<<"Recip. condition number = %e\n"<< rcond<<std::endl;
515  SCformat* Lstore = (SCformat *) L.Store;
516  NCformat* Ustore = (NCformat *) U.Store;
517  dinfo<<"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
518  dinfo<<"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
519  dinfo<<"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - nSuperLUCol<<std::endl;
520  QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
521  dinfo<<"L\\U MB "<<memusage.for_lu/1e6<<" \ttotal MB needed "<<memusage.total_needed/1e6
522  <<" \texpansions ";
523  std::cout<<stat.expansions<<std::endl;
524 
525  } else if ( info > 0 && lwork == -1 ) { // Memory allocation failed
526  dinfo<<"** Estimated memory: "<< info - nSuperLUCol<<std::endl;
527  }
528  if ( options.PrintStat ) StatPrint(&stat);
529  }
530  StatFree(&stat);
531  /*
532  NCformat* Ustore = (NCformat *) U.Store;
533  int k=0;
534  dPrint_CompCol_Matrix("U", &U);
535  for(int i=0; i < U.ncol; ++i, ++k){
536  std::cout<<i<<": ";
537  for(int c=Ustore->colptr[i]; c < Ustore->colptr[i+1]; ++c)
538  //if(Ustore->rowind[c]==i)
539  std::cout<<Ustore->rowind[c]<<"->"<<((double*)Ustore->nzval)[c]<<" ";
540  if(k==0){
541  //
542  k=-1;
543  }std::cout<<std::endl;
544  }
545  dPrint_SuperNode_Matrix("L", &L);
546  for(int i=0; i < U.ncol; ++i, ++k){
547  std::cout<<i<<": ";
548  for(int c=Ustore->colptr[i]; c < Ustore->colptr[i+1]; ++c)
549  //if(Ustore->rowind[c]==i)
550  std::cout<<Ustore->rowind[c]<<"->"<<((double*)Ustore->nzval)[c]<<" ";
551  if(k==0){
552  //
553  k=-1;
554  }std::cout<<std::endl;
555  } */
556  options.Fact = FACTORED;
557  }
558 
559  template<typename M>
560  void SuperLU<M>
562  {
563  if (mat.N() != b.dim())
564  DUNE_THROW(ISTLError, "Size of right-hand-side vector b does not match the number of matrix rows!");
565  if (mat.M() != x.dim())
566  DUNE_THROW(ISTLError, "Size of solution vector x does not match the number of matrix columns!");
567  if (mat.M()+mat.N()==0)
568  DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
569 
570  SuperMatrix* mB = &B;
571  SuperMatrix* mX = &X;
572  SuperMatrix rB, rX;
573  if (reusevector) {
574  if(first) {
575  SuperLUDenseMatChooser<T>::create(&B, (int)mat.N(), 1, reinterpret_cast<T*>(&b[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
576  SuperLUDenseMatChooser<T>::create(&X, (int)mat.N(), 1, reinterpret_cast<T*>(&x[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
577  first=false;
578  }else{
579  ((DNformat*)B.Store)->nzval=&b[0];
580  ((DNformat*)X.Store)->nzval=&x[0];
581  }
582  } else {
583  SuperLUDenseMatChooser<T>::create(&rB, (int)mat.N(), 1, reinterpret_cast<T*>(&b[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
584  SuperLUDenseMatChooser<T>::create(&rX, (int)mat.N(), 1, reinterpret_cast<T*>(&x[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
585  mB = &rB;
586  mX = &rX;
587  }
588  typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
589  int info;
590  mem_usage_t memusage;
591  SuperLUStat_t stat;
592  /* Initialize the statistics variables. */
593  StatInit(&stat);
594  /*
595  range_type d=b;
596  a->usmv(-1, x, d);
597 
598  double def0=d.two_norm();
599  */
600  options.IterRefine=SLU_DOUBLE;
601 
602  SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
603  &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
604  &memusage, &stat, &info);
605 
606  res.iterations=1;
607 
608  /*
609  if(options.Equil==YES)
610  // undo scaling of right hand side
611  std::transform(reinterpret_cast<T*>(&b[0]),reinterpret_cast<T*>(&b[0])+mat.M(),
612  C, reinterpret_cast<T*>(&d[0]), std::divides<T>());
613  else
614  d=b;
615  a->usmv(-1, x, d);
616  res.reduction=d.two_norm()/def0;
617  res.conv_rate = res.reduction;
618  res.converged=(res.reduction<1e-10||d.two_norm()<1e-18);
619  */
620  res.converged=true;
621 
622  if(verbose) {
623 
624  dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
625 
626  auto nSuperLUCol = static_cast<SuperMatrix&>(mat).ncol;
627 
628  if ( info == 0 || info == nSuperLUCol+1 ) {
629 
630  if ( options.IterRefine ) {
631  std::cout<<"Iterative Refinement: steps="
632  <<stat.RefineSteps<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
633  }else
634  std::cout<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
635  } else if ( info > 0 && lwork == -1 ) { // Memory allocation failed
636  std::cout<<"** Estimated memory: "<< info - nSuperLUCol<<" bytes"<<std::endl;
637  }
638 
639  if ( options.PrintStat ) StatPrint(&stat);
640  }
641  StatFree(&stat);
642  if (!reusevector) {
643  SUPERLU_FREE(rB.Store);
644  SUPERLU_FREE(rX.Store);
645  }
646  }
647 
648  template<typename M>
649  void SuperLU<M>
650  ::apply(T* x, T* b)
651  {
652  if(mat.N()+mat.M()==0)
653  DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
654 
655  SuperMatrix* mB = &B;
656  SuperMatrix* mX = &X;
657  SuperMatrix rB, rX;
658  if (reusevector) {
659  if(first) {
660  SuperLUDenseMatChooser<T>::create(&B, mat.N(), 1, b, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
661  SuperLUDenseMatChooser<T>::create(&X, mat.N(), 1, x, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
662  first=false;
663  }else{
664  ((DNformat*) B.Store)->nzval=b;
665  ((DNformat*)X.Store)->nzval=x;
666  }
667  } else {
668  SuperLUDenseMatChooser<T>::create(&rB, mat.N(), 1, b, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
669  SuperLUDenseMatChooser<T>::create(&rX, mat.N(), 1, x, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
670  mB = &rB;
671  mX = &rX;
672  }
673 
674  typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
675  int info;
676  mem_usage_t memusage;
677  SuperLUStat_t stat;
678  /* Initialize the statistics variables. */
679  StatInit(&stat);
680 
681  options.IterRefine=SLU_DOUBLE;
682 
683  SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
684  &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
685  &memusage, &stat, &info);
686 
687  if(verbose) {
688  dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
689 
690  auto nSuperLUCol = static_cast<SuperMatrix&>(mat).ncol;
691 
692  if ( info == 0 || info == nSuperLUCol+1 ) { // Factorization has succeeded
693 
694  if ( options.IterRefine ) {
695  dinfo<<"Iterative Refinement: steps="
696  <<stat.RefineSteps<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
697  }else
698  dinfo<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
699  } else if ( info > 0 && lwork == -1 ) { // Memory allocation failed
700  dinfo<<"** Estimated memory: "<< info - nSuperLUCol<<" bytes"<<std::endl;
701  }
702  if ( options.PrintStat ) StatPrint(&stat);
703  }
704 
705  StatFree(&stat);
706  if (!reusevector) {
707  SUPERLU_FREE(rB.Store);
708  SUPERLU_FREE(rX.Store);
709  }
710  }
713  template<typename T, typename A>
715  {
716  enum { value=true};
717  };
718 
719  template<typename T, typename A>
721  {
722  enum { value = true };
723  };
724 
725  struct SuperLUCreator {
726  template<class> struct isValidBlock : std::false_type{};
727  template<int k> struct isValidBlock<Dune::FieldVector<double,k>> : std::true_type{};
728  template<int k> struct isValidBlock<Dune::FieldVector<std::complex<double>,k>> : std::true_type{};
729  template<typename TL, typename M>
730  std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
731  typename Dune::TypeListElement<2, TL>::type>>
732  operator() (TL /*tl*/, const M& mat, const Dune::ParameterTree& config,
733  std::enable_if_t<isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
734  {
735  int verbose = config.get("verbose", 0);
736  return std::make_shared<Dune::SuperLU<M>>(mat,verbose);
737  }
738 
739  // second version with SFINAE to validate the template parameters of SuperLU
740  template<typename TL, typename M>
741  std::shared_ptr<Dune::InverseOperator<typename Dune::TypeListElement<1, TL>::type,
742  typename Dune::TypeListElement<2, TL>::type>>
743  operator() (TL /*tl*/, const M& /*mat*/, const Dune::ParameterTree& /*config*/,
744  std::enable_if_t<!isValidBlock<typename Dune::TypeListElement<1, TL>::type::block_type>::value,int> = 0) const
745  {
746  DUNE_THROW(UnsupportedType,
747  "Unsupported Type in SuperLU (only double and std::complex<double> supported)");
748  }
749  };
750  template<> struct SuperLUCreator::isValidBlock<double> : std::true_type{};
751  template<> struct SuperLUCreator::isValidBlock<std::complex<double>> : std::true_type{};
752 
754 } // end namespace DUNE
755 
756 // undefine macros from SuperLU's slu_util.h
757 #undef FIRSTCOL_OF_SNODE
758 #undef NO_MARKER
759 #undef NUM_TEMPV
760 #undef USER_ABORT
761 #undef USER_MALLOC
762 #undef SUPERLU_MALLOC
763 #undef USER_FREE
764 #undef SUPERLU_FREE
765 #undef CHECK_MALLOC
766 #undef SUPERLU_MAX
767 #undef SUPERLU_MIN
768 #undef L_SUB_START
769 #undef L_SUB
770 #undef L_NZ_START
771 #undef L_FST_SUPC
772 #undef U_NZ_START
773 #undef U_SUB
774 #undef TRUE
775 #undef FALSE
776 #undef EMPTY
777 #undef NODROP
778 #undef DROP_BASIC
779 #undef DROP_PROWS
780 #undef DROP_COLUMN
781 #undef DROP_AREA
782 #undef DROP_SECONDARY
783 #undef DROP_DYNAMIC
784 #undef DROP_INTERP
785 #undef MILU_ALPHA
786 
787 #endif // HAVE_SUPERLU
788 #endif // DUNE_SUPERLU_HH
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Implementations of the inverse operator interface.
Templates characterizing the type of a solver.
void setSubMatrix(const Matrix &mat, const S &rowIndexSet)
Definition: superlu.hh:455
void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: superlu.hh:561
DUNE_REGISTER_DIRECT_SOLVER("ldl", Dune::LDLCreator())
void setVerbosity(bool v)
Definition: superlu.hh:435
void free()
free allocated space.
Definition: superlu.hh:401
~SuperLU()
Definition: superlu.hh:394
SuperLU()
Empty default constructor.
Definition: superlu.hh:430
void setMatrix(const Matrix &mat)
Initialize data from given matrix.
Definition: superlu.hh:441
Matrix & mat
Definition: matrixmatrix.hh:345
Definition: allocator.hh:9
PropertyMapTypeSelector< Amg::VertexVisitedTag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > >::Type get([[maybe_unused]] const Amg::VertexVisitedTag &tag, Amg::PropertiesGraph< G, Amg::VertexProperties, EP, VM, EM > &graph)
Definition: dependency.hh:291
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
derive error class from the base class in common
Definition: istlexception.hh:17
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:753
Definition: overlappingschwarz.hh:692
A generic dynamic dense matrix.
Definition: matrix.hh:559
size_type M() const
Return the number of columns.
Definition: matrix.hh:698
size_type N() const
Return the number of rows.
Definition: matrix.hh:693
Statistics about the application of an inverse operator.
Definition: solver.hh:46
int iterations
Number of iterations.
Definition: solver.hh:65
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Abstract base class for all solvers.
Definition: solver.hh:97
Category
Definition: solvercategory.hh:21
Definition: solverregistry.hh:75
Definition: solvertype.hh:14
@ value
Whether this is a direct solver.
Definition: solvertype.hh:22
Definition: solvertype.hh:28
@ value
whether the solver internally uses column compressed storage
Definition: solvertype.hh:34
Definition: superlu.hh:43
Definition: superlu.hh:47
Definition: superlu.hh:51
Definition: superlu.hh:55
SuperLu Solver.
Definition: superlu.hh:269
void apply(domain_type &x, range_type &b, [[maybe_unused]] double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:340
SuperLUMatrix::size_type nnz() const
Definition: superlu.hh:353
typename Impl::SuperLUVectorChooser< M >::range_type range_type
The type of the range of the solver.
Definition: superlu.hh:282
M matrix_type
Definition: superlu.hh:274
SuperMatrixInitializer< Matrix > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:278
M Matrix
The matrix type.
Definition: superlu.hh:273
typename Impl::SuperLUVectorChooser< M >::domain_type domain_type
The type of the domain of the solver.
Definition: superlu.hh:280
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: superlu.hh:285
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:276
SuperLU(const Matrix &mat, const ParameterTree &config)
Constructs the SuperLU solver.
Definition: superlu.hh:318
const char * name()
Definition: superlu.hh:369
Definition: superlu.hh:725
std::shared_ptr< Dune::InverseOperator< typename Dune::TypeListElement< 1, TL >::type, typename Dune::TypeListElement< 2, TL >::type > > operator()(TL, const M &mat, const Dune::ParameterTree &config, std::enable_if_t< isValidBlock< typename Dune::TypeListElement< 1, TL >::type::block_type >::value, int >=0) const
Definition: superlu.hh:732
Definition: superlu.hh:726
Definition: supermatrix.hh:130
Definition: supermatrix.hh:177