dune-istl
2.2.0
|
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