DUNE-FEM (unstable)

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