Reference documentation for deal.II version GIT relicensing-14-gd371fb467a 2024-02-26 20:50:01+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\}}\)
Loading...
Searching...
No Matches
tridiagonal_matrix.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2005 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
15
18#include <deal.II/lac/vector.h>
19
20#include <complex>
21
23
24using namespace LAPACKSupport;
25
26template <typename number>
28 : diagonal(size, 0.)
29 , left((symmetric ? 0 : size), 0.)
30 , right(size, 0.)
31 , is_symmetric(symmetric)
32 , state(matrix)
33{}
34
35
36
37template <typename number>
38void
40{
41 is_symmetric = symmetric;
42 diagonal.resize(size);
43 right.resize(size);
44 left.resize(symmetric ? 0 : size);
45 state = matrix;
46}
47
48
49
50template <typename number>
51bool
53{
54 Assert(state == matrix, ExcState(state));
55
56 typename std::vector<number>::const_iterator i;
57 typename std::vector<number>::const_iterator e;
58
59 e = diagonal.end();
60 for (i = diagonal.begin(); i != e; ++i)
61 if (std::abs(*i) != 0.)
62 return false;
63
64 e = left.end();
65 for (i = left.begin(); i != e; ++i)
66 if (std::abs(*i) != 0.)
67 return false;
68
69 e = right.end();
70 for (i = right.begin(); i != e; ++i)
71 if (std::abs(*i) != 0.)
72 return false;
73 return true;
74}
75
76
77
78template <typename number>
79void
81 const Vector<number> &v,
82 const bool adding) const
83{
84 Assert(state == matrix, ExcState(state));
85
86 Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n()));
87 Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
88
89 if (n() == 0)
90 return;
91
92 // The actual loop skips the first and last row
93 const size_type e = n() - 1;
94 // Let iterators point to the first entry of each diagonal
95 typename std::vector<number>::const_iterator d = diagonal.begin();
96 typename std::vector<number>::const_iterator r = right.begin();
97 // The left diagonal starts one later or is equal to the right
98 // one for symmetric storage
99 typename std::vector<number>::const_iterator l = left.begin();
100 if (is_symmetric)
101 l = r;
102 else
103 ++l;
104
105 if (adding)
106 {
107 // Treat first row separately
108 w(0) += (*d) * v(0) + (*r) * v(1);
109 ++d;
110 ++r;
111 // All rows with three entries
112 for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
113 w(i) += (*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1);
114 // Last row is special again
115 w(e) += (*l) * v(e - 1) + (*d) * v(e);
116 }
117 else
118 {
119 w(0) = (*d) * v(0) + (*r) * v(1);
120 ++d;
121 ++r;
122 for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
123 w(i) = (*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1);
124 w(e) = (*l) * v(e - 1) + (*d) * v(e);
125 }
126}
127
128
129template <typename number>
130void
132 const Vector<number> &v) const
133{
134 vmult(w, v, /*adding = */ true);
135}
136
137
138
139template <typename number>
140void
142 const Vector<number> &v,
143 const bool adding) const
144{
145 Assert(state == matrix, ExcState(state));
146
147 Assert(w.size() == n(), ExcDimensionMismatch(w.size(), n()));
148 Assert(v.size() == n(), ExcDimensionMismatch(v.size(), n()));
149
150 if (n() == 0)
151 return;
152
153 const size_type e = n() - 1;
154 typename std::vector<number>::const_iterator d = diagonal.begin();
155 typename std::vector<number>::const_iterator r = right.begin();
156 typename std::vector<number>::const_iterator l = left.begin();
157 if (is_symmetric)
158 l = r;
159 else
160 ++l;
161
162 if (adding)
163 {
164 w(0) += (*d) * v(0) + (*l) * v(1);
165 ++d;
166 ++l;
167 for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
168 w(i) += (*l) * v(i + 1) + (*d) * v(i) + (*r) * v(i - 1);
169 w(e) += (*d) * v(e) + (*r) * v(e - 1);
170 }
171 else
172 {
173 w(0) = (*d) * v(0) + (*l) * v(1);
174 ++d;
175 ++l;
176 for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
177 w(i) = (*l) * v(i + 1) + (*d) * v(i) + (*r) * v(i - 1);
178 w(e) = (*d) * v(e) + (*r) * v(e - 1);
179 }
180}
181
182
183
184template <typename number>
185void
187 const Vector<number> &v) const
188{
189 Tvmult(w, v, true);
190}
191
192
193
194template <typename number>
195number
197 const Vector<number> &v) const
198{
199 Assert(state == matrix, ExcState(state));
200
201 const size_type e = n() - 1;
202 typename std::vector<number>::const_iterator d = diagonal.begin();
203 typename std::vector<number>::const_iterator r = right.begin();
204 typename std::vector<number>::const_iterator l = left.begin();
205 if (is_symmetric)
206 l = r;
207 else
208 ++l;
209
210 number result = w(0) * ((*d) * v(0) + (*r) * v(1));
211 ++d;
212 ++r;
213 for (size_type i = 1; i < e; ++i, ++d, ++r, ++l)
214 result += w(i) * ((*l) * v(i - 1) + (*d) * v(i) + (*r) * v(i + 1));
215 result += w(e) * ((*l) * v(e - 1) + (*d) * v(e));
216 return result;
217}
218
219
220
221template <typename number>
222number
224{
225 return matrix_scalar_product(v, v);
226}
227
228
229
230template <typename number>
231void
233{
234#ifdef DEAL_II_WITH_LAPACK
235 Assert(state == matrix, ExcState(state));
236 Assert(is_symmetric, ExcNotImplemented());
237
238 const types::blas_int nn = n();
239 types::blas_int info;
240 stev(&N,
241 &nn,
242 diagonal.data(),
243 right.data(),
244 static_cast<number *>(nullptr),
245 &one,
246 static_cast<number *>(nullptr),
247 &info);
248 Assert(info == 0, ExcInternalError());
249
251#else
252 AssertThrow(false, ExcNeedsLAPACK());
253#endif
254}
255
256
257
258template <typename number>
259number
261{
263 AssertIndexRange(i, n());
264 return diagonal[i];
265}
266
267
268
269template class TridiagonalMatrix<float>;
270template class TridiagonalMatrix<double>;
271#ifdef DEAL_II_WITH_COMPLEX_VALUES
274#endif
275
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
virtual size_type size() const override
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:502
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:503
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 > &)