Branch data Line data Source code
1 : : /*
2 : : * tridiag.cpp - tridiagonal matrix template class implementation
3 : : *
4 : : * Copyright (C) 2005, 2008 Stefan Jahn <stefan@lkcc.org>
5 : : *
6 : : * This is free software; you can redistribute it and/or modify
7 : : * it under the terms of the GNU General Public License as published by
8 : : * the Free Software Foundation; either version 2, or (at your option)
9 : : * any later version.
10 : : *
11 : : * This software is distributed in the hope that it will be useful,
12 : : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : : * GNU General Public License for more details.
15 : : *
16 : : * You should have received a copy of the GNU General Public License
17 : : * along with this package; see the file COPYING. If not, write to
18 : : * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
19 : : * Boston, MA 02110-1301, USA.
20 : : *
21 : : * $Id$
22 : : *
23 : : */
24 : :
25 : : #if HAVE_CONFIG_H
26 : : # include <config.h>
27 : : #endif
28 : :
29 : : #include <assert.h>
30 : : #include <stdio.h>
31 : : #include <stdlib.h>
32 : : #include <string.h>
33 : : #include <cmath>
34 : :
35 : : #include "compat.h"
36 : : #include "complex.h"
37 : : #include "tvector.h"
38 : :
39 : : namespace qucs {
40 : :
41 : : // Constructor creates an instance of the tridiag class.
42 : : template <class nr_type_t>
43 : 0 : tridiag<nr_type_t>::tridiag () {
44 : 0 : abov = belo = diag = offdiag = rhs = NULL;
45 : 0 : type = TRIDIAG_UNKNOWN;
46 : 0 : }
47 : :
48 : : /* The copy constructor creates a new instance based on the given
49 : : tridiag object. */
50 : : template <class nr_type_t>
51 : : tridiag<nr_type_t>::tridiag (const tridiag & t) {
52 : : abov = t.abov;
53 : : belo = t.belo;
54 : : diag = t.diag;
55 : : offdiag = t.offdiag;
56 : : rhs = t.rhs;
57 : : type = t.type;
58 : : }
59 : :
60 : : /* The assignment copy constructor creates a new instance based on the
61 : : given tridiag object. */
62 : : template <class nr_type_t>
63 : : const tridiag<nr_type_t>&
64 : : tridiag<nr_type_t>::operator=(const tridiag<nr_type_t> & t) {
65 : : if (&t != this) {
66 : : abov = t.abov;
67 : : belo = t.belo;
68 : : diag = t.diag;
69 : : offdiag = t.offdiag;
70 : : rhs = t.rhs;
71 : : type = t.type;
72 : : }
73 : : return *this;
74 : : }
75 : :
76 : : // Destructor deletes a tridiag object.
77 : : template <class nr_type_t>
78 : 0 : tridiag<nr_type_t>::~tridiag () {
79 : 0 : }
80 : :
81 : : // Set the diagonal vector of the tridiagonal matrix.
82 : : template <class nr_type_t>
83 : 0 : void tridiag<nr_type_t>::setDiagonal (tvector<nr_type_t> * v) {
84 : 0 : diag = v;
85 : 0 : }
86 : :
87 : : // Set the off-diagonal vector of the tridiagonal matrix.
88 : : template <class nr_type_t>
89 : 0 : void tridiag<nr_type_t>::setOffDiagonal (tvector<nr_type_t> * v) {
90 : 0 : abov = belo = offdiag = v;
91 : 0 : }
92 : :
93 : : // Set the above off-diagonal vector of the tridiagonal matrix.
94 : : template <class nr_type_t>
95 : : void tridiag<nr_type_t>::setA (tvector<nr_type_t> * v) {
96 : : abov = v;
97 : : }
98 : :
99 : : // Set the below off-diagonal vector of the tridiagonal matrix.
100 : : template <class nr_type_t>
101 : : void tridiag<nr_type_t>::setB (tvector<nr_type_t> * v) {
102 : : belo = v;
103 : : }
104 : :
105 : : // Set the right hand side vector of the equation system.
106 : : template <class nr_type_t>
107 : 0 : void tridiag<nr_type_t>::setRHS (tvector<nr_type_t> * v) {
108 : 0 : rhs = v;
109 : 0 : }
110 : :
111 : : /* Depending on the type of tridiagonal matrix the function runs one
112 : : of the solvers specialized on this type. The solvers are in-place
113 : : algorithms meaning the right hand side is replaced by the solution
114 : : and the original matrix entries are destroyed during solving the
115 : : system. */
116 : : template <class nr_type_t>
117 : 0 : void tridiag<nr_type_t>::solve (void) {
118 [ # # # # : 0 : switch (type) {
# ]
119 : : case TRIDIAG_NONSYM:
120 : 0 : solve_ns (); break;
121 : : case TRIDIAG_SYM:
122 : 0 : solve_s (); break;
123 : : case TRIDIAG_NONSYM_CYCLIC:
124 : 0 : solve_ns_cyc (); break;
125 : : case TRIDIAG_SYM_CYCLIC:
126 : 0 : solve_s_cyc (); break;
127 : : }
128 : 0 : }
129 : :
130 : : /* This function solves a tridiagonal equation system given
131 : : by
132 : : diag[0] abov[0] 0 ..... 0
133 : : belo[1] diag[1] abov[1] ..... 0
134 : : 0 belo[2] diag[2]
135 : : 0 0 belo[3] ..... abov[n-2]
136 : : ... ...
137 : : 0 ... belo[n-1] diag[n-1]
138 : : */
139 : : template <class nr_type_t>
140 : 0 : void tridiag<nr_type_t>::solve_ns (void) {
141 : 0 : d = al = diag->getData ();
142 : 0 : f = ga = abov->getData ();
143 : 0 : e = belo->getData ();
144 : 0 : b = c = x = rhs->getData ();
145 : 0 : int i, n = diag->getSize ();
146 : :
147 : : // factorize A = LU
148 : 0 : al[0] = d[0];
149 : 0 : ga[0] = f[0] / al[0];
150 [ # # ]: 0 : for (i = 1; i < n - 1; i++) {
151 : 0 : al[i] = d[i] - e[i] * ga[i-1];
152 : 0 : ga[i] = f[i] / al[i];
153 : : }
154 : 0 : al[n-1] = d[n-1] - e[n-1] * ga[n-2];
155 : :
156 : : // update b = Lc
157 : 0 : c[0] = b[0] / d[0];
158 [ # # ]: 0 : for (i = 1; i < n; i++) {
159 : 0 : c[i] = (b[i] - e[i] * c[i-1]) / al[i];
160 : : }
161 : :
162 : : // back substitution Rx = c
163 : 0 : x[n-1] = c[n-1];
164 [ # # ]: 0 : for (i = n - 2; i >= 0; i--) {
165 : 0 : x[i] = c[i] - ga[i] * x[i+1];
166 : : }
167 : 0 : }
168 : :
169 : : /* This function solves a cyclically tridiagonal equation system given
170 : : by
171 : : diag[0] abov[0] 0 ..... belo[0]
172 : : belo[1] diag[1] abov[1] ..... 0
173 : : 0 belo[2] diag[2]
174 : : 0 0 belo[3] ..... abov[n-2]
175 : : ... ...
176 : : abov[n-1] ... belo[n-1] diag[n-1]
177 : : */
178 : : template <class nr_type_t>
179 : 0 : void tridiag<nr_type_t>::solve_ns_cyc (void) {
180 : 0 : d = al = diag->getData ();
181 : 0 : f = ga = abov->getData ();
182 : 0 : e = be = belo->getData ();
183 : 0 : b = x = c = rhs->getData ();
184 : 0 : int i, n = diag->getSize ();
185 : 0 : de = new nr_type_t[n];
186 : 0 : ep = new nr_type_t[n];
187 : :
188 : : // factorize A = LU
189 : 0 : al[0] = d[0];
190 : 0 : ga[0] = f[0] / al[0];
191 : 0 : de[0] = e[0] / al[0];
192 [ # # ]: 0 : for (i = 1; i < n - 2; i++) {
193 : 0 : al[i] = d[i] - e[i] * ga[i-1];
194 : 0 : ga[i] = f[i] / al[i];
195 : 0 : be[i] = e[i];
196 : 0 : de[i] = -be[i] * de[i-1] / al[i];
197 : : }
198 : 0 : al[n-2] = d[n-2] - e[n-2] * ga[n-3];
199 : 0 : be[n-2] = e[n-2];
200 : 0 : ep[2] = f[n-1];
201 [ # # ]: 0 : for (i = 3; i < n; i++) {
202 : 0 : ep[i] = -ep[i-1] * ga[i-3];
203 : : }
204 : 0 : ga[n-2] = (f[n-2] - be[n-2] * de[n-3]) / al[n-2];
205 : 0 : be[n-1] = e[n-1] - ep[n-1] * ga[n-3];
206 : 0 : al[n-1] = d[n-1] - be[n-1] * ga[n-2];
207 [ # # ]: 0 : for (i = 2; i < n; i++) {
208 : 0 : al[n-1] -= ep[i] * de[i-2];
209 : : }
210 : :
211 : : // update Lc = b
212 : 0 : c[0] = b[0] / al[0];
213 [ # # ]: 0 : for (i = 1; i < n - 1; i++) {
214 : 0 : c[i] = (b[i] - c[i-1] * be[i]) / al[i];
215 : : }
216 : 0 : c[n-1] = b[n-1] - be[n-1] * c[n-2];
217 [ # # ]: 0 : for (i = 2; i < n; i++) {
218 : 0 : c[n-1] -= ep[i] * c[i-2];
219 : : }
220 : 0 : c[n-1] /= al[n-1];
221 : :
222 : : // back substitution Rx = c
223 : 0 : x[n-1] = c[n-1];
224 : 0 : x[n-2] = c[n-2] - ga[n-2] * x[n-1];
225 [ # # ]: 0 : for (i = n - 3; i >= 0; i--) {
226 : 0 : x[i] = c[i] - ga[i] * x[i+1] - de[i] * x[n-1];
227 : : }
228 : :
229 [ # # ]: 0 : delete[] de;
230 [ # # ]: 0 : delete[] ep;
231 : 0 : }
232 : :
233 : : /* This function solves a symmetric tridiagonal strongly nonsingular
234 : : equation system given by
235 : : diag[0] offdiag[0] 0 ..... 0
236 : : offdiag[0] diag[1] offdiag[1] ..... 0
237 : : 0 offdiag[1] diag[2]
238 : : 0 0 offdiag[2] ..... offdiag[n-2]
239 : : ... ...
240 : : 0 ... offdiag[n-2] diag[n-1]
241 : : */
242 : : template <class nr_type_t>
243 : 0 : void tridiag<nr_type_t>::solve_s (void) {
244 : 0 : d = al = diag->getData ();
245 : 0 : f = ga = offdiag->getData ();
246 : 0 : b = z = x = c = rhs->getData ();
247 : : nr_type_t t;
248 : 0 : int i, n = diag->getSize ();
249 : 0 : de = new nr_type_t[n];
250 : :
251 : : // factorize A = LDL'
252 : 0 : al[0] = d[0];
253 : 0 : t = f[0];
254 : 0 : ga[0] = t / al[0];
255 [ # # ]: 0 : for (i = 1; i < n - 1; i++) {
256 : 0 : al[i] = d[i] - t * ga[i-1];
257 : 0 : t = f[i];
258 : 0 : ga[i] = t / al[i];
259 : : }
260 : 0 : al[n-1] = d[n-1] - t * ga[n-2];
261 : :
262 : : // update L'z = b and Dc = z
263 : 0 : z[0] = b[0];
264 [ # # ]: 0 : for (i = 1; i < n; i++) {
265 : 0 : z[i] = b[i] - ga[i-1] * z[i-1];
266 : : }
267 [ # # ]: 0 : for (i = 0; i < n; i++) {
268 : 0 : c[i] = z[i] / al[i];
269 : : }
270 : :
271 : : // back substitution L'x = c
272 : 0 : x[n-1] = c[n-1];
273 [ # # ]: 0 : for (i = n-2; i >= 0; i--) {
274 : 0 : x[i] = c[i] - ga[i] * x[i+1];
275 : : }
276 : :
277 [ # # ]: 0 : delete[] de;
278 : 0 : }
279 : :
280 : : /* This function solves a symmetric cyclically tridiagonal strongly
281 : : nonsingular equation system given by
282 : : diag[0] offdiag[0] 0 ..... offdiag[n-1]
283 : : offdiag[0] diag[1] offdiag[1] ..... 0
284 : : 0 offdiag[1] diag[2]
285 : : 0 0 offdiag[2] ..... offdiag[n-2]
286 : : ... ...
287 : : offdiag[n-1] ... offdiag[n-2] diag[n-1]
288 : : */
289 : : template <class nr_type_t>
290 : 0 : void tridiag<nr_type_t>::solve_s_cyc (void) {
291 : 0 : d = al = diag->getData ();
292 : 0 : f = ga = offdiag->getData ();
293 : 0 : b = c = z = x = rhs->getData ();
294 : : nr_type_t t;
295 : 0 : int i, n = diag->getSize ();
296 : 0 : de = new nr_type_t[n];
297 : :
298 : : // factorize A = LDL'
299 : 0 : al[0] = d[0];
300 : 0 : t = f[0];
301 : 0 : ga[0] = t / al[0];
302 : 0 : de[0] = f[n-1] / al[0];
303 [ # # ]: 0 : for (i = 1; i < n - 2; i++) {
304 : 0 : al[i] = d[i] - t * ga[i-1];
305 : 0 : de[i] = -de[i-1] * t / al[i];
306 : 0 : t = f[i];
307 : 0 : ga[i] = t / al[i];
308 : : }
309 : 0 : al[n-2] = d[n-2] - t * ga[n-3];
310 : 0 : ga[n-2] = (f[n-2] - t * de[n-3]) / al[n-2];
311 : 0 : al[n-1] = d[n-1] - al[n-2] * ga[n-2] * ga[n-2];
312 [ # # ]: 0 : for (i = 0; i < n - 2; i++) {
313 : 0 : al[n-1] -= al[i] * de[i] * de[i];
314 : : }
315 : :
316 : : // update Lz = b and Dc = z
317 : 0 : z[0] = b[0];
318 [ # # ]: 0 : for (i = 1; i < n - 1; i++) {
319 : 0 : z[i] = b[i] - ga[i-1] * z[i-1];
320 : : }
321 : 0 : z[n-1] = b[n-1] - ga[n-2] * z[n-2];
322 [ # # ]: 0 : for (i = 0; i < n - 2; i++) {
323 : 0 : z[n-1] -= de[i] * z[i];
324 : : }
325 [ # # ]: 0 : for (i = 0; i < n; i++) {
326 : 0 : c[i] = z[i] / al[i];
327 : : }
328 : :
329 : : // back substitution L'x = c
330 : 0 : x[n-1] = c[n-1];
331 : 0 : x[n-2] = c[n-2] - ga[n-2] * x[n-1];
332 [ # # ]: 0 : for (i = n - 3; i >= 0; i--) {
333 : 0 : x[i] = c[i] - ga[i] * x[i+1] - de[i] * x[n-1];
334 : : }
335 : :
336 [ # # ]: 0 : delete[] de;
337 : 0 : }
338 : :
339 : : } // namespace qucs
|
