Reference documentation for deal.II version 9.5.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\}}\)
Loading...
Searching...
No Matches
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
25using namespace LAPACKSupport;
26
27template <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
38template <typename number>
39void
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
51template <typename number>
52bool
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
79template <typename number>
80void
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
130template <typename number>
131void
133 const Vector<number> &v) const
134{
135 vmult(w, v, /*adding = */ true);
136}
137
138
139
140template <typename number>
141void
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
185template <typename number>
186void
188 const Vector<number> &v) const
189{
190 Tvmult(w, v, true);
191}
192
193
194
195template <typename number>
196number
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
222template <typename number>
223number
225{
226 return matrix_scalar_product(v, v);
227}
228
229
230
231template <typename number>
232void
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
259template <typename number>
260number
262{
264 AssertIndexRange(i, n());
265 return diagonal[i];
266}
267
268
269
270template class TridiagonalMatrix<float>;
271template 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
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
size_type size() const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcNeedsLAPACK()
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcState(State arg1)
#define AssertThrow(cond, exc)
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
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)