Dune Core Modules (2.4.2)

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#ifdef TRUE
8#undef TRUE
9#endif
10#ifdef FALSE
11#undef FALSE
12#endif
13#ifndef SUPERLU_NTYPE
14#define SUPERLU_NTYPE 1
15#endif
16#ifdef SUPERLU_POST_2005_VERSION
17
18#if SUPERLU_NTYPE==0
19#include "slu_sdefs.h"
20#endif
21
22#if SUPERLU_NTYPE==1
23#include "slu_ddefs.h"
24#endif
25
26#if SUPERLU_NTYPE==2
27#include "slu_cdefs.h"
28#endif
29
30#if SUPERLU_NTYPE>=3
31#include "slu_zdefs.h"
32#endif
33
34#else
35
36#if SUPERLU_NTYPE==0
37#include "ssp_defs.h"
38#endif
39
40#if SUPERLU_NTYPE==1
41#include "dsp_defs.h"
42#warning Support for SuperLU older than SuperLU 3.0 from August 2005 is deprecated.
43#endif
44
45#if SUPERLU_NTYPE==2
46#include "csp_defs.h"
47#endif
48
49#if SUPERLU_NTYPE>=3
50#include "zsp_defs.h"
51#endif
52
53#endif
54#include "solvers.hh"
55#include "supermatrix.hh"
56#include <algorithm>
57#include <functional>
58#include "bcrsmatrix.hh"
59#include "bvector.hh"
60#include "istlexception.hh"
61#include <dune/common/fmatrix.hh>
62#include <dune/common/fvector.hh>
63#include <dune/common/stdstreams.hh>
64#include <dune/istl/solvertype.hh>
65
66namespace Dune
67{
68
79 template<class Matrix>
80 class SuperLU
81 {};
82
83 template<class M, class T, class TM, class TD, class TA>
84 class SeqOverlappingSchwarz;
85
86 template<class T, bool tag>
87 struct SeqOverlappingSchwarzAssemblerHelper;
88
89 template<class T>
90 struct SuperLUSolveChooser
91 {};
92
93 template<class T>
94 struct SuperLUDenseMatChooser
95 {};
96
97 template<class T>
98 struct SuperLUQueryChooser
99 {};
100
101 template<class T>
102 struct QuerySpaceChooser
103 {};
104
105#if SUPERLU_NTYPE==0
106 template<>
107 struct SuperLUDenseMatChooser<float>
108 {
109 static void create(SuperMatrix *mat, int n, int m, float *dat, int n1,
110 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
111 {
112 sCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
113
114 }
115
116 static void destroy(SuperMatrix*)
117 {}
118
119 };
120 template<>
121 struct SuperLUSolveChooser<float>
122 {
123 static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
124 char *equed, float *R, float *C, SuperMatrix *L, SuperMatrix *U,
125 void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
126 float *rpg, float *rcond, float *ferr, float *berr,
127 mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
128 {
129#if SUPERLU_MIN_VERSION_5
130 GlobalLU_t gLU;
131 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
132 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
133 &gLU, memusage, stat, info);
134#else
135 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
136 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
137 memusage, stat, info);
138#endif
139 }
140 };
141
142 template<>
143 struct QuerySpaceChooser<float>
144 {
145 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
146 {
147 sQuerySpace(L,U,memusage);
148 }
149 };
150
151#endif
152
153#if SUPERLU_NTYPE==1
154
155 template<>
156 struct SuperLUDenseMatChooser<double>
157 {
158 static void create(SuperMatrix *mat, int n, int m, double *dat, int n1,
159 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
160 {
161 dCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
162
163 }
164
165 static void destroy(SuperMatrix * /* mat */)
166 {}
167 };
168 template<>
169 struct SuperLUSolveChooser<double>
170 {
171 static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
172 char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U,
173 void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
174 double *rpg, double *rcond, double *ferr, double *berr,
175 mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
176 {
177#if SUPERLU_MIN_VERSION_5
178 GlobalLU_t gLU;
179 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
180 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
181 &gLU, memusage, stat, info);
182#else
183 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
184 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
185 memusage, stat, info);
186#endif
187 }
188 };
189
190 template<>
191 struct QuerySpaceChooser<double>
192 {
193 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
194 {
195 dQuerySpace(L,U,memusage);
196 }
197 };
198#endif
199
200#if SUPERLU_NTYPE>=3
201 template<>
202 struct SuperLUDenseMatChooser<std::complex<double> >
203 {
204 static void create(SuperMatrix *mat, int n, int m, std::complex<double> *dat, int n1,
205 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
206 {
207 zCreate_Dense_Matrix(mat, n, m, reinterpret_cast<doublecomplex*>(dat), n1, stype, dtype, mtype);
208
209 }
210
211 static void destroy(SuperMatrix*)
212 {}
213 };
214
215 template<>
216 struct SuperLUSolveChooser<std::complex<double> >
217 {
218 static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
219 char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U,
220 void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
221 double *rpg, double *rcond, double *ferr, double *berr,
222 mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
223 {
224#if SUPERLU_MIN_VERSION_5
225 GlobalLU_t gLU;
226 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
227 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
228 &gLU, memusage, stat, info);
229#else
230 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
231 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
232 memusage, stat, info);
233#endif
234 }
235 };
236
237 template<>
238 struct QuerySpaceChooser<std::complex<double> >
239 {
240 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
241 {
242 zQuerySpace(L,U,memusage);
243 }
244 };
245#endif
246
247#if SUPERLU_NTYPE==2
248 template<>
249 struct SuperLUDenseMatChooser<std::complex<float> >
250 {
251 static void create(SuperMatrix *mat, int n, int m, std::complex<float> *dat, int n1,
252 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
253 {
254 cCreate_Dense_Matrix(mat, n, m, reinterpret_cast< ::complex*>(dat), n1, stype, dtype, mtype);
255
256 }
257
258 static void destroy(SuperMatrix* /* mat */)
259 {}
260 };
261
262 template<>
263 struct SuperLUSolveChooser<std::complex<float> >
264 {
265 static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
266 char *equed, float *R, float *C, SuperMatrix *L, SuperMatrix *U,
267 void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
268 float *rpg, float *rcond, float *ferr, float *berr,
269 mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
270 {
271#if SUPERLU_MIN_VERSION_5
272 GlobalLU_t gLU;
273 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
274 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
275 &gLU, memusage, stat, info);
276#else
277 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
278 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
279 memusage, stat, info);
280#endif
281 }
282 };
283
284 template<>
285 struct QuerySpaceChooser<std::complex<float> >
286 {
287 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
288 {
289 cQuerySpace(L,U,memusage);
290 }
291 };
292#endif
293
307 template<typename T, typename A, int n, int m>
308 class SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A > >
309 : public InverseOperator<
310 BlockVector<FieldVector<T,m>,
311 typename A::template rebind<FieldVector<T,m> >::other>,
312 BlockVector<FieldVector<T,n>,
313 typename A::template rebind<FieldVector<T,n> >::other> >
314 {
315 public:
317 typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
318 typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> matrix_type;
320 typedef Dune::SuperLUMatrix<Matrix> SuperLUMatrix;
322 typedef SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> > MatrixInitializer;
324 typedef Dune::BlockVector<
325 FieldVector<T,m>,
326 typename A::template rebind<FieldVector<T,m> >::other> domain_type;
328 typedef Dune::BlockVector<
329 FieldVector<T,n>,
330 typename A::template rebind<FieldVector<T,n> >::other> range_type;
345 explicit SuperLU(const Matrix& mat, bool verbose=false,
346 bool reusevector=true);
353 SuperLU();
354
355 ~SuperLU();
356
360 void apply(domain_type& x, range_type& b, InverseOperatorResult& res);
361
365 void apply (domain_type& x, range_type& b, double reduction, InverseOperatorResult& res)
366 {
367 DUNE_UNUSED_PARAMETER(reduction);
368 apply(x,b,res);
369 }
370
374 void apply(T* x, T* b);
375
377 void setMatrix(const Matrix& mat);
378
379 typename SuperLUMatrix::size_type nnz() const
380 {
381 return mat.nnz();
382 }
383
384 template<class S>
385 void setSubMatrix(const Matrix& mat, const S& rowIndexSet);
386
387 void setVerbosity(bool v);
388
393 void free();
394
395 const char* name() { return "SuperLU"; }
396 private:
397 friend class std::mem_fun_ref_t<void,SuperLU>;
398 template<class M,class X, class TM, class TD, class T1>
399 friend class SeqOverlappingSchwarz;
400 friend struct SeqOverlappingSchwarzAssemblerHelper<SuperLU<Matrix>,true>;
401
402 SuperLUMatrix& getInternalMatrix() { return mat; }
403
405 void decompose();
406
407 SuperLUMatrix mat;
408 SuperMatrix L, U, B, X;
409 int *perm_c, *perm_r, *etree;
410 typename GetSuperLUType<T>::float_type *R, *C;
411 T *bstore;
412 superlu_options_t options;
413 char equed;
414 void *work;
415 int lwork;
416 bool first, verbose, reusevector;
417 };
418
419 template<typename T, typename A, int n, int m>
420 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
421 ::~SuperLU()
422 {
423 if(mat.N()+mat.M()>0)
424 free();
425 }
426
427 template<typename T, typename A, int n, int m>
428 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::free()
429 {
430 delete[] perm_c;
431 delete[] perm_r;
432 delete[] etree;
433 delete[] R;
434 delete[] C;
435 if(lwork>=0) {
436 Destroy_SuperNode_Matrix(&L);
437 Destroy_CompCol_Matrix(&U);
438 }
439 lwork=0;
440 if(!first && reusevector) {
441 SUPERLU_FREE(B.Store);
442 SUPERLU_FREE(X.Store);
443 }
444 mat.free();
445 }
446
447 template<typename T, typename A, int n, int m>
448 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
449 ::SuperLU(const Matrix& mat_, bool verbose_, bool reusevector_)
450 : work(0), lwork(0), first(true), verbose(verbose_),
451 reusevector(reusevector_)
452 {
453 setMatrix(mat_);
454
455 }
456 template<typename T, typename A, int n, int m>
457 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperLU()
458 : work(0), lwork(0),verbose(false),
459 reusevector(false)
460 {}
461 template<typename T, typename A, int n, int m>
462 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setVerbosity(bool v)
463 {
464 verbose=v;
465 }
466
467 template<typename T, typename A, int n, int m>
468 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setMatrix(const Matrix& mat_)
469 {
470 if(mat.N()+mat.M()>0) {
471 free();
472 }
473 lwork=0;
474 work=0;
475 //a=&mat_;
476 mat=mat_;
477 decompose();
478 }
479
480 template<typename T, typename A, int n, int m>
481 template<class S>
482 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(const Matrix& mat_,
483 const S& mrs)
484 {
485 if(mat.N()+mat.M()>0) {
486 free();
487 }
488 lwork=0;
489 work=0;
490 //a=&mat_;
491 mat.setMatrix(mat_,mrs);
492 decompose();
493 }
494
495 template<typename T, typename A, int n, int m>
496 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::decompose()
497 {
498
499 first = true;
500 perm_c = new int[mat.M()];
501 perm_r = new int[mat.N()];
502 etree = new int[mat.M()];
503 R = new typename GetSuperLUType<T>::float_type[mat.N()];
504 C = new typename GetSuperLUType<T>::float_type[mat.M()];
505
506 set_default_options(&options);
507 // Do the factorization
508 B.ncol=0;
509 B.Stype=SLU_DN;
510 B.Dtype=GetSuperLUType<T>::type;
511 B.Mtype= SLU_GE;
512 DNformat fakeFormat;
513 fakeFormat.lda=mat.N();
514 B.Store=&fakeFormat;
515 X.Stype=SLU_DN;
516 X.Dtype=GetSuperLUType<T>::type;
517 X.Mtype= SLU_GE;
518 X.ncol=0;
519 X.Store=&fakeFormat;
520
521 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr=1e10;
522 int info;
523 mem_usage_t memusage;
524 SuperLUStat_t stat;
525
526 StatInit(&stat);
527 SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
528 &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
529 &berr, &memusage, &stat, &info);
530
531 if(verbose) {
532 dinfo<<"LU factorization: dgssvx() returns info "<< info<<std::endl;
533
534 if ( info == 0 || info == n+1 ) {
535
536 if ( options.PivotGrowth )
537 dinfo<<"Recip. pivot growth = "<<rpg<<std::endl;
538 if ( options.ConditionNumber )
539 dinfo<<"Recip. condition number = %e\n"<< rcond<<std::endl;
540 SCformat* Lstore = (SCformat *) L.Store;
541 NCformat* Ustore = (NCformat *) U.Store;
542 dinfo<<"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
543 dinfo<<"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
544 dinfo<<"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - n<<std::endl;
545 QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
546 dinfo<<"L\\U MB "<<memusage.for_lu/1e6<<" \ttotal MB needed "<<memusage.total_needed/1e6
547 <<" \texpansions ";
548
549#ifdef HAVE_MEM_USAGE_T_EXPANSIONS
550 std::cout<<memusage.expansions<<std::endl;
551#else
552 std::cout<<stat.expansions<<std::endl;
553#endif
554 } else if ( info > 0 && lwork == -1 ) {
555 dinfo<<"** Estimated memory: "<< info - n<<std::endl;
556 }
557 if ( options.PrintStat ) StatPrint(&stat);
558 }
559 StatFree(&stat);
560 /*
561 NCformat* Ustore = (NCformat *) U.Store;
562 int k=0;
563 dPrint_CompCol_Matrix("U", &U);
564 for(int i=0; i < U.ncol; ++i, ++k){
565 std::cout<<i<<": ";
566 for(int c=Ustore->colptr[i]; c < Ustore->colptr[i+1]; ++c)
567 //if(Ustore->rowind[c]==i)
568 std::cout<<Ustore->rowind[c]<<"->"<<((double*)Ustore->nzval)[c]<<" ";
569 if(k==0){
570 //
571 k=-1;
572 }std::cout<<std::endl;
573 }
574 dPrint_SuperNode_Matrix("L", &L);
575 for(int i=0; i < U.ncol; ++i, ++k){
576 std::cout<<i<<": ";
577 for(int c=Ustore->colptr[i]; c < Ustore->colptr[i+1]; ++c)
578 //if(Ustore->rowind[c]==i)
579 std::cout<<Ustore->rowind[c]<<"->"<<((double*)Ustore->nzval)[c]<<" ";
580 if(k==0){
581 //
582 k=-1;
583 }std::cout<<std::endl;
584 } */
585 options.Fact = FACTORED;
586 }
587
588 template<typename T, typename A, int n, int m>
589 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
590 ::apply(domain_type& x, range_type& b, InverseOperatorResult& res)
591 {
592 if (mat.N() != b.dim())
593 DUNE_THROW(ISTLError, "Size of right-hand-side vector b does not match the number of matrix rows!");
594 if (mat.M() != x.dim())
595 DUNE_THROW(ISTLError, "Size of solution vector x does not match the number of matrix columns!");
596 if (mat.M()+mat.N()==0)
597 DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
598
599 SuperMatrix* mB = &B;
600 SuperMatrix* mX = &X;
601 SuperMatrix rB, rX;
602 if (reusevector) {
603 if(first) {
604 SuperLUDenseMatChooser<T>::create(&B, (int)mat.N(), 1, reinterpret_cast<T*>(&b[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
605 SuperLUDenseMatChooser<T>::create(&X, (int)mat.N(), 1, reinterpret_cast<T*>(&x[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
606 first=false;
607 }else{
608 ((DNformat*)B.Store)->nzval=&b[0];
609 ((DNformat*)X.Store)->nzval=&x[0];
610 }
611 } else {
612 SuperLUDenseMatChooser<T>::create(&rB, (int)mat.N(), 1, reinterpret_cast<T*>(&b[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
613 SuperLUDenseMatChooser<T>::create(&rX, (int)mat.N(), 1, reinterpret_cast<T*>(&x[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
614 mB = &rB;
615 mX = &rX;
616 }
617 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
618 int info;
619 mem_usage_t memusage;
620 SuperLUStat_t stat;
621 /* Initialize the statistics variables. */
622 StatInit(&stat);
623 /*
624 range_type d=b;
625 a->usmv(-1, x, d);
626
627 double def0=d.two_norm();
628 */
629#ifdef SUPERLU_MIN_VERSION_4_3
630 options.IterRefine=SLU_DOUBLE;
631#else
632 options.IterRefine=DOUBLE;
633#endif
634
635 SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
636 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
637 &memusage, &stat, &info);
638
639 res.iterations=1;
640
641 /*
642 if(options.Equil==YES)
643 // undo scaling of right hand side
644 std::transform(reinterpret_cast<T*>(&b[0]),reinterpret_cast<T*>(&b[0])+mat.M(),
645 C, reinterpret_cast<T*>(&d[0]), std::divides<T>());
646 else
647 d=b;
648 a->usmv(-1, x, d);
649 res.reduction=d.two_norm()/def0;
650 res.conv_rate = res.reduction;
651 res.converged=(res.reduction<1e-10||d.two_norm()<1e-18);
652 */
653 res.converged=true;
654
655 if(verbose) {
656
657 dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
658
659 if ( info == 0 || info == n+1 ) {
660
661 if ( options.IterRefine ) {
662 std::cout<<"Iterative Refinement: steps="
663 <<stat.RefineSteps<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
664 }else
665 std::cout<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
666 } else if ( info > 0 && lwork == -1 ) {
667 std::cout<<"** Estimated memory: "<< info - n<<" bytes"<<std::endl;
668 }
669
670 if ( options.PrintStat ) StatPrint(&stat);
671 }
672 StatFree(&stat);
673 if (!reusevector) {
674 SUPERLU_FREE(rB.Store);
675 SUPERLU_FREE(rX.Store);
676 }
677 }
678
679 template<typename T, typename A, int n, int m>
680 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
681 ::apply(T* x, T* b)
682 {
683 if(mat.N()+mat.M()==0)
684 DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
685
686 SuperMatrix* mB = &B;
687 SuperMatrix* mX = &X;
688 SuperMatrix rB, rX;
689 if (reusevector) {
690 if(first) {
691 SuperLUDenseMatChooser<T>::create(&B, mat.N(), 1, b, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
692 SuperLUDenseMatChooser<T>::create(&X, mat.N(), 1, x, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
693 first=false;
694 }else{
695 ((DNformat*) B.Store)->nzval=b;
696 ((DNformat*)X.Store)->nzval=x;
697 }
698 } else {
699 SuperLUDenseMatChooser<T>::create(&rB, mat.N(), 1, b, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
700 SuperLUDenseMatChooser<T>::create(&rX, mat.N(), 1, x, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
701 mB = &rB;
702 mX = &rX;
703 }
704
705 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
706 int info;
707 mem_usage_t memusage;
708 SuperLUStat_t stat;
709 /* Initialize the statistics variables. */
710 StatInit(&stat);
711
712#ifdef SUPERLU_MIN_VERSION_4_3
713 options.IterRefine=SLU_DOUBLE;
714#else
715 options.IterRefine=DOUBLE;
716#endif
717
718 SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
719 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
720 &memusage, &stat, &info);
721
722 if(verbose) {
723 dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
724
725 if ( info == 0 || info == n+1 ) {
726
727 if ( options.IterRefine ) {
728 dinfo<<"Iterative Refinement: steps="
729 <<stat.RefineSteps<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
730 }else
731 dinfo<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
732 } else if ( info > 0 && lwork == -1 ) {
733 dinfo<<"** Estimated memory: "<< info - n<<" bytes"<<std::endl;
734 }
735 if ( options.PrintStat ) StatPrint(&stat);
736 }
737
738 StatFree(&stat);
739 if (!reusevector) {
740 SUPERLU_FREE(rB.Store);
741 SUPERLU_FREE(rX.Store);
742 }
743 }
746 template<typename T, typename A, int n, int m>
747 struct IsDirectSolver<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
748 {
749 enum { value=true};
750 };
751
752 template<typename T, typename A, int n, int m>
753 struct StoresColumnCompressed<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
754 {
755 enum { value = true };
756 };
757}
758
759// undefine macros from SuperLU's slu_util.h
760#undef FIRSTCOL_OF_SNODE
761#undef NO_MARKER
762#undef NUM_TEMPV
763#undef USER_ABORT
764#undef USER_MALLOC
765#undef SUPERLU_MALLOC
766#undef USER_FREE
767#undef SUPERLU_FREE
768#undef CHECK_MALLOC
769#undef SUPERLU_MAX
770#undef SUPERLU_MIN
771#undef L_SUB_START
772#undef L_SUB
773#undef L_NZ_START
774#undef L_FST_SUPC
775#undef U_NZ_START
776#undef U_SUB
777#undef TRUE
778#undef FALSE
779#undef EMPTY
780#undef NODROP
781#undef DROP_BASIC
782#undef DROP_PROWS
783#undef DROP_COLUMN
784#undef DROP_AREA
785#undef DROP_SECONDARY
786#undef DROP_DYNAMIC
787#undef DROP_INTERP
788#undef MILU_ALPHA
789
790#endif // HAVE_SUPERLU
791#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:413
Dune::BlockVector
A vector of blocks with memory management.
Definition: bvector.hh:253
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:67
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:94
Dune::ISTLError
derive error class from the base class in common
Definition: istlexception.hh:16
Dune::InverseOperator
Abstract base class for all solvers.
Definition: solver.hh:79
Dune::Matrix
A generic dynamic dense matrix.
Definition: matrix.hh:25
Dune::Matrix::M
size_type M() const
Return the number of columns.
Definition: matrix.hh:165
Dune::Matrix::N
size_type N() const
Return the number of rows.
Definition: matrix.hh:160
Dune::SuperLU< BCRSMatrix< FieldMatrix< T, n, m >, A > >::range_type
Dune::BlockVector< FieldVector< T, n >, typename A::template rebind< FieldVector< T, n > >::other > range_type
The type of the range of the solver.
Definition: superlu.hh:330
Dune::SuperLU< BCRSMatrix< FieldMatrix< T, n, m >, A > >::SuperLUMatrix
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:320
Dune::SuperLU< BCRSMatrix< FieldMatrix< T, n, m >, A > >::domain_type
Dune::BlockVector< FieldVector< T, m >, typename A::template rebind< FieldVector< T, m > >::other > domain_type
The type of the domain of the solver.
Definition: superlu.hh:326
Dune::SuperLU< BCRSMatrix< FieldMatrix< T, n, m >, A > >::apply
void apply(domain_type &x, range_type &b, double reduction, InverseOperatorResult &res)
apply inverse operator, with given convergence criteria.
Definition: superlu.hh:365
Dune::SuperLU< BCRSMatrix< FieldMatrix< T, n, m >, A > >::MatrixInitializer
SuperMatrixInitializer< BCRSMatrix< FieldMatrix< T, n, m >, A > > MatrixInitializer
Type of an associated initializer class.
Definition: superlu.hh:322
Dune::block_vector_unmanaged::dim
size_type dim() const
dimension of the vector space
Definition: bvector.hh:223
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:243
Dune::dinfo
DInfoType dinfo(std::cout)
Stream for informative output.
Definition: stdstreams.hh:138
Dune
Dune namespace.
Definition: alignment.hh:10
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:32
Dune::InverseOperatorResult::iterations
int iterations
Number of iterations.
Definition: solver.hh:50
Dune::InverseOperatorResult::converged
bool converged
True if convergence criterion has been met.
Definition: solver.hh:56
DUNE_UNUSED_PARAMETER
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)