Qucs-core
0.0.19
|
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