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