Reference documentation for deal.II version GIT e22cb6a53e 2023-09-21 14:00:02+00:00
\(\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 - 2022 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 {
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 {
55  Assert(state == matrix, ExcState(state));
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 {
85  Assert(state == matrix, ExcState(state));
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 and last row
94  const size_type e = n() - 1;
95  // Let iterators point to the first entry of each diagonal
96  typename std::vector<number>::const_iterator d = diagonal.begin();
97  typename std::vector<number>::const_iterator r = right.begin();
98  // The left diagonal starts one later or is equal to the right
99  // one for symmetric storage
100  typename std::vector<number>::const_iterator l = left.begin();
101  if (is_symmetric)
102  l = r;
103  else
104  ++l;
105 
106  if (adding)
107  {
108  // Treat first row separately
109  w(0) += (*d) * v(0) + (*r) * v(1);
110  ++d;
111  ++r;
112  // All rows with three entries
113  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
114  w(i) += (*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1);
115  // Last row is special again
116  w(e) += (*l) * v(e - 1) + (*d) * v(e);
117  }
118  else
119  {
120  w(0) = (*d) * v(0) + (*r) * v(1);
121  ++d;
122  ++r;
123  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
124  w(i) = (*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1);
125  w(e) = (*l) * v(e - 1) + (*d) * v(e);
126  }
127 }
128 
129 
130 template <typename number>
131 void
133  const Vector<number> &v) const
134 {
135  vmult(w, v, /*adding = */ true);
136 }
137 
138 
139 
140 template <typename number>
141 void
143  const Vector<number> &v,
144  const bool adding) const
145 {
146  Assert(state == matrix, ExcState(state));
147 
148  Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n()));
149  Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
150 
151  if (n() == 0)
152  return;
153 
154  const size_type e = n() - 1;
155  typename std::vector<number>::const_iterator d = diagonal.begin();
156  typename std::vector<number>::const_iterator r = right.begin();
157  typename std::vector<number>::const_iterator l = left.begin();
158  if (is_symmetric)
159  l = r;
160  else
161  ++l;
162 
163  if (adding)
164  {
165  w(0) += (*d) * v(0) + (*l) * v(1);
166  ++d;
167  ++l;
168  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
169  w(i) += (*l) * v(i + 1) + (*d) * v(i) + (*r) * v(i - 1);
170  w(e) += (*d) * v(e) + (*r) * v(e - 1);
171  }
172  else
173  {
174  w(0) = (*d) * v(0) + (*l) * v(1);
175  ++d;
176  ++l;
177  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
178  w(i) = (*l) * v(i + 1) + (*d) * v(i) + (*r) * v(i - 1);
179  w(e) = (*d) * v(e) + (*r) * v(e - 1);
180  }
181 }
182 
183 
184 
185 template <typename number>
186 void
188  const Vector<number> &v) const
189 {
190  Tvmult(w, v, true);
191 }
192 
193 
194 
195 template <typename number>
196 number
198  const Vector<number> &v) const
199 {
200  Assert(state == matrix, ExcState(state));
201 
202  const size_type e = n() - 1;
203  typename std::vector<number>::const_iterator d = diagonal.begin();
204  typename std::vector<number>::const_iterator r = right.begin();
205  typename std::vector<number>::const_iterator l = left.begin();
206  if (is_symmetric)
207  l = r;
208  else
209  ++l;
210 
211  number result = w(0) * ((*d) * v(0) + (*r) * v(1));
212  ++d;
213  ++r;
214  for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
215  result += w(i) * ((*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1));
216  result += w(e) * ((*l) * v(e - 1) + (*d) * v(e));
217  return result;
218 }
219 
220 
221 
222 template <typename number>
223 number
225 {
226  return matrix_scalar_product(v, v);
227 }
228 
229 
230 
231 template <typename number>
232 void
234 {
235 #ifdef DEAL_II_WITH_LAPACK
236  Assert(state == matrix, ExcState(state));
237  Assert(is_symmetric, ExcNotImplemented());
238 
239  const types::blas_int nn = n();
240  types::blas_int info;
241  stev(&N,
242  &nn,
243  diagonal.data(),
244  right.data(),
245  static_cast<number *>(nullptr),
246  &one,
247  static_cast<number *>(nullptr),
248  &info);
249  Assert(info == 0, ExcInternalError());
250 
252 #else
253  AssertThrow(false, ExcNeedsLAPACK());
254 #endif
255 }
256 
257 
258 
259 template <typename number>
260 number
262 {
263  Assert(state == LAPACKSupport::eigenvalues, ExcState(state));
264  AssertIndexRange(i, n());
265  return diagonal[i];
266 }
267 
268 
269 
270 template class TridiagonalMatrix<float>;
271 template class TridiagonalMatrix<double>;
272 #ifdef DEAL_II_WITH_COMPLEX_VALUES
275 #endif
276 
void Tvmult_add(Vector< number > &w, const Vector< number > &v) const
void vmult(Vector< number > &w, const Vector< number > &v, const bool adding=false) const
number matrix_norm_square(const Vector< number > &v) const
void reinit(size_type n, bool symmetric=false)
void vmult_add(Vector< number > &w, const Vector< number > &v) const
types::global_dof_index size_type
void Tvmult(Vector< number > &w, const Vector< number > &v, const bool adding=false) const
number eigenvalue(const size_type i) const
TridiagonalMatrix(size_type n=0, bool symmetric=false)
number matrix_scalar_product(const Vector< number > &u, const Vector< number > &v) const
Definition: vector.h:110
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsLAPACK()
static ::ExceptionBase & ExcState(State arg1)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcNotImplemented()
#define AssertIndexRange(index, range)
Definition: exceptions.h:1857
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
void stev(const char *, const ::types::blas_int *, number1 *, number2 *, number3 *, const ::types::blas_int *, number4 *, ::types::blas_int *)
@ matrix
Contents is actually a matrix.
@ eigenvalues
Eigenvalue vector is filled.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
static const char N
static const types::blas_int one
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
int blas_int