dune-istl  2.2.0
btdmatrix.hh
Go to the documentation of this file.
1 #ifndef DUNE_BLOCK_TRIDIAGONAL_MATRIX_HH
2 #define DUNE_BLOCK_TRIDIAGONAL_MATRIX_HH
3 
4 #include <dune/common/fmatrix.hh>
6 
12 namespace Dune {
22  template <class B, class A=std::allocator<B> >
23 class BTDMatrix : public BCRSMatrix<B,A>
24 {
25 public:
26 
27  //===== type definitions and constants
28 
30  typedef typename B::field_type field_type;
31 
33  typedef B block_type;
34 
36  typedef A allocator_type;
37 
39  //typedef BCRSMatrix<B,A>::row_type row_type;
40 
42  typedef typename A::size_type size_type;
43 
45  enum {blocklevel = B::blocklevel+1};
46 
48  BTDMatrix() : BCRSMatrix<B,A>() {}
49 
50  explicit BTDMatrix(int size)
51  : BCRSMatrix<B,A>(size, size, BCRSMatrix<B,A>::random)
52  {
53  // special handling for 1x1 matrices
54  if (size==1) {
55 
56  this->BCRSMatrix<B,A>::setrowsize(0, 1);
58 
59  this->BCRSMatrix<B,A>::addindex(0, 0);
61 
62  return;
63  }
64 
65  // Set number of entries for each row
66  this->BCRSMatrix<B,A>::setrowsize(0, 2);
67 
68  for (int i=1; i<size-1; i++)
69  this->BCRSMatrix<B,A>::setrowsize(i, 3);
70 
71  this->BCRSMatrix<B,A>::setrowsize(size-1, 2);
72 
74 
75  // The actual entries for each row
76  this->BCRSMatrix<B,A>::addindex(0, 0);
77  this->BCRSMatrix<B,A>::addindex(0, 1);
78 
79  for (int i=1; i<size-1; i++) {
80  this->BCRSMatrix<B,A>::addindex(i, i-1);
81  this->BCRSMatrix<B,A>::addindex(i, i );
82  this->BCRSMatrix<B,A>::addindex(i, i+1);
83  }
84 
85  this->BCRSMatrix<B,A>::addindex(size-1, size-2);
86  this->BCRSMatrix<B,A>::addindex(size-1, size-1);
87 
89 
90  }
91 
93  BTDMatrix& operator= (const BTDMatrix& other) {
94  this->BCRSMatrix<B,A>::operator=(other);
95  return *this;
96  }
97 
101  return *this;
102  }
103 
109  template <class V>
110  void solve (V& x, const V& rhs) const {
111 
112  // special handling for 1x1 matrices. The generic algorithm doesn't work for them
113  if (this->N()==1) {
114  (*this)[0][0].solve(x[0],rhs[0]);
115  return;
116  }
117 
118  // Make copies of the rhs and the right matrix band
119  V d = rhs;
120  std::vector<block_type> c(this->N()-1);
121  for (size_t i=0; i<this->N()-1; i++)
122  c[i] = (*this)[i][i+1];
123 
124  /* Modify the coefficients. */
125  block_type a_00_inv;
126  FMatrixHelp::invertMatrix((*this)[0][0], a_00_inv);
127 
128  //c[0] /= (*this)[0][0]; /* Division by zero risk. */
129  block_type c_0_tmp = c[0];
130  FMatrixHelp::multMatrix(a_00_inv, c_0_tmp, c[0]);
131 
132  // d = a^{-1} d /* Division by zero would imply a singular matrix. */
133  typename V::block_type d_0_tmp = d[0];
134  (*this)[0][0].solve(d[0], d_0_tmp);
135 
136  for (unsigned int i = 1; i < this->N(); i++) {
137 
138  // id = ( a_ii - c_{i-1} a_{i, i-1} ) ^{-1}
139  block_type tmp;
140  FMatrixHelp::multMatrix(c[i-1], (*this)[i][i-1], tmp);
141  block_type id = (*this)[i][i];
142  id -= tmp;
143  id.invert(); /* Division by zero risk. */
144 
145  if (i<c.size()) {
146  // c[i] *= id
147  tmp = c[i];
148  FMatrixHelp::multMatrix(tmp, id, c[i]); /* Last value calculated is redundant. */
149  }
150 
151  // d[i] = (d[i] - d[i-1] * (*this)[i][i-1]) * id;
152  (*this)[i][i-1].mmv(d[i-1], d[i]);
153  typename V::block_type tmpVec = d[i];
154  id.mv(tmpVec, d[i]);
155  //d[i] *= id;
156 
157  }
158 
159  /* Now back substitute. */
160  x[this->N() - 1] = d[this->N() - 1];
161  for (int i = this->N() - 2; i >= 0; i--) {
162  //x[i] = d[i] - c[i] * x[i + 1];
163  x[i] = d[i];
164  c[i].mmv(x[i+1], x[i]);
165  }
166 
167  }
168 
169 private:
170 
171  // ////////////////////////////////////////////////////////////////////////////
172  // The following methods from the base class should now actually be called
173  // ////////////////////////////////////////////////////////////////////////////
174 
175  // createbegin and createend should be in there, too, but I can't get it to compile
176  // BCRSMatrix<B,A>::CreateIterator createbegin () {}
177  // BCRSMatrix<B,A>::CreateIterator createend () {}
178  void setrowsize (size_type i, size_type s) {}
179  void incrementrowsize (size_type i) {}
180  void endrowsizes () {}
181  void addindex (size_type row, size_type col) {}
182  void endindices () {}
183 };
186 } // end namespace Dune
187 
188 #endif