Dune Core Modules (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>
20#include <dune/istl/solverfactory.hh>
21
22namespace 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>
42 struct SuperLUSolveChooser
43 {};
44
45 template<class T>
46 struct SuperLUDenseMatChooser
47 {};
48
49 template<class T>
50 struct SuperLUQueryChooser
51 {};
52
53 template<class T>
54 struct QuerySpaceChooser
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;
276 typedef Dune::SuperLUMatrix<Matrix> SuperLUMatrix;
278 typedef SuperMatrixInitializer<Matrix> MatrixInitializer;
280 using domain_type = typename Impl::SuperLUVectorChooser<M>::domain_type;
282 using range_type = typename Impl::SuperLUVectorChooser<M>::range_type;
283
285 virtual SolverCategory::Category category() const
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
335 void apply(domain_type& x, range_type& b, InverseOperatorResult& res);
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;
373 friend struct SeqOverlappingSchwarzAssemblerHelper<SuperLU<Matrix>,true>;
374
375 SuperLUMatrix& getInternalMatrix() { return mat; }
376
378 void decompose();
379
380 SuperLUMatrix mat;
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>
394 ::~SuperLU()
395 {
396 if(mat.N()+mat.M()>0)
397 free();
398 }
399
400 template<typename M>
401 void SuperLU<M>::free()
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>
430 SuperLU<M>::SuperLU()
431 : work(0), lwork(0),verbose(false),
432 reusevector(false)
433 {}
434 template<typename M>
435 void SuperLU<M>::setVerbosity(bool v)
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>
455 void SuperLU<M>::setSubMatrix(const Matrix& mat_,
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>
561 ::apply(domain_type& x, range_type& b, InverseOperatorResult& res)
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>
714 struct IsDirectSolver<SuperLU<BCRSMatrix<T,A> > >
715 {
716 enum { value=true};
717 };
718
719 template<typename T, typename A>
720 struct StoresColumnCompressed<SuperLU<BCRSMatrix<T,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
753 DUNE_REGISTER_DIRECT_SOLVER("superlu", SuperLUCreator());
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
bcrsmatrix.hh
Implementation of the BCRSMatrix class.
bvector.hh
This file implements a vector space as a tensor product of a given vector space. The number of compon...
Dune::BCRSMatrix
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:464
Dune::ISTLError
derive error class from the base class in common
Definition: istlexception.hh:17
Dune::InverseOperator
Abstract base class for all solvers.
Definition: solver.hh:97
Dune::Matrix
A generic dynamic dense matrix.
Definition: matrix.hh:559
Dune::Matrix::M
size_type M() const
Return the number of columns.
Definition: matrix.hh:698
Dune::Matrix::N
size_type N() const
Return the number of rows.
Definition: matrix.hh:693
Dune::ParameterTree
Hierarchical structure of string parameters.
Definition: parametertree.hh:35
Dune::ParameterTree::get
std::string get(const std::string &key, const std::string &defaultValue) const
get value as string
Definition: parametertree.cc:183
Dune::SuperLU
SuperLu Solver.
Definition: superlu.hh:269
Dune::SuperLU::apply
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:340
Dune::SuperLU::range_type
typename Impl::SuperLUVectorChooser< M >::range_type range_type
The type of the range of the solver.
Definition: superlu.hh:282
Dune::SuperLU::MatrixInitializer
SuperMatrixInitializer< Matrix > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:278
Dune::SuperLU::Matrix
M Matrix
The matrix type.
Definition: superlu.hh:273
Dune::SuperLU::domain_type
typename Impl::SuperLUVectorChooser< M >::domain_type domain_type
The type of the domain of the solver.
Definition: superlu.hh:280
Dune::SuperLU::category
virtual SolverCategory::Category category() const
Category of the solver (see SolverCategory::Category)
Definition: superlu.hh:285
Dune::SuperLU::SuperLUMatrix
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:276
Dune::SuperLU::SuperLU
SuperLU(const Matrix &mat, const ParameterTree &config)
Constructs the SuperLU solver.
Definition: superlu.hh:318
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::SuperLU::apply
void apply(domain_type &x, range_type &b, InverseOperatorResult &res)
Apply inverse operator,.
Definition: superlu.hh:561
Dune::SuperLU::free
void free()
free allocated space.
Definition: superlu.hh:401
Dune::SuperLU::SuperLU
SuperLU()
Empty default constructor.
Definition: superlu.hh:430
Dune::SuperLU::setMatrix
void setMatrix(const Matrix &mat)
Initialize data from given matrix.
Definition: superlu.hh:441
Dune::dinfo
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
Dune
Dune namespace.
Definition: alignedallocator.hh:11
std
STL namespace.
solvers.hh
Implementations of the inverse operator interface.
solvertype.hh
Templates characterizing the type of a solver.
stdstreams.hh
Standard Dune debug streams.
Dune::InverseOperatorResult
Statistics about the application of an inverse operator.
Definition: solver.hh:46
Dune::InverseOperatorResult::iterations
int iterations
Number of iterations.
Definition: solver.hh:65
Dune::InverseOperatorResult::converged
bool converged
True if convergence criterion has been met.
Definition: solver.hh:71
Dune::SolverCategory::Category
Category
Definition: solvercategory.hh:21
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)