Dune Core Modules (2.4.1)

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 sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
130 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
131 memusage, stat, info);
132 }
133 };
134
135 template<>
136 struct QuerySpaceChooser<float>
137 {
138 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
139 {
140 sQuerySpace(L,U,memusage);
141 }
142 };
143
144#endif
145
146#if SUPERLU_NTYPE==1
147
148 template<>
149 struct SuperLUDenseMatChooser<double>
150 {
151 static void create(SuperMatrix *mat, int n, int m, double *dat, int n1,
152 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
153 {
154 dCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
155
156 }
157
158 static void destroy(SuperMatrix * /* mat */)
159 {}
160 };
161 template<>
162 struct SuperLUSolveChooser<double>
163 {
164 static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
165 char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U,
166 void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
167 double *rpg, double *rcond, double *ferr, double *berr,
168 mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
169 {
170 dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
171 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
172 memusage, stat, info);
173 }
174 };
175
176 template<>
177 struct QuerySpaceChooser<double>
178 {
179 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
180 {
181 dQuerySpace(L,U,memusage);
182 }
183 };
184#endif
185
186#if SUPERLU_NTYPE>=3
187 template<>
188 struct SuperLUDenseMatChooser<std::complex<double> >
189 {
190 static void create(SuperMatrix *mat, int n, int m, std::complex<double> *dat, int n1,
191 Stype_t stype, Dtype_t dtype, Mtype_t mtype)
192 {
193 zCreate_Dense_Matrix(mat, n, m, reinterpret_cast<doublecomplex*>(dat), n1, stype, dtype, mtype);
194
195 }
196
197 static void destroy(SuperMatrix*)
198 {}
199 };
200
201 template<>
202 struct SuperLUSolveChooser<std::complex<double> >
203 {
204 static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
205 char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U,
206 void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
207 double *rpg, double *rcond, double *ferr, double *berr,
208 mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
209 {
210 zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
211 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
212 memusage, stat, info);
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 cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
251 L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
252 memusage, stat, info);
253 }
254 };
255
256 template<>
257 struct QuerySpaceChooser<std::complex<float> >
258 {
259 static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
260 {
261 cQuerySpace(L,U,memusage);
262 }
263 };
264#endif
265
279 template<typename T, typename A, int n, int m>
280 class SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A > >
281 : public InverseOperator<
282 BlockVector<FieldVector<T,m>,
283 typename A::template rebind<FieldVector<T,m> >::other>,
284 BlockVector<FieldVector<T,n>,
285 typename A::template rebind<FieldVector<T,n> >::other> >
286 {
287 public:
289 typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
290 typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> matrix_type;
292 typedef Dune::SuperLUMatrix<Matrix> SuperLUMatrix;
294 typedef SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> > MatrixInitializer;
296 typedef Dune::BlockVector<
297 FieldVector<T,m>,
298 typename A::template rebind<FieldVector<T,m> >::other> domain_type;
300 typedef Dune::BlockVector<
301 FieldVector<T,n>,
302 typename A::template rebind<FieldVector<T,n> >::other> range_type;
317 explicit SuperLU(const Matrix& mat, bool verbose=false,
318 bool reusevector=true);
325 SuperLU();
326
327 ~SuperLU();
328
332 void apply(domain_type& x, range_type& b, InverseOperatorResult& res);
333
337 void apply (domain_type& x, range_type& b, double reduction, InverseOperatorResult& res)
338 {
339 DUNE_UNUSED_PARAMETER(reduction);
340 apply(x,b,res);
341 }
342
346 void apply(T* x, T* b);
347
349 void setMatrix(const Matrix& mat);
350
351 typename SuperLUMatrix::size_type nnz() const
352 {
353 return mat.nnz();
354 }
355
356 template<class S>
357 void setSubMatrix(const Matrix& mat, const S& rowIndexSet);
358
359 void setVerbosity(bool v);
360
365 void free();
366
367 const char* name() { return "SuperLU"; }
368 private:
369 friend class std::mem_fun_ref_t<void,SuperLU>;
370 template<class M,class X, class TM, class TD, class T1>
371 friend class SeqOverlappingSchwarz;
372 friend struct SeqOverlappingSchwarzAssemblerHelper<SuperLU<Matrix>,true>;
373
374 SuperLUMatrix& getInternalMatrix() { return mat; }
375
377 void decompose();
378
379 SuperLUMatrix mat;
380 SuperMatrix L, U, B, X;
381 int *perm_c, *perm_r, *etree;
382 typename GetSuperLUType<T>::float_type *R, *C;
383 T *bstore;
384 superlu_options_t options;
385 char equed;
386 void *work;
387 int lwork;
388 bool first, verbose, reusevector;
389 };
390
391 template<typename T, typename A, int n, int m>
392 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
393 ::~SuperLU()
394 {
395 if(mat.N()+mat.M()>0)
396 free();
397 }
398
399 template<typename T, typename A, int n, int m>
400 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::free()
401 {
402 delete[] perm_c;
403 delete[] perm_r;
404 delete[] etree;
405 delete[] R;
406 delete[] C;
407 if(lwork>=0) {
408 Destroy_SuperNode_Matrix(&L);
409 Destroy_CompCol_Matrix(&U);
410 }
411 lwork=0;
412 if(!first && reusevector) {
413 SUPERLU_FREE(B.Store);
414 SUPERLU_FREE(X.Store);
415 }
416 mat.free();
417 }
418
419 template<typename T, typename A, int n, int m>
420 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
421 ::SuperLU(const Matrix& mat_, bool verbose_, bool reusevector_)
422 : work(0), lwork(0), first(true), verbose(verbose_),
423 reusevector(reusevector_)
424 {
425 setMatrix(mat_);
426
427 }
428 template<typename T, typename A, int n, int m>
429 SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperLU()
430 : work(0), lwork(0),verbose(false),
431 reusevector(false)
432 {}
433 template<typename T, typename A, int n, int m>
434 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setVerbosity(bool v)
435 {
436 verbose=v;
437 }
438
439 template<typename T, typename A, int n, int m>
440 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setMatrix(const Matrix& mat_)
441 {
442 if(mat.N()+mat.M()>0) {
443 free();
444 }
445 lwork=0;
446 work=0;
447 //a=&mat_;
448 mat=mat_;
449 decompose();
450 }
451
452 template<typename T, typename A, int n, int m>
453 template<class S>
454 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(const Matrix& mat_,
455 const S& mrs)
456 {
457 if(mat.N()+mat.M()>0) {
458 free();
459 }
460 lwork=0;
461 work=0;
462 //a=&mat_;
463 mat.setMatrix(mat_,mrs);
464 decompose();
465 }
466
467 template<typename T, typename A, int n, int m>
468 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::decompose()
469 {
470
471 first = true;
472 perm_c = new int[mat.M()];
473 perm_r = new int[mat.N()];
474 etree = new int[mat.M()];
475 R = new typename GetSuperLUType<T>::float_type[mat.N()];
476 C = new typename GetSuperLUType<T>::float_type[mat.M()];
477
478 set_default_options(&options);
479 // Do the factorization
480 B.ncol=0;
481 B.Stype=SLU_DN;
482 B.Dtype=GetSuperLUType<T>::type;
483 B.Mtype= SLU_GE;
484 DNformat fakeFormat;
485 fakeFormat.lda=mat.N();
486 B.Store=&fakeFormat;
487 X.Stype=SLU_DN;
488 X.Dtype=GetSuperLUType<T>::type;
489 X.Mtype= SLU_GE;
490 X.ncol=0;
491 X.Store=&fakeFormat;
492
493 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr=1e10;
494 int info;
495 mem_usage_t memusage;
496 SuperLUStat_t stat;
497
498 StatInit(&stat);
499 SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
500 &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
501 &berr, &memusage, &stat, &info);
502
503 if(verbose) {
504 dinfo<<"LU factorization: dgssvx() returns info "<< info<<std::endl;
505
506 if ( info == 0 || info == n+1 ) {
507
508 if ( options.PivotGrowth )
509 dinfo<<"Recip. pivot growth = "<<rpg<<std::endl;
510 if ( options.ConditionNumber )
511 dinfo<<"Recip. condition number = %e\n"<< rcond<<std::endl;
512 SCformat* Lstore = (SCformat *) L.Store;
513 NCformat* Ustore = (NCformat *) U.Store;
514 dinfo<<"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
515 dinfo<<"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
516 dinfo<<"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - n<<std::endl;
517 QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
518 dinfo<<"L\\U MB "<<memusage.for_lu/1e6<<" \ttotal MB needed "<<memusage.total_needed/1e6
519 <<" \texpansions ";
520
521#ifdef HAVE_MEM_USAGE_T_EXPANSIONS
522 std::cout<<memusage.expansions<<std::endl;
523#else
524 std::cout<<stat.expansions<<std::endl;
525#endif
526 } else if ( info > 0 && lwork == -1 ) {
527 dinfo<<"** Estimated memory: "<< info - n<<std::endl;
528 }
529 if ( options.PrintStat ) StatPrint(&stat);
530 }
531 StatFree(&stat);
532 /*
533 NCformat* Ustore = (NCformat *) U.Store;
534 int k=0;
535 dPrint_CompCol_Matrix("U", &U);
536 for(int i=0; i < U.ncol; ++i, ++k){
537 std::cout<<i<<": ";
538 for(int c=Ustore->colptr[i]; c < Ustore->colptr[i+1]; ++c)
539 //if(Ustore->rowind[c]==i)
540 std::cout<<Ustore->rowind[c]<<"->"<<((double*)Ustore->nzval)[c]<<" ";
541 if(k==0){
542 //
543 k=-1;
544 }std::cout<<std::endl;
545 }
546 dPrint_SuperNode_Matrix("L", &L);
547 for(int i=0; i < U.ncol; ++i, ++k){
548 std::cout<<i<<": ";
549 for(int c=Ustore->colptr[i]; c < Ustore->colptr[i+1]; ++c)
550 //if(Ustore->rowind[c]==i)
551 std::cout<<Ustore->rowind[c]<<"->"<<((double*)Ustore->nzval)[c]<<" ";
552 if(k==0){
553 //
554 k=-1;
555 }std::cout<<std::endl;
556 } */
557 options.Fact = FACTORED;
558 }
559
560 template<typename T, typename A, int n, int m>
561 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
562 ::apply(domain_type& x, range_type& b, InverseOperatorResult& res)
563 {
564 if (mat.N() != b.dim())
565 DUNE_THROW(ISTLError, "Size of right-hand-side vector b does not match the number of matrix rows!");
566 if (mat.M() != x.dim())
567 DUNE_THROW(ISTLError, "Size of solution vector x does not match the number of matrix columns!");
568 if (mat.M()+mat.N()==0)
569 DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
570
571 SuperMatrix* mB = &B;
572 SuperMatrix* mX = &X;
573 SuperMatrix rB, rX;
574 if (reusevector) {
575 if(first) {
576 SuperLUDenseMatChooser<T>::create(&B, (int)mat.N(), 1, reinterpret_cast<T*>(&b[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
577 SuperLUDenseMatChooser<T>::create(&X, (int)mat.N(), 1, reinterpret_cast<T*>(&x[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
578 first=false;
579 }else{
580 ((DNformat*)B.Store)->nzval=&b[0];
581 ((DNformat*)X.Store)->nzval=&x[0];
582 }
583 } else {
584 SuperLUDenseMatChooser<T>::create(&rB, (int)mat.N(), 1, reinterpret_cast<T*>(&b[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
585 SuperLUDenseMatChooser<T>::create(&rX, (int)mat.N(), 1, reinterpret_cast<T*>(&x[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
586 mB = &rB;
587 mX = &rX;
588 }
589 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
590 int info;
591 mem_usage_t memusage;
592 SuperLUStat_t stat;
593 /* Initialize the statistics variables. */
594 StatInit(&stat);
595 /*
596 range_type d=b;
597 a->usmv(-1, x, d);
598
599 double def0=d.two_norm();
600 */
601#ifdef SUPERLU_MIN_VERSION_4_3
602 options.IterRefine=SLU_DOUBLE;
603#else
604 options.IterRefine=DOUBLE;
605#endif
606
607 SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
608 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
609 &memusage, &stat, &info);
610
611 res.iterations=1;
612
613 /*
614 if(options.Equil==YES)
615 // undo scaling of right hand side
616 std::transform(reinterpret_cast<T*>(&b[0]),reinterpret_cast<T*>(&b[0])+mat.M(),
617 C, reinterpret_cast<T*>(&d[0]), std::divides<T>());
618 else
619 d=b;
620 a->usmv(-1, x, d);
621 res.reduction=d.two_norm()/def0;
622 res.conv_rate = res.reduction;
623 res.converged=(res.reduction<1e-10||d.two_norm()<1e-18);
624 */
625 res.converged=true;
626
627 if(verbose) {
628
629 dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
630
631 if ( info == 0 || info == n+1 ) {
632
633 if ( options.IterRefine ) {
634 std::cout<<"Iterative Refinement: steps="
635 <<stat.RefineSteps<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
636 }else
637 std::cout<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
638 } else if ( info > 0 && lwork == -1 ) {
639 std::cout<<"** Estimated memory: "<< info - n<<" bytes"<<std::endl;
640 }
641
642 if ( options.PrintStat ) StatPrint(&stat);
643 }
644 StatFree(&stat);
645 if (!reusevector) {
646 SUPERLU_FREE(rB.Store);
647 SUPERLU_FREE(rX.Store);
648 }
649 }
650
651 template<typename T, typename A, int n, int m>
652 void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
653 ::apply(T* x, T* b)
654 {
655 if(mat.N()+mat.M()==0)
656 DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
657
658 SuperMatrix* mB = &B;
659 SuperMatrix* mX = &X;
660 SuperMatrix rB, rX;
661 if (reusevector) {
662 if(first) {
663 SuperLUDenseMatChooser<T>::create(&B, mat.N(), 1, b, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
664 SuperLUDenseMatChooser<T>::create(&X, mat.N(), 1, x, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
665 first=false;
666 }else{
667 ((DNformat*) B.Store)->nzval=b;
668 ((DNformat*)X.Store)->nzval=x;
669 }
670 } else {
671 SuperLUDenseMatChooser<T>::create(&rB, mat.N(), 1, b, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
672 SuperLUDenseMatChooser<T>::create(&rX, mat.N(), 1, x, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
673 mB = &rB;
674 mX = &rX;
675 }
676
677 typename GetSuperLUType<T>::float_type rpg, rcond, ferr=1e10, berr;
678 int info;
679 mem_usage_t memusage;
680 SuperLUStat_t stat;
681 /* Initialize the statistics variables. */
682 StatInit(&stat);
683
684#ifdef SUPERLU_MIN_VERSION_4_3
685 options.IterRefine=SLU_DOUBLE;
686#else
687 options.IterRefine=DOUBLE;
688#endif
689
690 SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
691 &L, &U, work, lwork, mB, mX, &rpg, &rcond, &ferr, &berr,
692 &memusage, &stat, &info);
693
694 if(verbose) {
695 dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
696
697 if ( info == 0 || info == n+1 ) {
698
699 if ( options.IterRefine ) {
700 dinfo<<"Iterative Refinement: steps="
701 <<stat.RefineSteps<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
702 }else
703 dinfo<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
704 } else if ( info > 0 && lwork == -1 ) {
705 dinfo<<"** Estimated memory: "<< info - n<<" bytes"<<std::endl;
706 }
707 if ( options.PrintStat ) StatPrint(&stat);
708 }
709
710 StatFree(&stat);
711 if (!reusevector) {
712 SUPERLU_FREE(rB.Store);
713 SUPERLU_FREE(rX.Store);
714 }
715 }
718 template<typename T, typename A, int n, int m>
719 struct IsDirectSolver<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
720 {
721 enum { value=true};
722 };
723
724 template<typename T, typename A, int n, int m>
725 struct StoresColumnCompressed<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
726 {
727 enum { value = true };
728 };
729}
730
731// undefine macros from SuperLU's slu_util.h
732#undef FIRSTCOL_OF_SNODE
733#undef NO_MARKER
734#undef NUM_TEMPV
735#undef USER_ABORT
736#undef USER_MALLOC
737#undef SUPERLU_MALLOC
738#undef USER_FREE
739#undef SUPERLU_FREE
740#undef CHECK_MALLOC
741#undef SUPERLU_MAX
742#undef SUPERLU_MIN
743#undef L_SUB_START
744#undef L_SUB
745#undef L_NZ_START
746#undef L_FST_SUPC
747#undef U_NZ_START
748#undef U_SUB
749#undef TRUE
750#undef FALSE
751#undef EMPTY
752#undef NODROP
753#undef DROP_BASIC
754#undef DROP_PROWS
755#undef DROP_COLUMN
756#undef DROP_AREA
757#undef DROP_SECONDARY
758#undef DROP_DYNAMIC
759#undef DROP_INTERP
760#undef MILU_ALPHA
761
762#endif // HAVE_SUPERLU
763#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:302
Dune::SuperLU< BCRSMatrix< FieldMatrix< T, n, m >, A > >::SuperLUMatrix
Dune::SuperLUMatrix< Matrix > SuperLUMatrix
The corresponding SuperLU Matrix type.
Definition: superlu.hh:292
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:298
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:337
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:294
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 24, 22:29, 2024)