Reference documentation for deal.II version Git 7026f387cc 2019-10-15 14:19:01 -0400
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
tridiagonal_matrix.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 #include <deal.II/lac/lapack_templates.h>
18 #include <deal.II/lac/tridiagonal_matrix.h>
19 #include <deal.II/lac/vector.h>
20 
21 #include <complex>
22 
23 DEAL_II_NAMESPACE_OPEN
24 
25 using namespace LAPACKSupport;
26 
27 template <typename number>
29  : diagonal(size, 0.)
30  , left((symmetric ? 0 : size), 0.)
31  , right(size, 0.)
32  , is_symmetric(symmetric)
33  , state(matrix)
34 {}
35 
36 
37 
38 template <typename number>
39 void
41 {
42  is_symmetric = symmetric;
43  diagonal.resize(size);
44  right.resize(size);
45  left.resize(symmetric ? 0 : size);
46  state = matrix;
47 }
48 
49 
50 
51 template <typename number>
52 bool
54 {
56 
57  typename std::vector<number>::const_iterator i;
58  typename std::vector<number>::const_iterator e;
59 
60  e = diagonal.end();
61  for (i = diagonal.begin(); i != e; ++i)
62  if (std::abs(*i) != 0.)
63  return false;
64 
65  e = left.end();
66  for (i = left.begin(); i != e; ++i)
67  if (std::abs(*i) != 0.)
68  return false;
69 
70  e = right.end();
71  for (i = right.begin(); i != e; ++i)
72  if (std::abs(*i) != 0.)
73  return false;
74  return true;
75 }
76 
77 
78 
79 template <typename number>
80 void
82  const Vector<number> &v,
83  const bool adding) const
84 {
86 
87  Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n()));
88  Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
89 
90  if (n() == 0)
91  return;
92 
93  // The actual loop skips the first
94  // and last row
95  const size_type e = n() - 1;
96  // Let iterators point to the first
97  // entry of each diagonal
98  typename std::vector<number>::const_iterator d = diagonal.begin();
99  typename std::vector<number>::const_iterator r = right.begin();
100  // The left diagonal starts one
101  // later or is equal to the right
102  // one for symmetric storage
103  typename std::vector<number>::const_iterator l = left.begin();
104  if (is_symmetric)
105  l = r;
106  else
107  ++l;
108 
109  if (adding)
110  {
111  // Treat first row separately
112  w(0) += (*d) * v(0) + (*r) * v(1);
113  ++d;
114  ++r;
115  // All rows with three entries
116  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
117  w(i) += (*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1);
118  // Last row is special again
119  w(e) += (*l) * v(e - 1) + (*d) * v(e);
120  }
121  else
122  {
123  w(0) = (*d) * v(0) + (*r) * v(1);
124  ++d;
125  ++r;
126  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
127  w(i) = (*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1);
128  w(e) = (*l) * v(e - 1) + (*d) * v(e);
129  }
130 }
131 
132 
133 template <typename number>
134 void
136  const Vector<number> &v) const
137 {
138  vmult(w, v, true);
139 }
140 
141 
142 
143 template <typename number>
144 void
146  const Vector<number> &v,
147  const bool adding) const
148 {
150 
151  Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n()));
152  Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
153 
154  if (n() == 0)
155  return;
156 
157  const size_type e = n() - 1;
158  typename std::vector<number>::const_iterator d = diagonal.begin();
159  typename std::vector<number>::const_iterator r = right.begin();
160  typename std::vector<number>::const_iterator l = left.begin();
161  if (is_symmetric)
162  l = r;
163  else
164  ++l;
165 
166  if (adding)
167  {
168  w(0) += (*d) * v(0) + (*l) * v(1);
169  ++d;
170  ++l;
171  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
172  w(i) += (*l) * v(i + 1) + (*d) * v(i) + (*r) * v(i - 1);
173  w(e) += (*d) * v(e) + (*r) * v(e - 1);
174  }
175  else
176  {
177  w(0) = (*d) * v(0) + (*l) * v(1);
178  ++d;
179  ++l;
180  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
181  w(i) = (*l) * v(i + 1) + (*d) * v(i) + (*r) * v(i - 1);
182  w(e) = (*d) * v(e) + (*r) * v(e - 1);
183  }
184 }
185 
186 
187 
188 template <typename number>
189 void
191  const Vector<number> &v) const
192 {
193  Tvmult(w, v, true);
194 }
195 
196 
197 
198 template <typename number>
199 number
201  const Vector<number> &v) const
202 {
204 
205  const size_type e = n() - 1;
206  typename std::vector<number>::const_iterator d = diagonal.begin();
207  typename std::vector<number>::const_iterator r = right.begin();
208  typename std::vector<number>::const_iterator l = left.begin();
209  if (is_symmetric)
210  l = r;
211  else
212  ++l;
213 
214  number result = w(0) * ((*d) * v(0) + (*r) * v(1));
215  ++d;
216  ++r;
217  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
218  result += w(i) * ((*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1));
219  result += w(e) * ((*l) * v(e - 1) + (*d) * v(e));
220  return result;
221 }
222 
223 
224 
225 template <typename number>
226 number
228 {
229  return matrix_scalar_product(v, v);
230 }
231 
232 
233 
234 template <typename number>
235 void
237 {
238 #ifdef DEAL_II_WITH_LAPACK
241 
242  const types::blas_int nn = n();
243  types::blas_int info;
244  stev(&N,
245  &nn,
246  diagonal.data(),
247  right.data(),
248  static_cast<number *>(nullptr),
249  &one,
250  static_cast<number *>(nullptr),
251  &info);
252  Assert(info == 0, ExcInternalError());
253 
255 #else
256  Assert(false, ExcNeedsLAPACK());
257 #endif
258 }
259 
260 
261 
262 template <typename number>
263 number
265 {
267  Assert(i < n(), ExcIndexRange(i, 0, n()));
268  return diagonal[i];
269 }
270 
271 
272 /*
273 template <typename number>
274 TridiagonalMatrix<number>::
275 {
276 }
277 
278 
279 */
280 
281 template class TridiagonalMatrix<float>;
282 template class TridiagonalMatrix<double>;
283 #ifdef DEAL_II_WITH_COMPLEX_VALUES
286 #endif
287 
288 DEAL_II_NAMESPACE_CLOSE
Matrix is symmetric.
Contents is actually a matrix.
std::vector< number > left
void vmult(Vector< number > &w, const Vector< number > &v, const bool adding=false) const
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void Tvmult_add(Vector< number > &w, const Vector< number > &v) const
Matrix is diagonal.
static ::ExceptionBase & ExcNeedsLAPACK()
static ::ExceptionBase & ExcState(State arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1407
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::vector< number > right
number matrix_scalar_product(const Vector< number > &u, const Vector< number > &v) const
TridiagonalMatrix(size_type n=0, bool symmetric=false)
number matrix_norm_square(const Vector< number > &v) const
size_type n() const
number eigenvalue(const size_type i) const
void Tvmult(Vector< number > &w, const Vector< number > &v, const bool adding=false) const
int blas_int
LAPACKSupport::State state
size_type size() const
std::vector< number > diagonal
static ::ExceptionBase & ExcNotImplemented()
void reinit(size_type n, bool symmetric=false)
Eigenvalue vector is filled.
void vmult_add(Vector< number > &w, const Vector< number > &v) const
types::global_dof_index size_type
static ::ExceptionBase & ExcInternalError()