My Project  0.0.16
QUCS Mapping
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
tmatrix.cpp
Go to the documentation of this file.
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: tmatrix.cpp 1825 2011-03-11 20:42:14Z ela $
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 <math.h>
34 
35 #include "compat.h"
36 #include "logging.h"
37 #include "complex.h"
38 #include "tmatrix.h"
39 
40 // Constructor creates an unnamed instance of the tmatrix class.
41 template <class nr_type_t>
43  rows = 0;
44  cols = 0;
45  data = NULL;
46 }
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>
52  rows = cols = s;
53  if (s > 0) {
54  data = new nr_type_t[s * s];
55  memset (data, 0, sizeof (nr_type_t) * s * s);
56  }
57  else data = NULL;
58 }
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>
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>
77  rows = m.rows;
78  cols = m.cols;
79  data = NULL;
80 
81  // copy tmatrix elements
82  if (rows > 0 && cols > 0) {
83  data = new nr_type_t[rows * cols];
84  memcpy (data, m.data, sizeof (nr_type_t) * rows * cols);
85  }
86 }
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>&
93  if (&m != this) {
94  rows = m.rows;
95  cols = m.cols;
96  if (data) { delete[] data; data = NULL; }
97  if (rows > 0 && cols > 0) {
98  data = new nr_type_t[rows * cols];
99  memcpy (data, m.data, sizeof (nr_type_t) * rows * cols);
100  }
101  }
102  return *this;
103 }
104 
105 // Destructor deletes a tmatrix object.
106 template <class nr_type_t>
108  if (data) delete[] data;
109 }
110 
111 // Returns the tmatrix element at the given row and column.
112 template <class nr_type_t>
113 nr_type_t tmatrix<nr_type_t>::get (int r, int c) {
114  assert (r >= 0 && r < rows && c >= 0 && c < cols);
115  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 void tmatrix<nr_type_t>::set (int r, int c, nr_type_t z) {
121  assert (r >= 0 && r < rows && c >= 0 && c < cols);
122  data[r * cols + c] = z;
123 }
124 
125 // Sets all the tmatrix elements to the given value.
126 template <class nr_type_t>
127 void tmatrix<nr_type_t>::set (nr_type_t z) {
128  for (int i = 0; i < rows * cols; i++) data[i] = z;
129 }
130 
131 // The function returns the given row in a tvector.
132 template <class nr_type_t>
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>
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>
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>
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 void tmatrix<nr_type_t>::exchangeRows (int r1, int r2) {
174  assert (r1 >= 0 && r2 >= 0 && r1 < rows && r2 < rows);
175  nr_type_t * s = new nr_type_t[cols];
176  int len = sizeof (nr_type_t) * cols;
177  memcpy (s, &data[r1 * cols], len);
178  memcpy (&data[r1 * cols], &data[r2 * cols], len);
179  memcpy (&data[r2 * cols], s, len);
180  delete[] s;
181 }
182 
183 // The function swaps the given columns with each other.
184 template <class nr_type_t>
185 void tmatrix<nr_type_t>::exchangeCols (int c1, int c2) {
186  assert (c1 >= 0 && c2 >= 0 && c1 < cols && c2 < cols);
187  nr_type_t s;
188  for (int r = 0; r < rows * cols; r += cols) {
189  s = data[r + c1];
190  data[r + c1] = data[r + c2];
191  data[r + c2] = s;
192  }
193 }
194 
195 // Compute inverse matrix of the given matrix by Gauss-Jordan elimination.
196 template <class nr_type_t>
198  nr_double_t MaxPivot;
199  nr_type_t f;
202  int i, c, r, pivot, n = a.getCols ();
203 
204  // create temporary matrix and the result matrix
205  b = tmatrix<nr_type_t> (a);
206  e = teye<nr_type_t> (n);
207 
208  // create the eye matrix in 'b' and the result in 'e'
209  for (i = 0; i < n; i++) {
210  // find maximum column value for pivoting
211  for (MaxPivot = 0, pivot = r = i; r < n; r++) {
212  if (abs (b.get (r, i)) > MaxPivot) {
213  MaxPivot = abs (b.get (r, i));
214  pivot = r;
215  }
216  }
217  // exchange rows if necessary
218  assert (MaxPivot != 0); // singular matrix
219  if (i != pivot) {
220  b.exchangeRows (i, pivot);
221  e.exchangeRows (i, pivot);
222  }
223 
224  // compute current row
225  f = b.get (i, i);
226  for (c = 0; c < n; c++) {
227  b.set (i, c, b.get (i, c) / f);
228  e.set (i, c, e.get (i, c) / f);
229  }
230 
231  // compute new rows and columns
232  for (r = 0; r < n; r++) {
233  if (r != i) {
234  f = b.get (r, i);
235  for (c = 0; c < n; c++) {
236  b.set (r, c, b.get (r, c) - f * b.get (i, c));
237  e.set (r, c, e.get (r, c) - f * e.get (i, c));
238  }
239  }
240  }
241  }
242  return e;
243 }
244 
245 // Create identity matrix with specified number of rows and columns.
246 template <class nr_type_t>
248  tmatrix<nr_type_t> res (n);
249  for (int r = 0; r < n; r++) res.set (r, r, 1);
250  return res;
251 }
252 
253 // Intrinsic matrix addition.
254 template <class nr_type_t>
256  assert (a.getRows () == rows && a.getCols () == cols);
257  nr_type_t * src = a.getData ();
258  nr_type_t * dst = data;
259  for (int i = 0; i < rows * cols; i++) *dst++ += *src++;
260  return *this;
261 }
262 
263 // Intrinsic matrix substraction.
264 template <class nr_type_t>
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>
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>
292  assert (a.getCols () == b.getSize ());
293  int r, c, n = a.getCols ();
294  nr_type_t z;
295  tvector<nr_type_t> res (n);
296 
297  for (r = 0; r < n; r++) {
298  for (c = 0, z = 0; c < n; c++) z += a.get (r, c) * b.get (c);
299  res.set (r, z);
300  }
301  return res;
302 }
303 
304 // Multiplication of vector (transposed) and matrix.
305 template <class nr_type_t>
307  assert (a.getSize () == b.getRows ());
308  int r, c, n = b.getRows ();
309  nr_type_t z;
310  tvector<nr_type_t> res (n);
311 
312  for (c = 0; c < n; c++) {
313  for (r = 0, z = 0; r < n; r++) z += a.get (r) * b.get (r, c);
314  res.set (c, z);
315  }
316  return res;
317 }
318 
319 // Transpose the matrix in place.
320 template <class nr_type_t>
322  nr_type_t v;
323  for (int r = 0; r < getRows (); r++)
324  for (int c = 0; c < r; c++) {
325  v = get (r, c);
326  set (r, c, get (c, r));
327  set (c, r, v);
328  }
329 }
330 
331 // Checks validity of matrix.
332 template <class nr_type_t>
334  for (int i = 0; i < rows * cols; i++)
335  if (!finite (real (data[i]))) return 0;
336  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 */