dune-istl
2.2.0
|
00001 #ifndef DUNE_BLOCK_TRIDIAGONAL_MATRIX_HH 00002 #define DUNE_BLOCK_TRIDIAGONAL_MATRIX_HH 00003 00004 #include <dune/common/fmatrix.hh> 00005 #include <dune/istl/bcrsmatrix.hh> 00006 00012 namespace Dune { 00022 template <class B, class A=std::allocator<B> > 00023 class BTDMatrix : public BCRSMatrix<B,A> 00024 { 00025 public: 00026 00027 //===== type definitions and constants 00028 00030 typedef typename B::field_type field_type; 00031 00033 typedef B block_type; 00034 00036 typedef A allocator_type; 00037 00039 //typedef BCRSMatrix<B,A>::row_type row_type; 00040 00042 typedef typename A::size_type size_type; 00043 00045 enum {blocklevel = B::blocklevel+1}; 00046 00048 BTDMatrix() : BCRSMatrix<B,A>() {} 00049 00050 explicit BTDMatrix(int size) 00051 : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random) 00052 { 00053 // special handling for 1x1 matrices 00054 if (size==1) { 00055 00056 this->BCRSMatrix<B,A>::setrowsize(0, 1); 00057 this->BCRSMatrix<B,A>::endrowsizes(); 00058 00059 this->BCRSMatrix<B,A>::addindex(0, 0); 00060 this->BCRSMatrix<B,A>::endindices(); 00061 00062 return; 00063 } 00064 00065 // Set number of entries for each row 00066 this->BCRSMatrix<B,A>::setrowsize(0, 2); 00067 00068 for (int i=1; i<size-1; i++) 00069 this->BCRSMatrix<B,A>::setrowsize(i, 3); 00070 00071 this->BCRSMatrix<B,A>::setrowsize(size-1, 2); 00072 00073 this->BCRSMatrix<B,A>::endrowsizes(); 00074 00075 // The actual entries for each row 00076 this->BCRSMatrix<B,A>::addindex(0, 0); 00077 this->BCRSMatrix<B,A>::addindex(0, 1); 00078 00079 for (int i=1; i<size-1; i++) { 00080 this->BCRSMatrix<B,A>::addindex(i, i-1); 00081 this->BCRSMatrix<B,A>::addindex(i, i ); 00082 this->BCRSMatrix<B,A>::addindex(i, i+1); 00083 } 00084 00085 this->BCRSMatrix<B,A>::addindex(size-1, size-2); 00086 this->BCRSMatrix<B,A>::addindex(size-1, size-1); 00087 00088 this->BCRSMatrix<B,A>::endindices(); 00089 00090 } 00091 00093 BTDMatrix& operator= (const BTDMatrix& other) { 00094 this->BCRSMatrix<B,A>::operator=(other); 00095 return *this; 00096 } 00097 00099 BTDMatrix& operator= (const field_type& k) { 00100 this->BCRSMatrix<B,A>::operator=(k); 00101 return *this; 00102 } 00103 00109 template <class V> 00110 void solve (V& x, const V& rhs) const { 00111 00112 // special handling for 1x1 matrices. The generic algorithm doesn't work for them 00113 if (this->N()==1) { 00114 (*this)[0][0].solve(x[0],rhs[0]); 00115 return; 00116 } 00117 00118 // Make copies of the rhs and the right matrix band 00119 V d = rhs; 00120 std::vector<block_type> c(this->N()-1); 00121 for (size_t i=0; i<this->N()-1; i++) 00122 c[i] = (*this)[i][i+1]; 00123 00124 /* Modify the coefficients. */ 00125 block_type a_00_inv; 00126 FMatrixHelp::invertMatrix((*this)[0][0], a_00_inv); 00127 00128 //c[0] /= (*this)[0][0]; /* Division by zero risk. */ 00129 block_type c_0_tmp = c[0]; 00130 FMatrixHelp::multMatrix(a_00_inv, c_0_tmp, c[0]); 00131 00132 // d = a^{-1} d /* Division by zero would imply a singular matrix. */ 00133 typename V::block_type d_0_tmp = d[0]; 00134 (*this)[0][0].solve(d[0], d_0_tmp); 00135 00136 for (unsigned int i = 1; i < this->N(); i++) { 00137 00138 // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1} 00139 block_type tmp; 00140 FMatrixHelp::multMatrix(c[i-1], (*this)[i][i-1], tmp); 00141 block_type id = (*this)[i][i]; 00142 id -= tmp; 00143 id.invert(); /* Division by zero risk. */ 00144 00145 if (i<c.size()) { 00146 // c[i] *= id 00147 tmp = c[i]; 00148 FMatrixHelp::multMatrix(tmp, id, c[i]); /* Last value calculated is redundant. */ 00149 } 00150 00151 // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id; 00152 (*this)[i][i-1].mmv(d[i-1], d[i]); 00153 typename V::block_type tmpVec = d[i]; 00154 id.mv(tmpVec, d[i]); 00155 //d[i] *= id; 00156 00157 } 00158 00159 /* Now back substitute. */ 00160 x[this->N() - 1] = d[this->N() - 1]; 00161 for (int i = this->N() - 2; i >= 0; i--) { 00162 //x[i] = d[i] - c[i] * x[i + 1]; 00163 x[i] = d[i]; 00164 c[i].mmv(x[i+1], x[i]); 00165 } 00166 00167 } 00168 00169 private: 00170 00171 // //////////////////////////////////////////////////////////////////////////// 00172 // The following methods from the base class should now actually be called 00173 // //////////////////////////////////////////////////////////////////////////// 00174 00175 // createbegin and createend should be in there, too, but I can't get it to compile 00176 // BCRSMatrix<B,A>::CreateIterator createbegin () {} 00177 // BCRSMatrix<B,A>::CreateIterator createend () {} 00178 void setrowsize (size_type i, size_type s) {} 00179 void incrementrowsize (size_type i) {} 00180 void endrowsizes () {} 00181 void addindex (size_type row, size_type col) {} 00182 void endindices () {} 00183 }; 00186 } // end namespace Dune 00187 00188 #endif