dune-istl  2.2.1
supermatrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=8 sw=2 sts=2:
3 #ifndef DUNE_SUPERLUMATRIX_HH
4 #define DUNE_SUPERLUMATRIX_HH
5 
6 #if HAVE_SUPERLU
7 #ifdef SUPERLU_POST_2005_VERSION
8 
9 #ifndef SUPERLU_NTYPE
10 #define SUPERLU_NTYPE 1
11 #endif
12 
13 #if SUPERLU_NTYPE==0
14 #include "slu_sdefs.h"
15 #endif
16 
17 #if SUPERLU_NTYPE==1
18 #include "slu_ddefs.h"
19 #endif
20 
21 #if SUPERLU_NTYPE==2
22 #include "slu_cdefs.h"
23 #endif
24 
25 #if SUPERLU_NTYPE>=3
26 #include "slu_zdefs.h"
27 #endif
28 
29 #else
30 
31 #if SUPERLU_NTYPE==0
32 #include "ssp_defs.h"
33 #endif
34 
35 #if SUPERLU_NTYPE==1
36 #include "dsp_defs.h"
37 #endif
38 
39 #if SUPERLU_NTYPE==2
40 #include "csp_defs.h"
41 #endif
42 
43 #if SUPERLU_NTYPE>=3
44 #include "zsp_defs.h"
45 #endif
46 
47 #endif
48 #include"bcrsmatrix.hh"
49 #include"bvector.hh"
50 #include<dune/common/fmatrix.hh>
51 #include<dune/common/fvector.hh>
52 #include<dune/common/typetraits.hh>
53 #include<limits>
54 
55 namespace Dune
56 {
57 
58  template<class T>
60  {
61  };
62 
63  template<class T>
65  {
66  };
67 
68 #if SUPERLU_NTYPE==0
69  template<>
71  {
72  static void create(SuperMatrix *mat, int n, int m, int offset,
73  float *values, int *rowindex, int* colindex,
74  Stype_t stype, Dtype_t dtype, Mtype_t mtype)
75  {
76  sCreate_CompCol_Matrix(mat, n, m, offset, values, rowindex, colindex,
77  stype, dtype, mtype);
78  }
79  };
80 
81  template<>
82  struct SuperMatrixPrinter<float>
83  {
84  static void print(char* name, SuperMatrix* mat)
85  {
86  sPrint_CompCol_Matrix(name, mat);
87  }
88  };
89 #endif
90 
91 #if SUPERLU_NTYPE==1
92  template<>
93  struct SuperMatrixCreateSparseChooser<double>
94  {
95  static void create(SuperMatrix *mat, int n, int m, int offset,
96  double *values, int *rowindex, int* colindex,
97  Stype_t stype, Dtype_t dtype, Mtype_t mtype)
98  {
99  dCreate_CompCol_Matrix(mat, n, m, offset, values, rowindex, colindex,
100  stype, dtype, mtype);
101  }
102  };
103 
104  template<>
105  struct SuperMatrixPrinter<double>
106  {
107  static void print(char* name, SuperMatrix* mat)
108  {
109  dPrint_CompCol_Matrix(name, mat);
110  }
111  };
112 #endif
113 
114 #if SUPERLU_NTYPE==2
115  template<>
116  struct SuperMatrixCreateSparseChooser<std::complex<float> >
117  {
118  static void create(SuperMatrix *mat, int n, int m, int offset,
119  std::complex<float> *values, int *rowindex, int* colindex,
120  Stype_t stype, Dtype_t dtype, Mtype_t mtype)
121  {
122  cCreate_CompCol_Matrix(mat, n, m, offset, reinterpret_cast< ::complex*>(values),
123  rowindex, colindex, stype, dtype, mtype);
124  }
125  };
126 
127  template<>
128  struct SuperMatrixPrinter<std::complex<float> >
129  {
130  static void print(char* name, SuperMatrix* mat)
131  {
132  cPrint_CompCol_Matrix(name, mat);
133  }
134  };
135 #endif
136 
137 #if SUPERLU_NTYPE>=3
138  template<>
139  struct SuperMatrixCreateSparseChooser<std::complex<double> >
140  {
141  static void create(SuperMatrix *mat, int n, int m, int offset,
142  std::complex<double> *values, int *rowindex, int* colindex,
143  Stype_t stype, Dtype_t dtype, Mtype_t mtype)
144  {
145  zCreate_CompCol_Matrix(mat, n, m, offset, reinterpret_cast<doublecomplex*>(values),
146  rowindex, colindex, stype, dtype, mtype);
147  }
148  };
149 
150  template<>
151  struct SuperMatrixPrinter<std::complex<double> >
152  {
153  static void print(char* name, SuperMatrix* mat)
154  {
155  zPrint_CompCol_Matrix(name, mat);
156  }
157  };
158 #endif
159 
165  template<class M>
167  {
168  public:
169  // @brief The type of the matrix.
170  typedef M Matrix;
171  // @brief The matrix row iterator type.
173 
179  : m_(m)
180  {}
181 
182  // @brief Get the row iterator at the first row.
184  {
185  return m_.begin();
186  }
187  //@brief Get the row iterator at the end of all rows.
189  {
190  return m_.end();
191  }
192  private:
193  const Matrix& m_;
194  };
195 
203  template<class M, class S>
205  {
206  public:
207  /* @brief the type of the matrix class. */
208  typedef M Matrix;
209  /* @brief the type of the set of valid row indices. */
210  typedef S RowIndexSet;
211 
217  MatrixRowSubset(const Matrix& m, const RowIndexSet& s)
218  : m_(m), s_(s)
219  {}
220 
221  const Matrix& matrix() const
222  {
223  return m_;
224  }
225 
226  const RowIndexSet& rowIndexSet() const
227  {
228  return s_;
229  }
230 
231  // @brief The matrix row iterator type.
233  : public ForwardIteratorFacade<const_iterator, const typename Matrix::row_type>
234  {
235  public:
236  const_iterator(typename Matrix::const_iterator firstRow,
237  typename RowIndexSet::const_iterator pos)
238  : firstRow_(firstRow), pos_(pos)
239  {}
240 
241 
242  const typename Matrix::row_type& dereference() const
243  {
244  return *(firstRow_+ *pos_);
245  }
246  bool equals(const const_iterator& o) const
247  {
248  return pos_==o.pos_;
249  }
250  void increment()
251  {
252  ++pos_;
253  }
254  typename RowIndexSet::value_type index() const
255  {
256  return *pos_;
257  }
258 
259  private:
260  // @brief Iterator pointing to the first row of the matrix.
261  typename Matrix::const_iterator firstRow_;
262  // @brief Iterator pointing to the current row index.
263  typename RowIndexSet::const_iterator pos_;
264  };
265 
266  // @brief Get the row iterator at the first row.
268  {
269  return const_iterator(m_.begin(), s_.begin());
270  }
271  //@brief Get the row iterator at the end of all rows.
273  {
274  return const_iterator(m_.begin(), s_.end());
275  }
276 
277  private:
278  // @brief The matrix for which we manage the row subset.
279  const Matrix& m_;
280  // @brief The set of row indices we manage.
281  const RowIndexSet& s_;
282 
283  };
284 
285  template<class T>
287  {
288  static const Dtype_t type;
289  };
290 
291  template<class T>
293  {
294  };
295 
296  template<class T>
297  const Dtype_t BaseGetSuperLUType<T>::type =
298  Dune::is_same<T,float>::value ? SLU_S :
299  ( Dune::is_same<T,std::complex<double> >::value ? SLU_Z :
300  ( Dune::is_same<T,std::complex<float> >::value ? SLU_C : SLU_D ));
301 
302  template<>
303  struct GetSuperLUType<double>
304  : public BaseGetSuperLUType<double>
305  {
306  typedef double float_type;
307  };
308 
309  template<>
310  struct GetSuperLUType<float>
311  : public BaseGetSuperLUType<float>
312  {
313  typedef float float_type;
314  };
315 
316  template<>
317  struct GetSuperLUType<std::complex<double> >
318  : public BaseGetSuperLUType<std::complex<double> >
319  {
320  typedef double float_type;
321  };
322 
323  template<>
324  struct GetSuperLUType<std::complex<float> >
325  : public BaseGetSuperLUType<std::complex<float> >
326  {
327  typedef float float_type;
328 
329  };
330 
335  template<class M>
337  {
338 
339  };
340 
341  template<class M>
343  {};
344 
345  template<class M, class X, class TM, class TD, class T1>
346  class SeqOverlappingSchwarz;
347 
348 
349  template<class T>
351 
352  template<class T>
353  class SuperLU;
354 
358  template<class B, class TA, int n, int m>
359  class SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
360  {
361  template<class M, class X, class TM, class TD, class T1>
362  friend class SeqOverlappingSchwarz;
363  friend struct SuperMatrixInitializer<BCRSMatrix<FieldMatrix<B,n,m>,TA> >;
364 
365  public:
368 
370 
371  typedef typename Matrix::size_type size_type;
372 
377  SuperLUMatrix(const Matrix& mat);
378 
379  SuperLUMatrix();
380 
382  ~SuperLUMatrix();
383 
385  operator SuperMatrix&()
386  {
387  return A;
388  }
389 
391  operator const SuperMatrix&()const
392  {
393  return A;
394  }
395  bool operator==(const Matrix& mat) const;
396 
401  size_type N() const
402  {
403  return N_;
404  }
405 
406  size_type nnz() const
407  {
408  return Nnz_/n/m;
409  }
410 
415  size_type M() const
416  {
417  return M_;
418  }
419 
420  SuperLUMatrix& operator=(const Matrix& mat);
421 
422  SuperLUMatrix& operator=(const SuperLUMatrix& mat);
423 
430  template<class S>
431  void setMatrix(const Matrix& mat, const S& mrs);
433  void free();
434  private:
436  void setMatrix(const Matrix& mat);
437 
438  int N_, M_, Nnz_;
439  B* values;
440  int* rowindex;
441  int* colstart;
442  SuperMatrix A;
443  };
444 
445  template<class T, class A, int n, int m>
446  void writeCompColMatrixToMatlab(const SuperLUMatrix<BCRSMatrix<FieldMatrix<T,n,m>,A> >& mat,
447  std::ostream& os)
448  {
449  const SuperMatrix a=static_cast<const SuperMatrix&>(mat);
450  NCformat *astore = (NCformat *) a.Store;
451  double* dp = (double*)astore->nzval;
452 
453  // remember old flags
454  std::ios_base::fmtflags oldflags = os.flags();
455  // set the output format
456  //os.setf(std::ios_base::scientific, std::ios_base::floatfield);
457  int oldprec = os.precision();
458  //os.precision(10);
459  //dPrint_CompCol_Matrix("A", const_cast<SuperMatrix*>(&a));
460 
461  os <<"[";
462  for(int row=0; row<a.nrow; ++row){
463  for(int col= 0; col < a.ncol; ++col){
464  // linear search for col
465  int i;
466  for(i=astore->colptr[col]; i < astore->colptr[col+1]; ++i)
467  if(astore->rowind[i]==row){
468  os<<dp[i]<<" ";
469  break;
470  }
471  if(i==astore->colptr[col+1])
472  // entry not present
473  os<<0<<" ";
474  }
475  if(row==a.nrow-1)
476  os<<"]";
477  os<<std::endl;
478  }
479  // reset the output format
480  os.flags(oldflags);
481  os.precision(oldprec);
482  }
483 
484 
485  template<class T, class A, int n, int m>
486  class SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >
487  {
488  template<class I, class S, class D>
490  public:
494  typedef typename Matrix::size_type size_type;
495 
497 
499 
501 
502  template<typename Iter>
503  void addRowNnz(const Iter& row)const;
504 
505  template<typename Iter, typename Set>
506  void addRowNnz(const Iter& row, const Set& s)const;
507 
508  void allocate();
509 
510  template<typename Iter>
511  void countEntries(const Iter& row, const CIter& col) const;
512 
513  void countEntries(size_type colidx)const;
514 
515  void calcColstart()const;
516 
517  template<typename Iter>
518  void copyValue(const Iter& row, const CIter& col) const;
519 
520  void copyValue(const CIter& col, size_type rowindex, size_type colidx)const;
521 
522  void createMatrix() const;
523 
524  private:
525 
526  void allocateMatrixStorage()const;
527 
528  void allocateMarker();
529 
531  size_type cols;
532  mutable size_type *marker;
533  };
534 
535  template<class T, class A, int n, int m>
537  : mat(&mat_), cols(mat_.N()), marker(0)
538  {
539  mat->Nnz_=0;
540  }
541 
542  template<class T, class A, int n, int m>
544  : mat(0), cols(0), marker(0)
545  {}
546 
547  template<class T, class A, int n, int m>
549  {
550  if(marker)
551  delete[] marker;
552  }
553 
554  template<class T, class A, int n, int m>
555  template<typename Iter>
556  void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::addRowNnz(const Iter& row)const
557  {
558  mat->Nnz_+=row->getsize();
559  }
560 
561  template<class T, class A, int n, int m>
562  template<typename Iter, typename Map>
564  const Map& indices)const
565  {
566  typedef typename Iter::value_type::const_iterator RIter;
567  typedef typename Map::const_iterator MIter;
568  MIter siter =indices.begin();
569  for(RIter entry=row->begin();entry!=row->end();++entry)
570  {
571  for(;siter!=indices.end() && *siter<entry.index(); ++siter);
572  if(siter==indices.end())
573  break;
574  if(*siter==entry.index())
575  // index is in subdomain
576  ++mat->Nnz_;
577  }
578  }
579 
580  template<class T, class A, int n, int m>
582  {
583  allocateMatrixStorage();
584  allocateMarker();
585  }
586 
587  template<class T, class A, int n, int m>
588  void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMatrixStorage()const
589  {
590  mat->Nnz_*=n*m;
591  // initialize data
592  mat->values=new T[mat->Nnz_];
593  mat->rowindex=new int[mat->Nnz_];
594  mat->colstart=new int[cols+1];
595  }
596 
597  template<class T, class A, int n, int m>
598  void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::allocateMarker()
599  {
600  marker = new typename Matrix::size_type[cols];
601 
602  for(size_type i=0; i < cols; ++i)
603  marker[i]=0;
604  }
605 
606  template<class T, class A, int n, int m>
607  template<typename Iter>
609  {
610  countEntries(col.index());
611 
612  }
613 
614  template<class T, class A, int n, int m>
616  {
617  for(size_type i=0; i < m; ++i){
618  assert(colindex*m+i<cols);
619  marker[colindex*m+i]+=n;
620  }
621 
622  }
623 
624  template<class T, class A, int n, int m>
626  {
627  mat->colstart[0]=0;
628  for(size_type i=0; i < cols; ++i){
629  assert(i<cols);
630  mat->colstart[i+1]=mat->colstart[i]+marker[i];
631  marker[i]=mat->colstart[i];
632  }
633  }
634 
635  template<class T, class A, int n, int m>
636  template<typename Iter>
637  void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const Iter& row, const CIter& col)const
638  {
639  copyValue(col, row.index(), col.index());
640  }
641 
642  template<class T, class A, int n, int m>
643  void SuperMatrixInitializer<BCRSMatrix<FieldMatrix<T,n,m>,A> >::copyValue(const CIter& col, size_type rowindex, size_type colindex)const
644  {
645  for(size_type i=0; i<n;i++){
646  for(size_type j=0; j<m; j++){
647  assert(colindex*m+j<cols-1 || (int)marker[colindex*m+j]<mat->colstart[colindex*m+j+1]);
648  assert((int)marker[colindex*m+j]<mat->Nnz_);
649  mat->rowindex[marker[colindex*m+j]]=rowindex*n+i;
650  mat->values[marker[colindex*m+j]]=(*col)[i][j];
651  ++marker[colindex*m+j]; // index for next entry in column
652  }
653  }
654  }
655 
656  template<class T, class A, int n, int m>
658  {
659  delete[] marker;
660  marker=0;
662  ::create(&mat->A, mat->N_, mat->M_, mat->colstart[cols],
663  mat->values, mat->rowindex, mat->colstart, SLU_NC,
664  static_cast<Dtype_t>(GetSuperLUType<T>::type), SLU_GE);
665  }
666 
667  template<class F, class MRS>
668  void copyToSuperMatrix(F& initializer, const MRS& mrs)
669  {
670  typedef typename MRS::const_iterator Iter;
671  typedef typename std::iterator_traits<Iter>::value_type::const_iterator CIter;
672  for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
673  initializer.addRowNnz(row);
674 
675  initializer.allocate();
676 
677  for(Iter row=mrs.begin(); row!= mrs.end(); ++row){
678 
679  for(CIter col=row->begin(); col != row->end(); ++col)
680  initializer.countEntries(row, col);
681  }
682 
683  initializer.calcColstart();
684 
685  for(Iter row=mrs.begin(); row!= mrs.end(); ++row){
686  for(CIter col=row->begin(); col != row->end(); ++col){
687  initializer.copyValue(row, col);
688  }
689 
690  }
691  initializer.createMatrix();
692  }
693 
694  template<class F, class M,class S>
695  void copyToSuperMatrix(F& initializer, const MatrixRowSubset<M,S>& mrs)
696  {
697  typedef MatrixRowSubset<M,S> MRS;
698  typedef typename MRS::RowIndexSet SIS;
699  typedef typename SIS::const_iterator SIter;
700  typedef typename MRS::const_iterator Iter;
701  typedef typename std::iterator_traits<Iter>::value_type row_type;
702  typedef typename row_type::const_iterator CIter;
703 
704  // Calculate upper Bound for nonzeros
705  for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
706  initializer.addRowNnz(row, mrs.rowIndexSet());
707 
708  initializer.allocate();
709 
710  typedef typename MRS::Matrix::size_type size_type;
711 
712  // A vector containing the corresponding indices in
713  // the to create submatrix.
714  // If an entry is the maximum of size_type then this index will not appear in
715  // the submatrix.
716  std::vector<size_type> subMatrixIndex(mrs.matrix().N(),
717  std::numeric_limits<size_type>::max());
718  size_type s=0;
719  for(SIter index = mrs.rowIndexSet().begin(); index!=mrs.rowIndexSet().end(); ++index)
720  subMatrixIndex[*index]=s++;
721 
722  for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
723  for(CIter col=row->begin(); col != row->end(); ++col){
724  if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
725  // This column is in our subset (use submatrix column index)
726  initializer.countEntries(subMatrixIndex[col.index()]);
727  }
728 
729  initializer.calcColstart();
730 
731  for(Iter row=mrs.begin(); row!= mrs.end(); ++row)
732  for(CIter col=row->begin(); col != row->end(); ++col){
733  if(subMatrixIndex[col.index()]!=std::numeric_limits<size_type>::max())
734  // This value is in our submatrix -> copy (use submatrix indices
735  initializer.copyValue(col, subMatrixIndex[row.index()], subMatrixIndex[col.index()]);
736  }
737  initializer.createMatrix();
738  }
739 
740 #ifndef DOXYGEN
741 
742  template<class B, class TA, int n, int m>
743  bool SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator==(const BCRSMatrix<FieldMatrix<B,n,m>,TA>& mat) const
744  {
745  const NCformat* S=static_cast<const NCformat *>(A.Store);
746  for(size_type col=0; col < M(); ++col){
747  for(int j=S->colptr[col]; j < S->colptr[col+1]; ++j){
748  int row=S->rowind[j];
749  if((mat[row/n][col/m])[row%n][col%m]!=reinterpret_cast<B*>(S->nzval)[j]){
750  std::cerr<<" bcrs["<<row/n<<"]["<<col/m<<"]["<<row%n<<"]["<<row%m
751  <<"]="<<(mat[row/n][col/m])[row%n][col%m]<<" super["<<row<<"]["<<col<<"]="<<reinterpret_cast<B*>(S->nzval)[j];
752 
753  return false;
754  }
755  }
756  }
757 
758  return true;
759  }
760 
761 #endif // DOYXGEN
762 
763  template<class B, class TA, int n, int m>
764  bool operator==(SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >& sla, BCRSMatrix<FieldMatrix<B,n,m>,TA>& a)
765  {
766  return a==sla;
767  }
768 
769 #ifndef DOXYGEN
770 
771  template<class B, class TA, int n, int m>
772  SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::SuperLUMatrix()
773  : N_(0), M_(0), Nnz_(0), values(0), rowindex(0), colstart(0)
774  {}
775 
776  template<class B, class TA, int n, int m>
777  SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
778  ::SuperLUMatrix(const Matrix& mat)
779  : N_(n*mat.N()), M_(m*mat.M()), Nnz_(n*m*mat.Nnz())
780  {
781 
782  }
783 
784  template<class B, class TA, int n, int m>
785  SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
786  SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const Matrix& mat)
787  {
788  if(N_+M_+Nnz_!=0)
789  free();
790  setMatrix(mat);
791  return *this;
792  }
793 
794  template<class B, class TA, int n, int m>
795  SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >&
796  SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::operator=(const SuperLUMatrix& mat)
797  {
798  if(N_+M_+Nnz_!=0)
799  free();
800  N_=mat.N_;
801  M_=mat.M_;
802  Nnz_= mat.Nnz_;
803  if(M_>0){
804  colstart=new int[M_+1];
805  for(int i=0; i<=M_; ++i)
806  colstart[i]=mat.colstart[i];
807  }
808 
809  if(Nnz_>0){
810  values = new B[Nnz_];
811  rowindex = new int[Nnz_];
812 
813  for(int i=0; i<Nnz_; ++i)
814  values[i]=mat.values[i];
815 
816  for(int i=0; i<Nnz_; ++i)
817  rowindex[i]=mat.rowindex[i];
818  }
819  if(M_+Nnz_>0)
820  dCreate_CompCol_Matrix(&A, N_, M_, Nnz_,
821  values, rowindex, colstart, SLU_NC, static_cast<Dtype_t>(GetSuperLUType<B>::type), SLU_GE);
822  return *this;
823  }
824 
825  template<class B, class TA, int n, int m>
826  void SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
827  ::setMatrix(const Matrix& mat)
828  {
829  N_=n*mat.N();
830  M_=m*mat.M();
831  SuperMatrixInitializer<Matrix> initializer(*this);
832 
833  copyToSuperMatrix(initializer, MatrixRowSet<Matrix>(mat));
834 
835 #ifdef DUNE_ISTL_WITH_CHECKING
836  char name[] = {'A',0};
837  if(N_<0)
838  SuperMatrixPrinter<B>::print(name,&A);
839  assert(*this==mat);
840 #endif
841  }
842 
843  template<class B, class TA, int n, int m>
844  template<class S>
845  void SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >
846  ::setMatrix(const Matrix& mat, const S& mrs)
847  {
848  if(N_+M_+Nnz_!=0)
849  free();
850  N_=mrs.size()*n;
851  M_=mrs.size()*m;
852  SuperMatrixInitializer<Matrix> initializer(*this);
853 
854  copyToSuperMatrix(initializer, MatrixRowSubset<Matrix,S>(mat,mrs));
855 
856 #ifdef DUNE_ISTL_WITH_CHECKING
857  char name[] = {'A',0};
858  if(N_<0)
859  SuperMatrixPrinter<B>::print(name,&A);
860 #endif
861  }
862 
863  template<class B, class TA, int n, int m>
864  SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::~SuperLUMatrix()
865  {
866  if(N_+M_+Nnz_!=0)
867  free();
868  }
869 
870  template<class B, class TA, int n, int m>
871  void SuperLUMatrix<BCRSMatrix<FieldMatrix<B,n,m>,TA> >::free()
872  {
873  delete[] values;
874  delete[] rowindex;
875  delete[] colstart;
876  SUPERLU_FREE(A.Store);
877  N_=M_=Nnz_=0;
878  }
879 
880 #endif // DOXYGEN
881 
882 }
883 #endif
884 #endif