Reference documentation for deal.II version 9.2.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tridiagonal_matrix.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2020 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 
19 #include <deal.II/lac/vector.h>
20 
21 #include <complex>
22 
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 {
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  AssertThrow(false, ExcNeedsLAPACK());
257 #endif
258 }
259 
260 
261 
262 template <typename number>
263 number
265 {
267  AssertIndexRange(i, 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 
void stev(const char *, const ::types::blas_int *, number1 *, number2 *, number3 *, const ::types::blas_int *, number4 *, ::types::blas_int *)
Matrix is symmetric.
Contents is actually a matrix.
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
std::vector< number > left
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
static const types::blas_int one
void vmult(Vector< number > &w, const Vector< number > &v, const bool adding=false) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
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:1419
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
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
TridiagonalMatrix(size_type n=0, bool symmetric=false)
number matrix_norm_square(const Vector< number > &v) const
size_type n() const
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
number eigenvalue(const size_type i) const
void Tvmult(Vector< number > &w, const Vector< number > &v, const bool adding=false) const
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
LAPACKSupport::State state
static const char N
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
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()