dune-istl  2.2.0
supermatrix.hh
Go to the documentation of this file.
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=8 sw=2 sts=2:
00003 #ifndef DUNE_SUPERLUMATRIX_HH
00004 #define DUNE_SUPERLUMATRIX_HH
00005 
00006 #if HAVE_SUPERLU
00007 #ifdef SUPERLU_POST_2005_VERSION
00008 
00009 #ifndef SUPERLU_NTYPE
00010 #define SUPERLU_NTYPE 1
00011 #endif
00012 
00013 #if SUPERLU_NTYPE==0
00014 #include "slu_sdefs.h"
00015 #endif
00016 
00017 #if SUPERLU_NTYPE==1
00018 #include "slu_ddefs.h"
00019 #endif
00020 
00021 #if SUPERLU_NTYPE==2
00022 #include "slu_cdefs.h"
00023 #endif
00024 
00025 #if SUPERLU_NTYPE>=3
00026 #include "slu_zdefs.h"
00027 #endif
00028 
00029 #else
00030 
00031 #if SUPERLU_NTYPE==0
00032 #include "ssp_defs.h"
00033 #endif
00034 
00035 #if SUPERLU_NTYPE==1
00036 #include "dsp_defs.h"
00037 #endif
00038 
00039 #if SUPERLU_NTYPE==2
00040 #include "csp_defs.h"
00041 #endif
00042 
00043 #if SUPERLU_NTYPE>=3
00044 #include "zsp_defs.h"
00045 #endif
00046 
00047 #endif
00048 #include"bcrsmatrix.hh"
00049 #include"bvector.hh"
00050 #include<dune/common/fmatrix.hh>
00051 #include<dune/common/fvector.hh>
00052 #include<limits>
00053 
00054 namespace Dune
00055 {
00056 
00057   template<class T>
00058   struct SuperMatrixCreateSparseChooser
00059   {
00060   };
00061   
00062   template<class T>
00063   struct SuperMatrixPrinter
00064   {
00065   };
00066 
00067 #if SUPERLU_NTYPE==0
00068   template<>
00069   struct SuperMatrixCreateSparseChooser<float>
00070   {
00071     static void create(SuperMatrix *mat, int n, int m, int offset,
00072                        float *values, int *rowindex, int* colindex,
00073                        Stype_t stype, Dtype_t dtype, Mtype_t mtype)
00074     {
00075       sCreate_CompCol_Matrix(mat, n, m, offset, values, rowindex, colindex,
00076                            stype, dtype, mtype);
00077     }
00078   };
00079   
00080   template<>
00081   struct SuperMatrixPrinter<float>
00082   {
00083     static void print(char* name, SuperMatrix* mat)
00084     {
00085       sPrint_CompCol_Matrix(name, mat);
00086     }
00087   };
00088 #endif
00089 
00090 #if SUPERLU_NTYPE==1
00091   template<>
00092   struct SuperMatrixCreateSparseChooser<double>
00093   {
00094     static void create(SuperMatrix *mat, int n, int m, int offset,
00095                        double *values, int *rowindex, int* colindex,
00096                        Stype_t stype, Dtype_t dtype, Mtype_t mtype)
00097     {
00098       dCreate_CompCol_Matrix(mat, n, m, offset, values, rowindex, colindex,
00099                            stype, dtype, mtype);
00100     }
00101   };
00102   
00103   template<>
00104   struct SuperMatrixPrinter<double>
00105   {
00106     static void print(char* name, SuperMatrix* mat)
00107     {
00108       dPrint_CompCol_Matrix(name, mat);
00109     }
00110   };
00111 #endif
00112 
00113 #if SUPERLU_NTYPE==2
00114   template<>
00115   struct SuperMatrixCreateSparseChooser<std::complex<float> >
00116   {
00117     static void create(SuperMatrix *mat, int n, int m, int offset,
00118                        std::complex<float> *values, int *rowindex, int* colindex,
00119                        Stype_t stype, Dtype_t dtype, Mtype_t mtype)
00120     {
00121       cCreate_CompCol_Matrix(mat, n, m, offset, reinterpret_cast< ::complex*>(values), 
00122                              rowindex, colindex, stype, dtype, mtype);
00123     }
00124   };
00125   
00126   template<>
00127   struct SuperMatrixPrinter<std::complex<float> >
00128   {
00129     static void print(char* name, SuperMatrix* mat)
00130     {
00131       cPrint_CompCol_Matrix(name, mat);
00132     }
00133   };
00134 #endif
00135 
00136 #if SUPERLU_NTYPE>=3
00137   template<>
00138   struct SuperMatrixCreateSparseChooser<std::complex<double> >
00139   {
00140     static void create(SuperMatrix *mat, int n, int m, int offset,
00141                        std::complex<double> *values, int *rowindex, int* colindex,
00142                        Stype_t stype, Dtype_t dtype, Mtype_t mtype)
00143     {
00144       zCreate_CompCol_Matrix(mat, n, m, offset, reinterpret_cast<doublecomplex*>(values),
00145                              rowindex, colindex, stype, dtype, mtype);
00146     }
00147   };
00148   
00149   template<>
00150   struct SuperMatrixPrinter<std::complex<double> >
00151   {
00152     static void print(char* name, SuperMatrix* mat)
00153     {
00154       zPrint_CompCol_Matrix(name, mat);
00155     }
00156   };
00157 #endif
00158     
00164   template<class M>
00165   class MatrixRowSet
00166   {
00167   public:
00168     // @brief The type of the matrix.
00169     typedef M Matrix;
00170     // @brief The matrix row iterator type.
00171     typedef typename Matrix::ConstRowIterator const_iterator;
00172     
00177     MatrixRowSet(const Matrix& m)
00178       : m_(m)
00179     {}
00180     
00181     // @brief Get the row iterator at the first row.
00182     const_iterator begin() const
00183     {
00184       return m_.begin();
00185     }
00186     //@brief Get the row iterator at the end of all rows.
00187     const_iterator end() const
00188     {
00189       return m_.end();
00190     }
00191   private:
00192     const Matrix& m_;
00193   };
00194   
00202   template<class M, class S>
00203   class MatrixRowSubset
00204   {
00205   public:
00206     /* @brief the type of the matrix class. */
00207     typedef M Matrix;
00208     /* @brief the type of the set of valid row indices. */
00209     typedef S RowIndexSet;
00210     
00216     MatrixRowSubset(const Matrix& m, const RowIndexSet& s)
00217       : m_(m), s_(s)
00218     {}
00219 
00220     const Matrix& matrix() const
00221     {
00222       return m_;
00223     }
00224     
00225     const RowIndexSet& rowIndexSet() const
00226     {
00227       return s_;
00228     }
00229     
00230     // @brief The matrix row iterator type.
00231     class const_iterator
00232       : public ForwardIteratorFacade<const_iterator, const typename Matrix::row_type>
00233     {
00234     public:
00235       const_iterator(typename Matrix::const_iterator firstRow,
00236                      typename RowIndexSet::const_iterator pos)
00237         : firstRow_(firstRow), pos_(pos)
00238       {}
00239       
00240       
00241       const typename Matrix::row_type& dereference() const
00242       {
00243         return *(firstRow_+ *pos_);        
00244       }
00245       bool equals(const const_iterator& o) const
00246       {
00247         return pos_==o.pos_;
00248       }
00249       void increment()
00250       {
00251         ++pos_;
00252       }
00253       typename RowIndexSet::value_type index() const
00254       {
00255         return *pos_;
00256       }
00257       
00258     private:
00259       // @brief Iterator pointing to the first row of the matrix.
00260       typename Matrix::const_iterator firstRow_;
00261       // @brief Iterator pointing to the current row index.
00262       typename RowIndexSet::const_iterator pos_; 
00263     };
00264 
00265     // @brief Get the row iterator at the first row.
00266     const_iterator begin() const
00267     {
00268       return const_iterator(m_.begin(), s_.begin());
00269     }
00270     //@brief Get the row iterator at the end of all rows.
00271     const_iterator end() const
00272     {
00273       return const_iterator(m_.begin(), s_.end());
00274     }
00275     
00276   private:
00277     // @brief The matrix for which we manage the row subset.
00278     const Matrix& m_;
00279     // @brief The set of row indices we manage.
00280     const RowIndexSet& s_;
00281     
00282   };
00283   
00284   template<class T>
00285   struct GetSuperLUType
00286   {
00287   };
00288 
00289   template<>
00290   struct GetSuperLUType<double>
00291   {
00292     static const Dtype_t type;
00293     typedef double float_type;
00294   
00295   };
00296   const Dtype_t GetSuperLUType<double>::type=SLU_D;
00297   
00298   template<>
00299   struct GetSuperLUType<float>
00300   {
00301     static const Dtype_t type;
00302     typedef float float_type;
00303   };
00304 
00305   const Dtype_t GetSuperLUType<float>::type=SLU_S;
00306 
00307   template<>
00308   struct GetSuperLUType<std::complex<double> >
00309   {
00310     static const Dtype_t type;
00311     typedef double float_type;
00312   };
00313 
00314   const Dtype_t GetSuperLUType<std::complex<double> >::type=SLU_Z;
00315 
00316   template<>
00317   struct GetSuperLUType<std::complex<float> >
00318   {
00319     static const Dtype_t type;
00320     typedef float float_type;
00321   
00322   };
00323 
00324   const Dtype_t GetSuperLUType<std::complex<float> >::type=SLU_C;
00329   template<class M>
00330   struct SuperLUMatrix
00331   {
00332   
00333   };
00334 
00335   template<class M>
00336   struct SuperMatrixInitializer
00337   {};
00338 
00339   template<class M, class X, class TM, class TD, class T1>
00340   class SeqOverlappingSchwarz;
00341 
00342   
00343   template<class T>
00344   class SeqOverlappingSchwarzAssembler;
00345 
00346   template<class T>
00347   class SuperLU;
00348 
00352   template<class B, class TA, int n, int m>
00353   class SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
00354   {
00355     template<class M, class X, class TM, class TD, class T1>
00356     friend class SeqOverlappingSchwarz;
00357     friend class SuperMatrixInitializer<BCRSMatrix<FieldMatrix<B,n,m>,TA> >;
00358 
00359   public:
00361     typedef BCRSMatrix<FieldMatrix<B,n,m>,TA> Matrix;
00362 
00363     friend class SeqOverlappingSchwarzAssembler<SuperLU<Matrix> >;
00364     
00365     typedef typename Matrix::size_type size_type;
00366     
00371     SuperLUMatrix(const Matrix& mat);
00372 
00373     SuperLUMatrix();
00374     
00376     ~SuperLUMatrix();
00377     
00379     operator SuperMatrix&()
00380     {
00381       return A;
00382     }
00383     
00385     operator const SuperMatrix&()const
00386     {
00387       return A;
00388     }
00389     bool operator==(const Matrix& mat) const;
00390     
00395     size_type N() const
00396     {
00397       return N_;
00398     }
00399 
00400     size_type nnz() const
00401     {
00402       return Nnz_/n/m;
00403     }
00404     
00409     size_type M() const
00410     {
00411       return M_;
00412     }
00413     
00414     SuperLUMatrix& operator=(const Matrix& mat);
00415     
00416     SuperLUMatrix& operator=(const SuperLUMatrix& mat);
00417 
00424     template<class S>
00425     void setMatrix(const Matrix& mat, const S& mrs);
00427     void free();
00428   private: 
00430     void setMatrix(const Matrix& mat);
00431 
00432     int N_, M_, Nnz_;
00433     B* values;
00434     int* rowindex;
00435     int* colstart;
00436     SuperMatrix A;
00437   };
00438 
00439   template<class T, class A, int n, int m>
00440   void writeCompColMatrixToMatlab(const SuperLUMatrix<BCRSMatrix<FieldMatrix<T,n,m>,A> >& mat,
00441                            std::ostream& os)
00442   {
00443     const SuperMatrix a=static_cast<const SuperMatrix&>(mat);
00444     NCformat *astore = (NCformat *) a.Store;
00445     double* dp = (double*)astore->nzval;
00446 
00447     // remember old flags
00448     std::ios_base::fmtflags oldflags = os.flags();
00449     // set the output format
00450     //os.setf(std::ios_base::scientific, std::ios_base::floatfield);
00451     int oldprec = os.precision();
00452     //os.precision(10);
00453     //dPrint_CompCol_Matrix("A", const_cast<SuperMatrix*>(&a));
00454     
00455     os <<"[";
00456     for(int row=0; row<a.nrow; ++row){
00457       for(int col= 0; col < a.ncol; ++col){
00458         // linear search for col
00459         int i;
00460         for(i=astore->colptr[col]; i < astore->colptr[col+1]; ++i)
00461           if(astore->rowind[i]==row){
00462             os<<dp[i]<<" ";
00463             break;
00464           }
00465         if(i==astore->colptr[col+1])
00466           // entry not present
00467           os<<0<<" ";
00468       }
00469       if(row==a.nrow-1)
00470         os<<"]";
00471       os<<std::endl;
00472     }
00473     // reset the output format
00474     os.flags(oldflags);
00475     os.precision(oldprec);
00476   }
00477   
00478   
00479   template<class T, class A, int n, int m>
00480   class SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >
00481   {
00482     template<class I, class S, class D>
00483     friend class OverlappingSchwarzInitializer;
00484   public:
00485     typedef Dune::BCRSMatrix<FieldMatrix<T,n,m>,A> Matrix;
00486     typedef Dune::SuperLUMatrix<Matrix> SuperLUMatrix;
00487     typedef typename Matrix::row_type::const_iterator CIter;
00488     typedef typename Matrix::size_type size_type;
00489     
00490     SuperMatrixInitializer(SuperLUMatrix& lum);
00491 
00492     SuperMatrixInitializer();
00493     
00494     ~SuperMatrixInitializer();
00495     
00496     template<typename Iter>
00497     void addRowNnz(const Iter& row)const;
00498 
00499     template<typename Iter, typename Set>
00500     void addRowNnz(const Iter& row, const Set& s)const;
00501 
00502     void allocate();
00503         
00504     template<typename Iter>
00505     void countEntries(const Iter& row, const CIter& col) const;
00506     
00507     void countEntries(size_type colidx)const;
00508 
00509     void calcColstart()const;
00510 
00511     template<typename Iter>
00512     void copyValue(const Iter& row, const CIter& col) const;
00513 
00514     void copyValue(const CIter& col, size_type rowindex, size_type colidx)const;
00515 
00516     void createMatrix() const;
00517     
00518   private:
00519     
00520     void allocateMatrixStorage()const;
00521 
00522     void allocateMarker();
00523 
00524     SuperLUMatrix* mat;
00525     size_type cols;
00526     mutable size_type *marker;
00527   };
00528   
00529   template<class T, class A, int n, int m>
00530   SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperMatrixInitializer(SuperLUMatrix& mat_)
00531     : mat(&mat_), cols(mat_.N()), marker(0)
00532   {
00533     mat->Nnz_=0;
00534   }
00535 
00536   template<class T, class A, int n, int m>
00537   SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::SuperMatrixInitializer()
00538     : mat(0), cols(0), marker(0)
00539   {}
00540   
00541   template<class T, class A, int n, int m>
00542   SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::~SuperMatrixInitializer()
00543   {
00544     if(marker)
00545       delete[] marker;
00546   }
00547   
00548   template<class T, class A, int n, int m>
00549   template<typename Iter>
00550   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row)const
00551   {
00552     mat->Nnz_+=row->getsize();
00553   }
00554 
00555   template<class T, class A, int n, int m>
00556   template<typename Iter, typename Map>
00557   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row,
00558                                                                             const Map& indices)const
00559   {
00560     typedef typename  Iter::value_type::const_iterator RIter;
00561     typedef typename Map::const_iterator MIter;
00562     MIter siter =indices.begin();
00563     for(RIter entry=row->begin();entry!=row->end();++entry)
00564       {
00565         for(;siter!=indices.end() && *siter<entry.index(); ++siter);
00566         if(siter==indices.end())
00567           break;
00568         if(*siter==entry.index())
00569           // index is in subdomain
00570           ++mat->Nnz_;
00571       }
00572   }
00573 
00574   template<class T, class A, int n, int m>
00575   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocate()
00576   {
00577     allocateMatrixStorage();
00578     allocateMarker();
00579   }
00580   
00581   template<class T, class A, int n, int m>
00582   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMatrixStorage()const
00583   {
00584     mat->Nnz_*=n*m;
00585     // initialize data
00586     mat->values=new T[mat->Nnz_];
00587     mat->rowindex=new int[mat->Nnz_];
00588     mat->colstart=new int[cols+1];
00589   }
00590 
00591   template<class T, class A, int n, int m>
00592   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMarker()
00593   {          
00594     marker = new typename Matrix::size_type[cols];
00595 
00596     for(size_type i=0; i < cols; ++i)
00597       marker[i]=0;
00598   }
00599   
00600   template<class T, class A, int n, int m>
00601   template<typename Iter>
00602   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::countEntries(const Iter& row, const CIter& col)const
00603   {
00604     countEntries(col.index());
00605     
00606   }
00607   
00608   template<class T, class A, int n, int m>
00609   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::countEntries(size_type colindex)const
00610   {
00611     for(size_type i=0; i < m; ++i){
00612       assert(colindex*m+i<cols);
00613       marker[colindex*m+i]+=n;
00614     }
00615     
00616   }
00617 
00618   template<class T, class A, int n, int m>
00619   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::calcColstart()const
00620   {
00621     mat->colstart[0]=0;
00622     for(size_type i=0; i < cols; ++i){
00623       assert(i<cols);
00624         mat->colstart[i+1]=mat->colstart[i]+marker[i];
00625         marker[i]=mat->colstart[i];
00626     }
00627   }
00628   
00629   template<class T, class A, int n, int m>
00630   template<typename Iter>
00631   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const Iter& row, const CIter& col)const
00632   {
00633     copyValue(col, row.index(), col.index());
00634   }
00635 
00636   template<class T, class A, int n, int m>
00637   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const CIter& col, size_type rowindex, size_type colindex)const
00638   {
00639     for(size_type i=0; i<n;i++){
00640       for(size_type j=0; j<m; j++){
00641         assert(colindex*m+j<cols-1 || (int)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
00642         assert((int)marker[colindex*m+j]<mat->Nnz_);
00643         mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
00644         mat->values[marker[colindex*m+j]]=(*col)[i][j];
00645         ++marker[colindex*m+j]; // index for next entry in column
00646       }
00647     }
00648   }
00649   
00650   template<class T, class A, int n, int m>
00651   void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::createMatrix() const
00652   {
00653     delete[] marker;
00654     marker=0;
00655     SuperMatrixCreateSparseChooser<T>
00656       ::create(&mat->A, mat->N_, mat->M_, mat->colstart[cols], 
00657                mat->values, mat->rowindex, mat->colstart, SLU_NC, 
00658                static_cast<Dtype_t>(GetSuperLUType<T>::type), SLU_GE);
00659   }
00660 
00661   template<class F, class MRS>
00662   void copyToSuperMatrix(F& initializer, const MRS& mrs)
00663   {
00664     typedef typename MRS::const_iterator Iter;
00665     typedef typename  std::iterator_traits<Iter>::value_type::const_iterator CIter;
00666     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
00667       initializer.addRowNnz(row);
00668     
00669     initializer.allocate();
00670     
00671     for(Iter row=mrs.begin(); row!= mrs.end(); ++row){
00672       
00673       for(CIter col=row->begin(); col != row->end(); ++col)
00674         initializer.countEntries(row, col);
00675     }
00676     
00677     initializer.calcColstart();
00678     
00679     for(Iter row=mrs.begin(); row!= mrs.end(); ++row){
00680       for(CIter col=row->begin(); col != row->end(); ++col){
00681           initializer.copyValue(row, col);
00682       }
00683       
00684     }
00685     initializer.createMatrix();
00686   }
00687 
00688   template<class F, class M,class S>
00689   void copyToSuperMatrix(F& initializer, const MatrixRowSubset<M,S>& mrs)
00690   {
00691     typedef MatrixRowSubset<M,S> MRS;
00692     typedef typename MRS::RowIndexSet SIS;
00693     typedef typename SIS::const_iterator SIter;
00694     typedef typename MRS::const_iterator Iter;
00695     typedef typename std::iterator_traits<Iter>::value_type row_type;
00696     typedef typename row_type::const_iterator CIter;
00697 
00698     // Calculate upper Bound for nonzeros
00699     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
00700       initializer.addRowNnz(row, mrs.rowIndexSet());
00701 
00702     initializer.allocate();
00703 
00704     typedef typename MRS::Matrix::size_type size_type;
00705 
00706     // A vector containing the corresponding indices in
00707     // the to create submatrix.
00708     // If an entry is the maximum of size_type then this index will not appear in 
00709     // the submatrix.
00710     std::vector<size_type> subMatrixIndex(mrs.matrix().N(), 
00711                                          std::numeric_limits<size_type>::max());
00712     size_type s=0;
00713     for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
00714       subMatrixIndex[*index]=s++;
00715     
00716     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
00717       for(CIter col=row->begin(); col != row->end(); ++col){
00718         if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
00719           // This column is in our subset (use submatrix column index)
00720           initializer.countEntries(subMatrixIndex[col.index()]);
00721     }
00722     
00723     initializer.calcColstart();
00724     
00725     for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
00726       for(CIter col=row->begin(); col != row->end(); ++col){
00727         if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
00728           // This value is in our submatrix -> copy (use submatrix indices
00729           initializer.copyValue(col, subMatrixIndex[row.index()], subMatrixIndex[col.index()]);
00730     }
00731     initializer.createMatrix();
00732   }
00733 
00734 #ifndef DOXYGEN
00735   
00736   template<class B, class TA, int n, int m>
00737   bool  SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator==(const BCRSMatrix<FieldMatrix<B,n,m>,TA>& mat) const
00738   {
00739     const NCformat* S=static_cast<const NCformat *>(A.Store);
00740     for(size_type col=0; col < M(); ++col){
00741       for(int j=S->colptr[col]; j < S->colptr[col+1]; ++j){
00742         int row=S->rowind[j];
00743         if((mat[row/n][col/m])[row%n][col%m]!=reinterpret_cast<B*>(S->nzval)[j]){
00744           std::cerr<<" bcrs["<<row/n<<"]["<<col/m<<"]["<<row%n<<"]["<<row%m
00745                    <<"]="<<(mat[row/n][col/m])[row%n][col%m]<<" super["<<row<<"]["<<col<<"]="<<reinterpret_cast<B*>(S->nzval)[j];
00746           
00747           return false;
00748         }
00749       }
00750     }
00751     
00752     return true;
00753   }
00754   
00755 #endif // DOYXGEN
00756 
00757   template<class B, class TA, int n, int m>
00758   bool operator==(SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >& sla, BCRSMatrix<FieldMatrix<B,n,m>,TA>& a)
00759   {
00760     return a==sla;
00761   }
00762 
00763 #ifndef DOXYGEN
00764   
00765   template<class B, class TA, int n, int m>
00766   SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::SuperLUMatrix()
00767     : N_(0), M_(0), Nnz_(0), values(0), rowindex(0), colstart(0)
00768   {}
00769 
00770   template<class B, class TA, int n, int m>
00771   SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
00772   ::SuperLUMatrix(const Matrix& mat)
00773     : N_(n*mat.N()), M_(m*mat.M()), Nnz_(n*m*mat.Nnz())
00774   {
00775     
00776   }
00777 
00778   template<class B, class TA, int n, int m>
00779   SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
00780   SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const Matrix& mat)
00781   {
00782     if(N_+M_+Nnz_!=0)
00783       free();
00784     setMatrix(mat);
00785     return *this;
00786   }
00787 
00788   template<class B, class TA, int n, int m>
00789   SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
00790   SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const SuperLUMatrix& mat)
00791   {
00792     if(N_+M_+Nnz_!=0)
00793       free();
00794     N_=mat.N_;
00795     M_=mat.M_;
00796     Nnz_= mat.Nnz_;
00797     if(M_>0){
00798       colstart=new int[M_+1];
00799       for(int i=0; i<=M_; ++i)
00800         colstart[i]=mat.colstart[i];
00801     }
00802     
00803     if(Nnz_>0){
00804       values = new B[Nnz_];
00805       rowindex = new int[Nnz_];
00806       
00807       for(int i=0; i<Nnz_; ++i)
00808         values[i]=mat.values[i];
00809       
00810       for(int i=0; i<Nnz_; ++i)
00811         rowindex[i]=mat.rowindex[i];
00812     }
00813     if(M_+Nnz_>0)
00814       dCreate_CompCol_Matrix(&A, N_, M_, Nnz_, 
00815                            values, rowindex, colstart, SLU_NC, static_cast<Dtype_t>(GetSuperLUType<B>::type), SLU_GE);
00816     return *this;
00817   }
00818 
00819   template<class B, class TA, int n, int m>
00820   void SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
00821   ::setMatrix(const Matrix& mat)
00822   {
00823     N_=n*mat.N();
00824     M_=m*mat.M();
00825     SuperMatrixInitializer<Matrix> initializer(*this);
00826     
00827     copyToSuperMatrix(initializer, MatrixRowSet<Matrix>(mat));
00828     
00829 #ifdef DUNE_ISTL_WITH_CHECKING
00830     char name[] = {'A',0};
00831     if(N_<0)
00832       SuperMatrixPrinter<B>::print(name,&A);
00833     assert(*this==mat);
00834 #endif
00835   }
00836 
00837    template<class B, class TA, int n, int m>
00838    template<class S>
00839   void SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
00840    ::setMatrix(const Matrix& mat, const S& mrs)
00841   {
00842     if(N_+M_+Nnz_!=0)
00843       free();
00844     N_=mrs.size()*n;
00845     M_=mrs.size()*m;
00846     SuperMatrixInitializer<Matrix> initializer(*this);
00847     
00848     copyToSuperMatrix(initializer, MatrixRowSubset<Matrix,S>(mat,mrs));
00849     
00850 #ifdef DUNE_ISTL_WITH_CHECKING
00851     char name[] = {'A',0};
00852     if(N_<0)
00853       SuperMatrixPrinter<B>::print(name,&A);
00854 #endif
00855   }
00856 
00857   template<class B, class TA, int n, int m>
00858   SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::~SuperLUMatrix()
00859   {
00860     if(N_+M_+Nnz_!=0)
00861       free();
00862   }
00863 
00864   template<class B, class TA, int n, int m>
00865   void SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::free()
00866   {
00867     delete[] values;
00868     delete[] rowindex;
00869     delete[] colstart;
00870     SUPERLU_FREE(A.Store);
00871     N_=M_=Nnz_=0;
00872   }
00873 
00874 #endif // DOXYGEN
00875 
00876 }
00877 #endif  
00878 #endif