Branch data Line data Source code
1 : : /*
2 : : * tmatrix.cpp - simple matrix template class implementation
3 : : *
4 : : * Copyright (C) 2004, 2005, 2006, 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 : : #include "qucs_typedefs.h"
26 : :
27 : : #include <assert.h>
28 : : #include <stdio.h>
29 : : #include <stdlib.h>
30 : : #include <string.h>
31 : : #include <cmath>
32 : :
33 : : #include "compat.h"
34 : : #include "logging.h"
35 : : #include "complex.h"
36 : : #include "tmatrix.h"
37 : :
38 : : namespace qucs {
39 : :
40 : : // Constructor creates an unnamed instance of the tmatrix class.
41 : : template <class nr_type_t>
42 : 0 : tmatrix<nr_type_t>::tmatrix () {
43 : 0 : rows = 0;
44 : 0 : cols = 0;
45 : 0 : data = NULL;
46 : 0 : }
47 : :
48 : : /* Constructor creates an unnamed instance of the tmatrix class with a
49 : : certain number of rows and columns. Creates a square tmatrix. */
50 : : template <class nr_type_t>
51 : 454776 : tmatrix<nr_type_t>::tmatrix (int s) {
52 : 454776 : rows = cols = s;
53 [ + - ]: 454776 : if (s > 0) {
54 [ + + ]: 697660 : data = new nr_type_t[s * s];
55 : 454776 : memset (data, 0, sizeof (nr_type_t) * s * s);
56 : : }
57 : 0 : else data = NULL;
58 : 454776 : }
59 : :
60 : : /* Constructor creates an unnamed instance of the tmatrix class with a
61 : : certain number of rows and columns. */
62 : : template <class nr_type_t>
63 : : tmatrix<nr_type_t>::tmatrix (int r, int c) {
64 : : rows = r;
65 : : cols = c;
66 : : if (r > 0 && c > 0) {
67 : : data = new nr_type_t[r * c];
68 : : memset (data, 0, sizeof (nr_type_t) * r * c);
69 : : }
70 : : else data = NULL;
71 : : }
72 : :
73 : : /* The copy constructor creates a new instance based on the given
74 : : tmatrix object. */
75 : : template <class nr_type_t>
76 : 25685 : tmatrix<nr_type_t>::tmatrix (const tmatrix & m) {
77 : 25685 : rows = m.rows;
78 : 25685 : cols = m.cols;
79 : 25685 : data = NULL;
80 : :
81 : : // copy tmatrix elements
82 [ + - ][ + - ]: 25685 : if (rows > 0 && cols > 0) {
83 [ + + ]: 2042821 : data = new nr_type_t[rows * cols];
84 : 25685 : memcpy (data, m.data, sizeof (nr_type_t) * rows * cols);
85 : : }
86 : 25685 : }
87 : :
88 : : /* The assignment copy constructor creates a new instance based on the
89 : : given tmatrix object. */
90 : : template <class nr_type_t>
91 : : const tmatrix<nr_type_t>&
92 : 2 : tmatrix<nr_type_t>::operator=(const tmatrix<nr_type_t> & m) {
93 [ + - ]: 2 : if (&m != this) {
94 : 2 : rows = m.rows;
95 : 2 : cols = m.cols;
96 [ + - ][ + - ]: 2 : if (data) { delete[] data; data = NULL; }
97 [ + - ][ + - ]: 2 : if (rows > 0 && cols > 0) {
98 [ + + ]: 3174 : data = new nr_type_t[rows * cols];
99 : 2 : memcpy (data, m.data, sizeof (nr_type_t) * rows * cols);
100 : : }
101 : : }
102 : 2 : return *this;
103 : : }
104 : :
105 : : // Destructor deletes a tmatrix object.
106 : : template <class nr_type_t>
107 : 480461 : tmatrix<nr_type_t>::~tmatrix () {
108 [ + - ][ + - ]: 480461 : if (data) delete[] data;
109 : 480461 : }
110 : :
111 : : // Returns the tmatrix element at the given row and column.
112 : : template <class nr_type_t>
113 : 2212918 : nr_type_t tmatrix<nr_type_t>::get (int r, int c) {
114 [ + - ][ + - ]: 2212918 : assert (r >= 0 && r < rows && c >= 0 && c < cols);
[ + - ][ - + ]
[ - + ]
115 : 2212918 : return data[r * cols + c];
116 : : }
117 : :
118 : : // Sets the tmatrix element at the given row and column.
119 : : template <class nr_type_t>
120 : 58758696 : void tmatrix<nr_type_t>::set (int r, int c, nr_type_t z) {
121 [ + - ][ + - ]: 58758696 : assert (r >= 0 && r < rows && c >= 0 && c < cols);
[ + - ][ - + ]
[ - + ]
122 : 58758696 : data[r * cols + c] = z;
123 : 58758696 : }
124 : :
125 : : // Sets all the tmatrix elements to the given value.
126 : : template <class nr_type_t>
127 : 22 : void tmatrix<nr_type_t>::set (nr_type_t z) {
128 [ + + ]: 5654 : for (int i = 0; i < rows * cols; i++) data[i] = z;
129 : 22 : }
130 : :
131 : : // The function returns the given row in a tvector.
132 : : template <class nr_type_t>
133 : : tvector<nr_type_t> tmatrix<nr_type_t>::getRow (int r) {
134 : : assert (r >= 0 && r < rows);
135 : : tvector<nr_type_t> res (cols);
136 : : nr_type_t * dst = res.getData ();
137 : : nr_type_t * src = &data[r * cols];
138 : : memcpy (dst, src, sizeof (nr_type_t) * cols);
139 : : return res;
140 : : }
141 : :
142 : : // Puts the given tvector into the given row of the tmatrix instance.
143 : : template <class nr_type_t>
144 : : void tmatrix<nr_type_t>::setRow (int r, tvector<nr_type_t> v) {
145 : : assert (r >= 0 && r < rows && v.getSize () == cols);
146 : : nr_type_t * dst = &data[r * cols];
147 : : nr_type_t * src = v.getData ();
148 : : memcpy (dst, src, sizeof (nr_type_t) * cols);
149 : : }
150 : :
151 : : // The function returns the given column in a tvector.
152 : : template <class nr_type_t>
153 : : tvector<nr_type_t> tmatrix<nr_type_t>::getCol (int c) {
154 : : assert (c >= 0 && c < cols);
155 : : tvector<nr_type_t> res (rows);
156 : : nr_type_t * dst = res.getData ();
157 : : nr_type_t * src = &data[c];
158 : : for (int r = 0; r < rows; r++, src += cols, dst++) *dst = *src;
159 : : return res;
160 : : }
161 : :
162 : : // Puts the given tvector into the given column of the tmatrix instance.
163 : : template <class nr_type_t>
164 : : void tmatrix<nr_type_t>::setCol (int c, tvector<nr_type_t> v) {
165 : : assert (c >= 0 && c < cols && v.getSize () == rows);
166 : : nr_type_t * dst = &data[c];
167 : : nr_type_t * src = v.getData ();
168 : : for (int r = 0; r < rows; r++, src++, dst += cols) *dst = *src;
169 : : }
170 : :
171 : : // The function swaps the given rows with each other.
172 : : template <class nr_type_t>
173 : 2224451 : void tmatrix<nr_type_t>::exchangeRows (int r1, int r2) {
174 [ + - ][ + - ]: 2224451 : assert (r1 >= 0 && r2 >= 0 && r1 < rows && r2 < rows);
[ + - ][ - + ]
[ - + ]
175 [ + + ]: 2682452 : nr_type_t * s = new nr_type_t[cols];
176 : 2224451 : int len = sizeof (nr_type_t) * cols;
177 : 2224451 : memcpy (s, &data[r1 * cols], len);
178 : 2224451 : memcpy (&data[r1 * cols], &data[r2 * cols], len);
179 : 2224451 : memcpy (&data[r2 * cols], s, len);
180 [ + - ][ + - ]: 2224451 : delete[] s;
181 : 2224451 : }
182 : :
183 : : // The function swaps the given columns with each other.
184 : : template <class nr_type_t>
185 : 0 : void tmatrix<nr_type_t>::exchangeCols (int c1, int c2) {
186 [ # # ][ # # ]: 0 : assert (c1 >= 0 && c2 >= 0 && c1 < cols && c2 < cols);
[ # # ][ # # ]
[ # # ]
187 : 0 : nr_type_t s;
188 [ # # ]: 0 : for (int r = 0; r < rows * cols; r += cols) {
189 : 0 : s = data[r + c1];
190 : 0 : data[r + c1] = data[r + c2];
191 : 0 : data[r + c2] = s;
192 : : }
193 : 0 : }
194 : :
195 : : // Compute inverse matrix of the given matrix by Gauss-Jordan elimination.
196 : : template <class nr_type_t>
197 : 0 : tmatrix<nr_type_t> inverse (tmatrix<nr_type_t> a) {
198 : : nr_double_t MaxPivot;
199 : 0 : nr_type_t f;
200 : 0 : tmatrix<nr_type_t> b;
201 : 0 : tmatrix<nr_type_t> e;
202 : 0 : int i, c, r, pivot, n = a.getCols ();
203 : :
204 : : // create temporary matrix and the result matrix
205 [ # # ][ # # ]: 0 : b = tmatrix<nr_type_t> (a);
206 [ # # ][ # # ]: 0 : e = teye<nr_type_t> (n);
207 : :
208 : : // create the eye matrix in 'b' and the result in 'e'
209 [ # # ]: 0 : for (i = 0; i < n; i++) {
210 : : // find maximum column value for pivoting
211 [ # # ]: 0 : for (MaxPivot = 0, pivot = r = i; r < n; r++) {
212 [ # # ]: 0 : if (abs (b.get (r, i)) > MaxPivot) {
[ # # # ]
[ # # ]
213 [ # # ][ # # ]: 0 : MaxPivot = abs (b.get (r, i));
214 : 0 : pivot = r;
215 : : }
216 : : }
217 : : // exchange rows if necessary
218 [ # # ]: 0 : assert (MaxPivot != 0); // singular matrix
219 [ # # ]: 0 : if (i != pivot) {
220 [ # # ]: 0 : b.exchangeRows (i, pivot);
221 [ # # ]: 0 : e.exchangeRows (i, pivot);
222 : : }
223 : :
224 : : // compute current row
225 [ # # ]: 0 : f = b.get (i, i);
226 [ # # ]: 0 : for (c = 0; c < n; c++) {
227 [ # # ][ # # ]: 0 : b.set (i, c, b.get (i, c) / f);
[ # # ]
228 [ # # ][ # # ]: 0 : e.set (i, c, e.get (i, c) / f);
[ # # ]
229 : : }
230 : :
231 : : // compute new rows and columns
232 [ # # ]: 0 : for (r = 0; r < n; r++) {
233 [ # # ]: 0 : if (r != i) {
234 [ # # ]: 0 : f = b.get (r, i);
235 [ # # ]: 0 : for (c = 0; c < n; c++) {
236 [ # # ][ # # ]: 0 : b.set (r, c, b.get (r, c) - f * b.get (i, c));
[ # # ][ # # ]
237 [ # # ][ # # ]: 0 : e.set (r, c, e.get (r, c) - f * e.get (i, c));
[ # # ][ # # ]
238 : : }
239 : : }
240 : : }
241 : : }
242 : 0 : return e;
243 : : }
244 : :
245 : : // Create identity matrix with specified number of rows and columns.
246 : : template <class nr_type_t>
247 : 0 : tmatrix<nr_type_t> teye (int n) {
248 : 0 : tmatrix<nr_type_t> res (n);
249 [ # # ][ # # ]: 0 : for (int r = 0; r < n; r++) res.set (r, r, 1);
[ # ][ # ]
250 : 0 : return res;
251 : : }
252 : :
253 : : // Intrinsic matrix addition.
254 : : template <class nr_type_t>
255 : 10 : tmatrix<nr_type_t> tmatrix<nr_type_t>::operator += (tmatrix<nr_type_t> a) {
256 [ + - ][ - + ]: 10 : assert (a.getRows () == rows && a.getCols () == cols);
[ - + ]
257 : 10 : nr_type_t * src = a.getData ();
258 : 10 : nr_type_t * dst = data;
259 [ + + ]: 2570 : for (int i = 0; i < rows * cols; i++) *dst++ += *src++;
260 : 10 : return *this;
261 : : }
262 : :
263 : : // Intrinsic matrix substraction.
264 : : template <class nr_type_t>
265 : : tmatrix<nr_type_t> tmatrix<nr_type_t>::operator -= (tmatrix<nr_type_t> a) {
266 : : assert (a.getRows () == rows && a.getCols () == cols);
267 : : nr_type_t * src = a.getData ();
268 : : nr_type_t * dst = data;
269 : : for (int i = 0; i < rows * cols; i++) *dst++ -= *src++;
270 : : return *this;
271 : : }
272 : :
273 : : // Matrix multiplication.
274 : : template <class nr_type_t>
275 : : tmatrix<nr_type_t> operator * (tmatrix<nr_type_t> a, tmatrix<nr_type_t> b) {
276 : : assert (a.getCols () == b.getRows ());
277 : : int r, c, i, n = a.getCols ();
278 : : nr_type_t z;
279 : : tmatrix<nr_type_t> res (a.getRows (), b.getCols ());
280 : : for (r = 0; r < a.getRows (); r++) {
281 : : for (c = 0; c < b.getCols (); c++) {
282 : : for (i = 0, z = 0; i < n; i++) z += a.get (r, i) * b.get (i, c);
283 : : res.set (r, c, z);
284 : : }
285 : : }
286 : : return res;
287 : : }
288 : :
289 : : // Multiplication of matrix and vector.
290 : : template <class nr_type_t>
291 : 0 : tvector<nr_type_t> operator * (tmatrix<nr_type_t> a, tvector<nr_type_t> b) {
292 [ # # ]: 0 : assert (a.getCols () == b.getSize ());
293 : 0 : int r, c, n = a.getCols ();
294 : 0 : nr_type_t z;
295 [ # # ]: 0 : tvector<nr_type_t> res (n);
296 : :
297 [ # # ][ # # ]: 0 : for (r = 0; r < n; r++) {
298 [ # # ][ # # ]: 0 : for (c = 0, z = 0; c < n; c++) z += a.get (r, c) * b.get (c);
[ # # ][ # # ]
[ # ][ # ]
[ # # ]
299 [ # # ]: 0 : res.set (r, z);
300 : : }
301 : 0 : return res;
302 : : }
303 : :
304 : : // Multiplication of vector (transposed) and matrix.
305 : : template <class nr_type_t>
306 : 25662 : tvector<nr_type_t> operator * (tvector<nr_type_t> a, tmatrix<nr_type_t> b) {
307 [ - + ]: 25662 : assert (a.getSize () == b.getRows ());
308 : 25662 : int r, c, n = b.getRows ();
309 : 25662 : nr_type_t z;
310 [ + - ]: 25662 : tvector<nr_type_t> res (n);
311 : :
312 [ + + ]: 251160 : for (c = 0; c < n; c++) {
313 [ + - ][ + - ]: 2233140 : for (r = 0, z = 0; r < n; r++) z += a.get (r) * b.get (r, c);
[ + - ][ + + ]
314 [ + - ]: 225498 : res.set (c, z);
315 : : }
316 : 25662 : return res;
317 : : }
318 : :
319 : : // Transpose the matrix in place.
320 : : template <class nr_type_t>
321 : 3003 : void tmatrix<nr_type_t>::transpose (void) {
322 : 3003 : nr_type_t v;
323 [ + + ]: 28665 : for (int r = 0; r < getRows (); r++)
324 [ + + ][ # # ]: 125580 : for (int c = 0; c < r; c++) {
325 [ + - ]: 99918 : v = get (r, c);
326 [ + - ][ + - ]: 99918 : set (r, c, get (c, r));
327 [ + - ]: 99918 : set (c, r, v);
328 : : }
329 : 3003 : }
330 : :
331 : : // Checks validity of matrix.
332 : : template <class nr_type_t>
333 : 214420 : int tmatrix<nr_type_t>::isFinite (void) {
334 [ + + ]: 20623147 : for (int i = 0; i < rows * cols; i++)
335 [ - + ]: 20408727 : if (!std::isfinite (real (data[i]))) return 0;
336 : 214420 : return 1;
337 : : }
338 : :
339 : : #ifdef DEBUG
340 : : // Debug function: Prints the matrix object.
341 : : template <class nr_type_t>
342 : : void tmatrix<nr_type_t>::print (bool realonly) {
343 : : for (int r = 0; r < rows; r++) {
344 : : for (int c = 0; c < cols; c++) {
345 : : if (realonly) {
346 : : fprintf (stderr, "%+.2e%s", (double) real (get (r, c)),
347 : : c != cols - 1 ? " " : "");
348 : : } else {
349 : : fprintf (stderr, "%+.2e%+.2ei%s", (double) real (get (r, c)),
350 : : (double) imag (get (r, c)), c != cols - 1 ? " " : "");
351 : : }
352 : : }
353 : : fprintf (stderr, ";\n");
354 : : }
355 : : }
356 : : #endif /* DEBUG */
357 : :
358 : : } // namespace qucs
|
