dune-istl  2.2.0
superlu.hh
Go to the documentation of this file.
00001 #ifndef DUNE_SUPERLU_HH
00002 #define DUNE_SUPERLU_HH
00003 
00004 #if HAVE_SUPERLU
00005 #ifdef TRUE
00006 #undef TRUE
00007 #endif
00008 #ifdef FALSE
00009 #undef FALSE
00010 #endif
00011 #ifndef SUPERLU_NTYPE
00012 #define  SUPERLU_NTYPE 1
00013 #endif
00014 #ifdef SUPERLU_POST_2005_VERSION
00015 
00016 #if SUPERLU_NTYPE==0
00017 #include "slu_sdefs.h"
00018 #endif
00019 
00020 #if SUPERLU_NTYPE==1
00021 #include "slu_ddefs.h"
00022 #endif
00023 
00024 #if SUPERLU_NTYPE==2
00025 #include "slu_cdefs.h"
00026 #endif
00027 
00028 #if SUPERLU_NTYPE>=3
00029 #include "slu_zdefs.h"
00030 #endif
00031 
00032 #else
00033 
00034 #if SUPERLU_NTYPE==0
00035 #include "ssp_defs.h"
00036 #endif
00037 
00038 #if SUPERLU_NTYPE==1
00039 #include "dsp_defs.h"
00040 #warning Support for SuperLU older than SuperLU 3.0 from August 2005 is deprecated.
00041 #endif
00042 
00043 #if SUPERLU_NTYPE==2
00044 #include "csp_defs.h"
00045 #endif
00046 
00047 #if SUPERLU_NTYPE>=3
00048 #include "zsp_defs.h"
00049 #endif
00050 
00051 #endif
00052 #include "solvers.hh"
00053 #include "supermatrix.hh"
00054 #include<algorithm>
00055 #include<functional>
00056 #include"bcrsmatrix.hh"
00057 #include"bvector.hh"
00058 #include"istlexception.hh"
00059 #include<dune/common/fmatrix.hh>
00060 #include<dune/common/fvector.hh>
00061 #include<dune/common/stdstreams.hh>
00062 #include <dune/istl/solvertype.hh>
00063 
00064 namespace Dune
00065 {
00066 
00077   template<class Matrix>
00078   class SuperLU 
00079   {
00080   };
00081 
00082   template<class M, class T, class TM, class TD, class TA>
00083   class SeqOverlappingSchwarz;
00084   
00085   template<class T>
00086   class SeqOverlappingSchwarzAssembler;
00087 
00088   template<class T>
00089   struct SuperLUSolveChooser
00090   {
00091   };
00092 
00093   template<class T>
00094   struct SuperLUDenseMatChooser
00095   {
00096   };
00097 
00098   template<class T>
00099   struct SuperLUQueryChooser
00100   {
00101   };
00102 
00103   template<class T>
00104   struct QuerySpaceChooser
00105   {
00106   };
00107   
00108 #if SUPERLU_NTYPE==0
00109   template<>
00110   struct SuperLUDenseMatChooser<float>
00111   {
00112     static void create(SuperMatrix *mat, int n, int m, float *dat, int n1,
00113                            Stype_t stype, Dtype_t dtype, Mtype_t mtype)
00114     {
00115       sCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
00116 
00117     }
00118     
00119     static void destroy(SuperMatrix *m)
00120     {
00121     }
00122     
00123   };
00124   template<>
00125   struct SuperLUSolveChooser<float>
00126   {
00127     static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
00128                       char *equed, float *R, float *C, SuperMatrix *L, SuperMatrix *U,
00129                       void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
00130                       float *rpg, float *rcond, float *ferr, float *berr,
00131                       mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
00132     { 
00133       sgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
00134              L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
00135              memusage, stat, info);
00136     }
00137   };
00138   
00139   template<>
00140   struct QuerySpaceChooser<float>
00141   {
00142     static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
00143     {
00144       sQuerySpace(L,U,memusage);
00145     }
00146   };
00147     
00148 #endif
00149 
00150 #if SUPERLU_NTYPE==1  
00151   
00152   template<>
00153   struct SuperLUDenseMatChooser<double>
00154   {
00155     static void create(SuperMatrix *mat, int n, int m, double *dat, int n1,
00156                            Stype_t stype, Dtype_t dtype, Mtype_t mtype)
00157     {
00158       dCreate_Dense_Matrix(mat, n, m, dat, n1, stype, dtype, mtype);
00159 
00160     }
00161     
00162     static void destroy(SuperMatrix *mat)
00163     {
00164     }
00165   };
00166   template<>
00167   struct SuperLUSolveChooser<double>
00168   {
00169     static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
00170                       char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U,
00171                       void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
00172                       double *rpg, double *rcond, double *ferr, double *berr,
00173                       mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
00174     {
00175       dgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
00176              L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
00177              memusage, stat, info);
00178     }
00179   };
00180 
00181   template<>
00182   struct QuerySpaceChooser<double>
00183   {
00184     static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
00185     {
00186       dQuerySpace(L,U,memusage);
00187     }
00188   };
00189 #endif
00190 
00191 #if SUPERLU_NTYPE>=3
00192   template<>
00193   struct SuperLUDenseMatChooser<std::complex<double> >
00194   {
00195     static void create(SuperMatrix *mat, int n, int m, std::complex<double> *dat, int n1,
00196                            Stype_t stype, Dtype_t dtype, Mtype_t mtype)
00197     {
00198       zCreate_Dense_Matrix(mat, n, m, reinterpret_cast<doublecomplex*>(dat), n1, stype, dtype, mtype);
00199 
00200     }
00201     
00202     static void destroy(SuperMatrix *mat)
00203     {
00204     }
00205   };
00206 
00207   template<>
00208   struct SuperLUSolveChooser<std::complex<double> >
00209   {
00210     static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
00211                       char *equed, double *R, double *C, SuperMatrix *L, SuperMatrix *U,
00212                       void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
00213                       double *rpg, double *rcond, double *ferr, double *berr,
00214                mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
00215     {
00216       zgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
00217              L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
00218              memusage, stat, info);
00219     }
00220   };
00221   
00222   template<>
00223   struct QuerySpaceChooser<std::complex<double> >
00224   {
00225     static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
00226     {
00227       zQuerySpace(L,U,memusage);
00228     }
00229   };
00230 #endif
00231 
00232 #if SUPERLU_NTYPE==2
00233   template<>
00234   struct SuperLUDenseMatChooser<std::complex<float> >
00235   {
00236     static void create(SuperMatrix *mat, int n, int m, std::complex<float> *dat, int n1,
00237                            Stype_t stype, Dtype_t dtype, Mtype_t mtype)
00238     {
00239       cCreate_Dense_Matrix(mat, n, m, reinterpret_cast< ::complex*>(dat), n1, stype, dtype, mtype);
00240 
00241     }
00242     
00243     static void destroy(SuperMatrix *mat)
00244     {
00245     }
00246   };
00247 
00248   template<>
00249   struct SuperLUSolveChooser<std::complex<float> >
00250   {
00251     static void solve(superlu_options_t *options, SuperMatrix *mat, int *perm_c, int *perm_r, int *etree,
00252                       char *equed, float *R, float *C, SuperMatrix *L, SuperMatrix *U,
00253                       void *work, int lwork, SuperMatrix *B, SuperMatrix *X,
00254                       float *rpg, float *rcond, float *ferr, float *berr,
00255                       mem_usage_t *memusage, SuperLUStat_t *stat, int *info)
00256     {
00257       cgssvx(options, mat, perm_c, perm_r, etree, equed, R, C,
00258              L, U, work, lwork, B, X, rpg, rcond, ferr, berr,
00259              memusage, stat, info);
00260     }
00261   };
00262 
00263   template<>
00264   struct QuerySpaceChooser<std::complex<float> >
00265   {
00266     static void querySpace(SuperMatrix* L, SuperMatrix* U, mem_usage_t* memusage)
00267     {
00268       cQuerySpace(L,U,memusage);
00269     }
00270   };
00271 #endif
00272 
00286   template<typename T, typename A, int n, int m>
00287   class SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A > >
00288     : public InverseOperator<
00289         BlockVector<FieldVector<T,m>,
00290                     typename A::template rebind<FieldVector<T,m> >::other>,
00291         BlockVector<FieldVector<T,n>,
00292                     typename A::template rebind<FieldVector<T,n> >::other> >
00293   {
00294   public:
00295     /* @brief The matrix type. */
00296     typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
00297     /* @brief The corresponding SuperLU Matrix type.*/
00298     typedef Dune::SuperLUMatrix<Matrix> SuperLUMatrix;
00300     typedef Dune::BlockVector<
00301       FieldVector<T,m>,
00302       typename A::template rebind<FieldVector<T,m> >::other> domain_type;
00304     typedef Dune::BlockVector<
00305       FieldVector<T,n>,
00306       typename A::template rebind<FieldVector<T,n> >::other> range_type;
00316     explicit SuperLU(const Matrix& mat, bool verbose=false);
00322     SuperLU();
00323 
00324     ~SuperLU();
00325 
00329     void apply(domain_type& x, range_type& b, InverseOperatorResult& res);
00330 
00334     void apply (domain_type& x, range_type& b, double reduction, InverseOperatorResult& res)
00335     {
00336       apply(x,b,res);
00337     }
00338     
00342     void apply(T* x, T* b);
00343     
00345     void setMatrix(const Matrix& mat);
00346 
00347     typename SuperLUMatrix::size_type nnz() const
00348     {
00349       return mat.nnz();
00350     }
00351     
00352     template<class S>
00353     void setSubMatrix(const Matrix& mat, const S& rowIndexSet);
00354 
00355     void setVerbosity(bool v);
00356     
00361     void free();    
00362   private:
00363     friend class std::mem_fun_ref_t<void,SuperLU>;
00364     template<class M,class X, class TM, class TD, class T1>
00365     friend class SeqOverlappingSchwarz;
00366     friend class SeqOverlappingSchwarzAssembler<SuperLU<Matrix> >;
00367     
00369     void decompose();
00370     
00371     SuperLUMatrix mat;
00372     SuperMatrix L, U, B, X;
00373     int *perm_c, *perm_r, *etree;
00374     typename GetSuperLUType<T>::float_type *R, *C;
00375     T *bstore;
00376     superlu_options_t options;
00377     char equed;
00378     void *work;
00379     int lwork;
00380     bool first, verbose;    
00381   };
00382   
00383   template<typename T, typename A, int n, int m>
00384   SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
00385   ::~SuperLU()
00386   {
00387     if(mat.N()+mat.M()>0)
00388       free();
00389   }
00390 
00391   template<typename T, typename A, int n, int m>
00392   void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::free()
00393   {
00394     delete[] perm_c;
00395     delete[] perm_r;
00396     delete[] etree;
00397     delete[] R;
00398     delete[] C;
00399     if(lwork>=0){
00400       Destroy_SuperNode_Matrix(&L);
00401       Destroy_CompCol_Matrix(&U);
00402     }
00403     lwork=0;
00404     if(!first){
00405       SUPERLU_FREE(B.Store);
00406       SUPERLU_FREE(X.Store);
00407     }
00408     mat.free();
00409   }
00410 
00411   template<typename T, typename A, int n, int m>
00412   SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
00413   ::SuperLU(const Matrix& mat_, bool verbose_)
00414     : work(0), lwork(0), first(true), verbose(verbose_)
00415   {
00416     setMatrix(mat_);
00417     
00418   }
00419   template<typename T, typename A, int n, int m>
00420   SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperLU()
00421     :    work(0), lwork(0),verbose(false)
00422   {}
00423   template<typename T, typename A, int n, int m>
00424   void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setVerbosity(bool v)
00425   {
00426     verbose=v;
00427   }
00428   
00429   template<typename T, typename A, int n, int m>
00430   void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setMatrix(const Matrix& mat_)
00431   {
00432     if(mat.N()+mat.M()>0){
00433       free();
00434     }
00435     lwork=0;
00436     work=0;
00437     //a=&mat_;
00438     mat=mat_;
00439     decompose();
00440   }
00441 
00442   template<typename T, typename A, int n, int m>
00443   template<class S>
00444   void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::setSubMatrix(const Matrix& mat_,
00445                                                                 const S& mrs)
00446   {
00447     if(mat.N()+mat.M()>0){
00448       free();
00449     }
00450     lwork=0;
00451     work=0;
00452     //a=&mat_;
00453     mat.setMatrix(mat_,mrs);
00454     decompose();    
00455   }
00456 
00457   template<typename T, typename A, int n, int m>
00458   void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::decompose()
00459   {
00460     
00461     first = true;
00462     perm_c = new int[mat.M()];
00463     perm_r = new int[mat.N()];
00464     etree  = new int[mat.M()];
00465     R = new typename GetSuperLUType<T>::float_type[mat.N()];
00466     C = new typename GetSuperLUType<T>::float_type[mat.M()];
00467     
00468     set_default_options(&options);
00469     // Do the factorization
00470     B.ncol=0;
00471     B.Stype=SLU_DN;
00472     B.Dtype=GetSuperLUType<T>::type;
00473     B.Mtype= SLU_GE;
00474     DNformat fakeFormat;
00475     fakeFormat.lda=mat.N();
00476     B.Store=&fakeFormat;
00477     X.Stype=SLU_DN;
00478     X.Dtype=GetSuperLUType<T>::type;
00479     X.Mtype= SLU_GE;
00480     X.ncol=0;
00481     X.Store=&fakeFormat;
00482     
00483     typename GetSuperLUType<T>::float_type rpg, rcond, ferr, berr;
00484     int info;
00485     mem_usage_t memusage;
00486     SuperLUStat_t stat;
00487 
00488     StatInit(&stat);
00489     SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
00490                                   &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr,
00491                                   &berr, &memusage, &stat, &info);
00492 
00493     if(verbose){
00494     dinfo<<"LU factorization: dgssvx() returns info "<< info<<std::endl;
00495 
00496     if ( info == 0 || info == n+1 ) {
00497 
00498       if ( options.PivotGrowth ) 
00499         dinfo<<"Recip. pivot growth = "<<rpg<<std::endl;
00500       if ( options.ConditionNumber )
00501         dinfo<<"Recip. condition number = %e\n"<< rcond<<std::endl;
00502       SCformat* Lstore = (SCformat *) L.Store;
00503       NCformat* Ustore = (NCformat *) U.Store;
00504       dinfo<<"No of nonzeros in factor L = "<< Lstore->nnz<<std::endl;
00505       dinfo<<"No of nonzeros in factor U = "<< Ustore->nnz<<std::endl;
00506       dinfo<<"No of nonzeros in L+U = "<< Lstore->nnz + Ustore->nnz - n<<std::endl;
00507       QuerySpaceChooser<T>::querySpace(&L, &U, &memusage);
00508       dinfo<<"L\\U MB "<<memusage.for_lu/1e6<<" \ttotal MB needed "<<memusage.total_needed/1e6
00509            <<" \texpansions ";
00510       
00511 #ifdef HAVE_MEM_USAGE_T_EXPANSIONS
00512       std::cout<<memusage.expansions<<std::endl;
00513 #else
00514       std::cout<<stat.expansions<<std::endl;
00515 #endif
00516     } else if ( info > 0 && lwork == -1 ) {
00517       dinfo<<"** Estimated memory: "<< info - n<<std::endl;
00518     }
00519     if ( options.PrintStat ) StatPrint(&stat);
00520     }
00521     StatFree(&stat);
00522     /*
00523     NCformat* Ustore = (NCformat *) U.Store;
00524     int k=0;
00525     dPrint_CompCol_Matrix("U", &U);
00526     for(int i=0; i < U.ncol; ++i, ++k){
00527       std::cout<<i<<": ";
00528       for(int c=Ustore->colptr[i]; c < Ustore->colptr[i+1]; ++c)
00529         //if(Ustore->rowind[c]==i)
00530         std::cout<<Ustore->rowind[c]<<"->"<<((double*)Ustore->nzval)[c]<<" ";
00531       if(k==0){
00532         //
00533         k=-1;
00534       }std::cout<<std::endl;
00535     }
00536     dPrint_SuperNode_Matrix("L", &L);
00537     for(int i=0; i < U.ncol; ++i, ++k){
00538       std::cout<<i<<": ";
00539       for(int c=Ustore->colptr[i]; c < Ustore->colptr[i+1]; ++c)
00540         //if(Ustore->rowind[c]==i)
00541         std::cout<<Ustore->rowind[c]<<"->"<<((double*)Ustore->nzval)[c]<<" ";
00542       if(k==0){
00543         //
00544         k=-1;
00545       }std::cout<<std::endl;
00546   } */
00547   options.Fact = FACTORED;
00548   }
00549   
00550   template<typename T, typename A, int n, int m>
00551   void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
00552   ::apply(domain_type& x, range_type& b, InverseOperatorResult& res)
00553   {
00554     if(mat.M()+mat.N()==0)
00555       DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
00556     
00557     if(first){
00558       SuperLUDenseMatChooser<T>::create(&B, (int)mat.N(), 1,  reinterpret_cast<T*>(&b[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
00559       SuperLUDenseMatChooser<T>::create(&X, (int)mat.N(), 1,  reinterpret_cast<T*>(&x[0]), (int)mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
00560       first=false;
00561     }else{
00562       ((DNformat*) B.Store)->nzval=&b[0];
00563       ((DNformat*)X.Store)->nzval=&x[0];
00564     }
00565     typename GetSuperLUType<T>::float_type rpg, rcond, ferr, berr;
00566     int info;
00567     mem_usage_t memusage;
00568     SuperLUStat_t stat;
00569     /* Initialize the statistics variables. */
00570     StatInit(&stat);
00571     /*
00572     range_type d=b;
00573     a->usmv(-1, x, d);
00574     
00575     double def0=d.two_norm();
00576     */
00577 #ifdef SUPERLU_MIN_VERSION_4_3
00578     options.IterRefine=SLU_DOUBLE;
00579 #else
00580     options.IterRefine=DOUBLE;
00581 #endif
00582     
00583     SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
00584                                    &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr, &berr,
00585                                    &memusage, &stat, &info);
00586 
00587     res.iterations=1;
00588 
00589     /*
00590     if(options.Equil==YES)    
00591       // undo scaling of right hand side
00592         std::transform(reinterpret_cast<T*>(&b[0]),reinterpret_cast<T*>(&b[0])+mat.M(),
00593                        C, reinterpret_cast<T*>(&d[0]), std::divides<T>());
00594     else
00595       d=b;
00596     a->usmv(-1, x, d);
00597     res.reduction=d.two_norm()/def0;
00598     res.conv_rate = res.reduction;
00599     res.converged=(res.reduction<1e-10||d.two_norm()<1e-18);
00600     */
00601     res.converged=true;
00602     
00603     if(verbose){
00604       
00605     dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
00606     
00607     if ( info == 0 || info == n+1 ) {
00608 
00609         if ( options.IterRefine ) {
00610           std::cout<<"Iterative Refinement: steps="
00611                <<stat.RefineSteps<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
00612         }else
00613           std::cout<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
00614     } else if ( info > 0 && lwork == -1 ) {
00615       std::cout<<"** Estimated memory: "<< info - n<<" bytes"<<std::endl;
00616     }
00617     
00618     if ( options.PrintStat ) StatPrint(&stat);
00619     }
00620     StatFree(&stat);
00621   }
00622 
00623   template<typename T, typename A, int n, int m>
00624   void SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >
00625   ::apply(T* x, T* b)
00626   {
00627     if(mat.N()+mat.M()==0)
00628       DUNE_THROW(ISTLError, "Matrix of SuperLU is null!");
00629     
00630     if(first){
00631       SuperLUDenseMatChooser<T>::create(&B, mat.N(), 1,  b, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
00632       SuperLUDenseMatChooser<T>::create(&X, mat.N(), 1,  x, mat.N(), SLU_DN, GetSuperLUType<T>::type, SLU_GE);
00633       first=false;
00634     }else{
00635       ((DNformat*) B.Store)->nzval=b;
00636       ((DNformat*)X.Store)->nzval=x;
00637     }
00638     
00639     typename GetSuperLUType<T>::float_type rpg, rcond, ferr, berr;
00640     int info;
00641     mem_usage_t memusage;
00642     SuperLUStat_t stat;
00643     /* Initialize the statistics variables. */
00644     StatInit(&stat);
00645 
00646 #ifdef SUPERLU_MIN_VERSION_4_3
00647     options.IterRefine=SLU_DOUBLE;
00648 #else
00649     options.IterRefine=DOUBLE;
00650 #endif
00651     
00652     SuperLUSolveChooser<T>::solve(&options, &static_cast<SuperMatrix&>(mat), perm_c, perm_r, etree, &equed, R, C,
00653                                  &L, &U, work, lwork, &B, &X, &rpg, &rcond, &ferr, &berr,
00654                                  &memusage, &stat, &info);
00655     
00656     if(verbose){
00657       dinfo<<"Triangular solve: dgssvx() returns info "<< info<<std::endl;
00658       
00659       if ( info == 0 || info == n+1 ) {
00660         
00661         if ( options.IterRefine ) {
00662           dinfo<<"Iterative Refinement: steps="
00663                <<stat.RefineSteps<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
00664         }else
00665           dinfo<<" FERR="<<ferr<<" BERR="<<berr<<std::endl;
00666       } else if ( info > 0 && lwork == -1 ) {
00667         dinfo<<"** Estimated memory: "<< info - n<<" bytes"<<std::endl;
00668       }
00669       if ( options.PrintStat ) StatPrint(&stat);
00670     }
00671     
00672     StatFree(&stat);
00673   }
00676   template<typename T, typename A, int n, int m>
00677   struct IsDirectSolver<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
00678   {
00679     enum{ value=true};
00680   };
00681 }
00682 
00683 #endif // HAVE_SUPERLU
00684 #endif // DUNE_SUPERLU_HH