Qucs-core  0.0.19
tmatrix.cpp
Go to the documentation of this file.
00001 /*
00002  * tmatrix.cpp - simple matrix template class implementation
00003  *
00004  * Copyright (C) 2004, 2005, 2006, 2008 Stefan Jahn <stefan@lkcc.org>
00005  *
00006  * This is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2, or (at your option)
00009  * any later version.
00010  *
00011  * This software is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this package; see the file COPYING.  If not, write to
00018  * the Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
00019  * Boston, MA 02110-1301, USA.
00020  *
00021  * $Id$
00022  *
00023  */
00024 
00025 #include "qucs_typedefs.h"
00026 
00027 #include <assert.h>
00028 #include <stdio.h>
00029 #include <stdlib.h>
00030 #include <string.h>
00031 #include <cmath>
00032 
00033 #include "compat.h"
00034 #include "logging.h"
00035 #include "complex.h"
00036 #include "tmatrix.h"
00037 
00038 namespace qucs {
00039 
00040 // Constructor creates an unnamed instance of the tmatrix class.
00041 template <class nr_type_t>
00042 tmatrix<nr_type_t>::tmatrix () {
00043   rows = 0;
00044   cols = 0;
00045   data = NULL;
00046 }
00047 
00048 /* Constructor creates an unnamed instance of the tmatrix class with a
00049    certain number of rows and columns.  Creates a square tmatrix.  */
00050 template <class nr_type_t>
00051 tmatrix<nr_type_t>::tmatrix (int s)  {
00052   rows = cols = s;
00053   if (s > 0) {
00054     data = new nr_type_t[s * s];
00055     memset (data, 0, sizeof (nr_type_t) * s * s);
00056   }
00057   else data = NULL;
00058 }
00059 
00060 /* Constructor creates an unnamed instance of the tmatrix class with a
00061    certain number of rows and columns.  */
00062 template <class nr_type_t>
00063 tmatrix<nr_type_t>::tmatrix (int r, int c)  {
00064   rows = r;
00065   cols = c;
00066   if (r > 0 && c > 0) {
00067     data = new nr_type_t[r * c];
00068     memset (data, 0, sizeof (nr_type_t) * r * c);
00069   }
00070   else data = NULL;
00071 }
00072 
00073 /* The copy constructor creates a new instance based on the given
00074    tmatrix object. */
00075 template <class nr_type_t>
00076 tmatrix<nr_type_t>::tmatrix (const tmatrix & m) {
00077   rows = m.rows;
00078   cols = m.cols;
00079   data = NULL;
00080 
00081   // copy tmatrix elements
00082   if (rows > 0 && cols > 0) {
00083     data = new nr_type_t[rows * cols];
00084     memcpy (data, m.data, sizeof (nr_type_t) * rows * cols);
00085   }
00086 }
00087 
00088 /* The assignment copy constructor creates a new instance based on the
00089    given tmatrix object. */
00090 template <class nr_type_t>
00091 const tmatrix<nr_type_t>&
00092 tmatrix<nr_type_t>::operator=(const tmatrix<nr_type_t> & m) {
00093   if (&m != this) {
00094     rows = m.rows;
00095     cols = m.cols;
00096     delete[] data;
00097     if (rows > 0 && cols > 0) {
00098       data = new nr_type_t[rows * cols];
00099       memcpy (data, m.data, sizeof (nr_type_t) * rows * cols);
00100     }
00101   }
00102   return *this;
00103 }
00104 
00105 // Destructor deletes a tmatrix object.
00106 template <class nr_type_t>
00107 tmatrix<nr_type_t>::~tmatrix () {
00108   delete[] data;
00109 }
00110 
00111 // Returns the tmatrix element at the given row and column.
00112 template <class nr_type_t>
00113 nr_type_t tmatrix<nr_type_t>::get (int r, int c) {
00114   assert (r >= 0 && r < rows && c >= 0 && c < cols);
00115   return data[r * cols + c];
00116 }
00117 
00118 // Sets the tmatrix element at the given row and column.
00119 template <class nr_type_t>
00120 void tmatrix<nr_type_t>::set (int r, int c, nr_type_t z) {
00121   assert (r >= 0 && r < rows && c >= 0 && c < cols);
00122   data[r * cols + c] = z;
00123 }
00124 
00125 // Sets all the tmatrix elements to the given value.
00126 template <class nr_type_t>
00127 void tmatrix<nr_type_t>::set (nr_type_t z) {
00128   for (int i = 0; i < rows * cols; i++) data[i] = z;
00129 }
00130 
00131 // The function returns the given row in a tvector.
00132 template <class nr_type_t>
00133 tvector<nr_type_t> tmatrix<nr_type_t>::getRow (int r) {
00134   assert (r >= 0 && r < rows);
00135   tvector<nr_type_t> res (cols);
00136   nr_type_t * dst = res.getData ();
00137   nr_type_t * src = &data[r * cols];
00138   memcpy (dst, src, sizeof (nr_type_t) * cols);
00139   return res;
00140 }
00141 
00142 // Puts the given tvector into the given row of the tmatrix instance.
00143 template <class nr_type_t>
00144 void tmatrix<nr_type_t>::setRow (int r, tvector<nr_type_t> v) {
00145   assert (r >= 0 && r < rows && v.getSize () == cols);
00146   nr_type_t * dst = &data[r * cols];
00147   nr_type_t * src = v.getData ();
00148   memcpy (dst, src, sizeof (nr_type_t) * cols);
00149 }
00150 
00151 // The function returns the given column in a tvector.
00152 template <class nr_type_t>
00153 tvector<nr_type_t> tmatrix<nr_type_t>::getCol (int c) {
00154   assert (c >= 0 && c < cols);
00155   tvector<nr_type_t> res (rows);
00156   nr_type_t * dst = res.getData ();
00157   nr_type_t * src = &data[c];
00158   for (int r = 0; r < rows; r++, src += cols, dst++) *dst = *src;
00159   return res;
00160 }
00161 
00162 // Puts the given tvector into the given column of the tmatrix instance.
00163 template <class nr_type_t>
00164 void tmatrix<nr_type_t>::setCol (int c, tvector<nr_type_t> v) {
00165   assert (c >= 0 && c < cols && v.getSize () == rows);
00166   nr_type_t * dst = &data[c];
00167   nr_type_t * src = v.getData ();
00168   for (int r = 0; r < rows; r++, src++, dst += cols) *dst = *src;
00169 }
00170 
00171 // The function swaps the given rows with each other.
00172 template <class nr_type_t>
00173 void tmatrix<nr_type_t>::exchangeRows (int r1, int r2) {
00174   assert (r1 >= 0 && r2 >= 0 && r1 < rows && r2 < rows);
00175   nr_type_t * s = new nr_type_t[cols];
00176   int len = sizeof (nr_type_t) * cols;
00177   memcpy (s, &data[r1 * cols], len);
00178   memcpy (&data[r1 * cols], &data[r2 * cols], len);
00179   memcpy (&data[r2 * cols], s, len);
00180   delete[] s;
00181 }
00182 
00183 // The function swaps the given columns with each other.
00184 template <class nr_type_t>
00185 void tmatrix<nr_type_t>::exchangeCols (int c1, int c2) {
00186   assert (c1 >= 0 && c2 >= 0 && c1 < cols && c2 < cols);
00187   nr_type_t s;
00188   for (int r = 0; r < rows * cols; r += cols) {
00189     s = data[r + c1];
00190     data[r + c1] = data[r + c2];
00191     data[r + c2] = s;
00192   }
00193 }
00194 
00195 // Compute inverse matrix of the given matrix by Gauss-Jordan elimination.
00196 template <class nr_type_t>
00197 tmatrix<nr_type_t> inverse (tmatrix<nr_type_t> a) {
00198   nr_double_t MaxPivot;
00199   nr_type_t f;
00200   tmatrix<nr_type_t> b;
00201   tmatrix<nr_type_t> e;
00202   int i, c, r, pivot, n = a.getCols ();
00203 
00204   // create temporary matrix and the result matrix
00205   b = tmatrix<nr_type_t> (a);
00206   e = teye<nr_type_t> (n);
00207 
00208   // create the eye matrix in 'b' and the result in 'e'
00209   for (i = 0; i < n; i++) {
00210     // find maximum column value for pivoting
00211     for (MaxPivot = 0, pivot = r = i; r < n; r++) {
00212       if (abs (b.get (r, i)) > MaxPivot) {
00213         MaxPivot = abs (b.get (r, i));
00214         pivot = r;
00215       }
00216     }
00217     // exchange rows if necessary
00218     assert (MaxPivot != 0); // singular matrix
00219     if (i != pivot) {
00220       b.exchangeRows (i, pivot);
00221       e.exchangeRows (i, pivot);
00222     }
00223 
00224     // compute current row
00225     f = b.get (i, i);
00226     for (c = 0; c < n; c++) {
00227       b.set (i, c, b.get (i, c) / f);
00228       e.set (i, c, e.get (i, c) / f);
00229     }
00230 
00231     // compute new rows and columns
00232     for (r = 0; r < n; r++) {
00233       if (r != i) {
00234         f = b.get (r, i);
00235         for (c = 0; c < n; c++) {
00236           b.set (r, c, b.get (r, c) - f * b.get (i, c));
00237           e.set (r, c, e.get (r, c) - f * e.get (i, c));
00238         }
00239       }
00240     }
00241   }
00242   return e;
00243 }
00244 
00245 // Create identity matrix with specified number of rows and columns.
00246 template <class nr_type_t>
00247 tmatrix<nr_type_t> teye (int n) {
00248   tmatrix<nr_type_t> res (n);
00249   for (int r = 0; r < n; r++) res.set (r, r, 1);
00250   return res;
00251 }
00252 
00253 // Intrinsic matrix addition.
00254 template <class nr_type_t>
00255 tmatrix<nr_type_t> tmatrix<nr_type_t>::operator += (tmatrix<nr_type_t> a) {
00256   assert (a.getRows () == rows && a.getCols () == cols);
00257   nr_type_t * src = a.getData ();
00258   nr_type_t * dst = data;
00259   for (int i = 0; i < rows * cols; i++) *dst++ += *src++;
00260   return *this;
00261 }
00262 
00263 // Intrinsic matrix substraction.
00264 template <class nr_type_t>
00265 tmatrix<nr_type_t> tmatrix<nr_type_t>::operator -= (tmatrix<nr_type_t> a) {
00266   assert (a.getRows () == rows && a.getCols () == cols);
00267   nr_type_t * src = a.getData ();
00268   nr_type_t * dst = data;
00269   for (int i = 0; i < rows * cols; i++) *dst++ -= *src++;
00270   return *this;
00271 }
00272 
00273 // Matrix multiplication.
00274 template <class nr_type_t>
00275 tmatrix<nr_type_t> operator * (tmatrix<nr_type_t> a, tmatrix<nr_type_t> b) {
00276   assert (a.getCols () == b.getRows ());
00277   int r, c, i, n = a.getCols ();
00278   nr_type_t z;
00279   tmatrix<nr_type_t> res (a.getRows (), b.getCols ());
00280   for (r = 0; r < a.getRows (); r++) {
00281     for (c = 0; c < b.getCols (); c++) {
00282       for (i = 0, z = 0; i < n; i++) z += a.get (r, i) * b.get (i, c);
00283       res.set (r, c, z);
00284     }
00285   }
00286   return res;
00287 }
00288 
00289 // Multiplication of matrix and vector.
00290 template <class nr_type_t>
00291 tvector<nr_type_t> operator * (tmatrix<nr_type_t> a, tvector<nr_type_t> b) {
00292   assert (a.getCols () == b.size ());
00293   int r, c, n = a.getCols ();
00294   nr_type_t z;
00295   tvector<nr_type_t> res (n);
00296 
00297   for (r = 0; r < n; r++) {
00298     for (c = 0, z = 0; c < n; c++) z += a.get (r, c) * b.get (c);
00299     res.set (r, z);
00300   }
00301   return res;
00302 }
00303 
00304 // Multiplication of vector (transposed) and matrix.
00305 template <class nr_type_t>
00306 tvector<nr_type_t> operator * (tvector<nr_type_t> a, tmatrix<nr_type_t> b) {
00307   assert (a.size () == b.getRows ());
00308   int r, c, n = b.getRows ();
00309   nr_type_t z;
00310   tvector<nr_type_t> res (n);
00311 
00312   for (c = 0; c < n; c++) {
00313     for (r = 0, z = 0; r < n; r++) z += a.get (r) * b.get (r, c);
00314     res.set (c, z);
00315   }
00316   return res;
00317 }
00318 
00319 // Transpose the matrix in place.
00320 template <class nr_type_t>
00321 void tmatrix<nr_type_t>::transpose (void) {
00322   nr_type_t v;
00323   for (int r = 0; r < getRows (); r++)
00324     for (int c = 0; c < r; c++) {
00325       v = get (r, c);
00326       set (r, c, get (c, r));
00327       set (c, r, v);
00328     }
00329 }
00330 
00331 // Checks validity of matrix.
00332 template <class nr_type_t>
00333 int tmatrix<nr_type_t>::isFinite (void) {
00334   for (int i = 0; i < rows * cols; i++)
00335     if (!std::isfinite (real (data[i]))) return 0;
00336   return 1;
00337 }
00338 
00339 #ifdef DEBUG
00340 // Debug function: Prints the matrix object.
00341 template <class nr_type_t>
00342 void tmatrix<nr_type_t>::print (bool realonly) {
00343   for (int r = 0; r < rows; r++) {
00344     for (int c = 0; c < cols; c++) {
00345       if (realonly) {
00346         fprintf (stderr, "%+.2e%s", (double) real (get (r, c)),
00347                  c != cols - 1 ? " " : "");
00348       } else {
00349         fprintf (stderr, "%+.2e%+.2ei%s", (double) real (get (r, c)),
00350                  (double) imag (get (r, c)), c != cols - 1 ? " " : "");
00351       }
00352     }
00353     fprintf (stderr, ";\n");
00354   }
00355 }
00356 #endif /* DEBUG */
00357 
00358 } // namespace qucs