1 /******************************************************************************\
2 * Copyright (c) 2001
3 *
4 * Author(s):
5 * Volker Fischer
6 *
7 * Description:
8 * c++ Mathematic Library (Matlib)
9 *
10 ******************************************************************************
11 *
12 * This program is free software; you can redistribute it and/or modify it under
13 * the terms of the GNU General Public License as published by the Free Software
14 * Foundation; either version 2 of the License, or (at your option) any later
15 * version.
16 *
17 * This program is distributed in the hope that it will be useful, but WITHOUT
18 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
19 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
20 * details.
21 *
22 * You should have received a copy of the GNU General Public License along with
23 * this program; if not, write to the Free Software Foundation, Inc.,
24 * 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
25 *
26 \******************************************************************************/
27
28 #ifndef _MATLIB_H_
29 #define _MATLIB_H_
30
31 #include <cstdlib>
32 #include <math.h>
33 #include <complex>
34 using namespace std;
35 #include "../GlobalDefinitions.h"
36
37
38 /* Definitions ****************************************************************/
39 /* Two different types: constant and temporary buffer */
40 enum EVecTy {VTY_CONST, VTY_TEMP};
41
42
43 /* These definitions save a lot of redundant code */
44 #define _VECOP(TYPE, LENGTH, FCT) const int iL = LENGTH; \
45 CMatlibVector<TYPE> vecRet(iL, VTY_TEMP); \
46 for (int i = 0; i < iL; i++) \
47 vecRet[i] = FCT; \
48 return vecRet
49
50 #define _VECOPCL(FCT) for (int i = 0; i < iVectorLength; i++) \
51 operator[](i) FCT; \
52 return *this
53
54 #define _MATOP(TYPE, LENGTHROW, LENGTHCOL, FCT) \
55 const int iRL = LENGTHROW; \
56 CMatlibMatrix<TYPE> matRet(iRL, LENGTHCOL, VTY_TEMP); \
57 for (int i = 0; i < iRL; i++) \
58 matRet[i] = FCT; \
59 return matRet
60
61 #define _MATOPCL(FCT) for (int i = 0; i < iRowSize; i++) \
62 operator[](i) FCT; \
63 return *this
64
65
66 /* In debug mode, test input parameters */
67 #ifdef _DEBUG_
68 #define _TESTRNGR(POS) if ((POS >= iVectorLength) || (POS < 0)) \
69 DebugError("MatLibrRead", "POS", POS, \
70 "iVectorLength", iVectorLength)
71 #define _TESTRNGW(POS) if ((POS >= iVectorLength) || (POS < 0)) \
72 DebugError("MatLibrWrite", "POS", POS, \
73 "iVectorLength", iVectorLength)
74 #define _TESTSIZE(INP) if (INP != iVectorLength) \
75 DebugError("SizeCheck", "INP", INP, \
76 "iVectorLength", iVectorLength)
77 #define _TESTRNGRM(POS) if ((POS >= iRowSize) || (POS < 0)) \
78 DebugError("MatLibrReadMatrix", "POS", POS, \
79 "iRowSize", iRowSize)
80 #define _TESTRNGWM(POS) if ((POS >= iRowSize) || (POS < 0)) \
81 DebugError("MatLibrWriteMatrix", "POS", POS, \
82 "iRowSize", iRowSize)
83 #define _TESTSIZEM(INP) if (INP != iRowSize) \
84 DebugError("MatLibOperatorMatrix=()", "INP", INP, \
85 "iRowSize", iRowSize)
86 #else
87
88 // On Visual c++ 2005 Express Edition there is a segmentation fault if these macros are empty
89 // TODO: FIX this with a better solution
90 #ifdef _MSC_VER
91 #define _TESTRNGR(POS) if (POS != POS) int idummy=0
92 #define _TESTRNGW(POS) if (POS != POS) int idummy=0
93 #define _TESTSIZE(INP) if (INP != INP) int idummy=0
94 #define _TESTRNGRM(POS) if (POS != POS) int idummy=0
95 #define _TESTRNGWM(POS) if (POS != POS) int idummy=0
96 #define _TESTSIZEM(INP) if (INP != INP) int idummy=0
97 #else
98 #define _TESTRNGR(POS)
99 #define _TESTRNGW(POS)
100 #define _TESTSIZE(INP)
101 #define _TESTRNGRM(POS)
102 #define _TESTRNGWM(POS)
103 #define _TESTSIZEM(INP)
104 #endif
105 #endif
106
107
108 /* Classes ********************************************************************/
109 /* Prototypes */
110 template<class T> class CMatlibVector;
111 template<class T> class CMatlibMatrix;
112
113 /* Here we can choose the precision of the Matlib calculations */
114 typedef _REAL CReal;
115 typedef complex<CReal> CComplex;
116 typedef CMatlibVector<CReal> CRealVector;
117 typedef CMatlibVector<CComplex> CComplexVector;
118 typedef CMatlibMatrix<CReal> CRealMatrix;
119 typedef CMatlibMatrix<CComplex> CComplexMatrix;
120
121
122 /******************************************************************************/
123 /* CMatlibVector class ********************************************************/
124 /******************************************************************************/
125 template<class T>
126 class CMatlibVector
127 {
128 public:
129 /* Construction, Destruction -------------------------------------------- */
CMatlibVector()130 CMatlibVector() : eVType(VTY_CONST), iVectorLength(0), pData(NULL) {}
131 CMatlibVector(const int iNLen, const EVecTy eNTy = VTY_CONST) :
eVType(eNTy)132 eVType(eNTy), iVectorLength(0), pData(NULL) {Init(iNLen);}
CMatlibVector(const int iNLen,const T tIniVal)133 CMatlibVector(const int iNLen, const T tIniVal) :
134 eVType(VTY_CONST), iVectorLength(0), pData(NULL) {Init(iNLen, tIniVal);}
135 CMatlibVector(CMatlibVector<T>& vecI);
136 CMatlibVector(const CMatlibVector<T>& vecI);
~CMatlibVector()137 virtual ~CMatlibVector() {if (pData != NULL) delete[] pData;}
138
CMatlibVector(const CMatlibVector<CReal> & fvReal,const CMatlibVector<CReal> & fvImag)139 CMatlibVector(const CMatlibVector<CReal>& fvReal, const CMatlibVector<CReal>& fvImag) :
140 eVType(VTY_CONST/*VTY_TEMP*/), iVectorLength(fvReal.GetSize()), pData(NULL)
141 {
142 /* Allocate data block for vector */
143 pData = new CComplex[iVectorLength];
144
145 /* Copy data from real-vectors in complex vector */
146 for (int i = 0; i < iVectorLength; i++)
147 pData[i] = CComplex(fvReal[i], fvImag[i]);
148 }
149
150 /* Operator[] (Regular indices!!!) */
151 inline T& operator[](int const iPos) const
152 {_TESTRNGR(iPos); return pData[iPos];}
153 inline T& operator[](int const iPos)
154 {_TESTRNGW(iPos); return pData[iPos];} // For use as l value
155
156 /* Operator() */
operator()157 inline T& operator()(int const iPos) const
158 {_TESTRNGR(iPos - 1); return pData[iPos - 1];}
operator()159 inline T& operator()(int const iPos)
160 {_TESTRNGW(iPos - 1); return pData[iPos - 1];} // For use as l value
161
162 CMatlibVector<T> operator()(const int iFrom, const int iTo) const;
163 CMatlibVector<T> operator()(const int iFrom, const int iStep, const int iTo) const;
164
GetSize()165 inline int GetSize() const {return iVectorLength;}
166 void Init(const int iIniLen, const T tIniVal = 0);
167 inline CMatlibVector<T>& Reset(const T tResVal = 0) {_VECOPCL(= tResVal);}
168 CMatlibVector<T>& PutIn(const int iFrom, const int iTo, CMatlibVector<T>& fvA);
169 CMatlibVector<T>& Merge(const CMatlibVector<T>& vecA, T& tB);
170 CMatlibVector<T>& Merge(const CMatlibVector<T>& vecA, const CMatlibVector<T>& vecB);
171 CMatlibVector<T>& Merge(const CMatlibVector<T>& vecA, const CMatlibVector<T>& vecB,
172 const CMatlibVector<T>& vecC);
173
174
175 /* operator= */
176 inline CMatlibVector<T>& operator=(const CMatlibVector<CReal>& vecI)
177 {
178 Init(vecI.GetSize());
179 for (int i = 0; i < iVectorLength; i++)
180 operator[](i) = vecI[i];
181
182 return *this;
183 }
184
185 inline CMatlibVector<CComplex>& operator=(const CMatlibVector<CComplex>& vecI)
186 {
187 Init(vecI.GetSize());
188 for (int i = 0; i < iVectorLength; i++)
189 operator[](i) = vecI[i];
190
191 return *this;
192 }
193
194 /* operator*= */
195 inline CMatlibVector<T>& operator*=(const CReal& rI)
196 {_VECOPCL(*= rI);}
197 inline CMatlibVector<CComplex>& operator*=(const CComplex& cI)
198 {_VECOPCL(*= cI);}
199 inline CMatlibVector<T>& operator*=(const CMatlibVector<CReal>& vecI)
200 {_VECOPCL(*= vecI[i]);}
201 inline CMatlibVector<CComplex>& operator*=(const CMatlibVector<CComplex>& vecI)
202 {_VECOPCL(*= vecI[i]);}
203
204 /* operator/= */
205 inline CMatlibVector<T>& operator/=(const CReal& rI)
206 {_VECOPCL(/= rI);}
207 inline CMatlibVector<CComplex>& operator/=(const CComplex& cI)
208 {_VECOPCL(/= cI);}
209 inline CMatlibVector<T>& operator/=(const CMatlibVector<CReal>& vecI)
210 {_VECOPCL(/= vecI[i]);}
211 inline CMatlibVector<CComplex>& operator/=(const CMatlibVector<CComplex>& vecI)
212 {_VECOPCL(/= vecI[i]);}
213
214 /* operator+= */
215 inline CMatlibVector<T>& operator+=(const CReal& rI)
216 {_VECOPCL(+= rI);}
217 inline CMatlibVector<CComplex>& operator+=(const CComplex& cI)
218 {_VECOPCL(+= cI);}
219 inline CMatlibVector<T>& operator+=(const CMatlibVector<CReal>& vecI)
220 {_VECOPCL(+= vecI[i]);}
221 inline CMatlibVector<CComplex>& operator+=(const CMatlibVector<CComplex>& vecI)
222 {_VECOPCL(+= vecI[i]);}
223
224 /* operator-= */
225 inline CMatlibVector<T>& operator-=(const CReal& rI)
226 {_VECOPCL(-= rI);}
227 inline CMatlibVector<CComplex>& operator-=(const CComplex& cI)
228 {_VECOPCL(-= cI);}
229 inline CMatlibVector<T>& operator-=(const CMatlibVector<CReal>& vecI)
230 {_VECOPCL(-= vecI[i]);}
231 inline CMatlibVector<CComplex>& operator-=(const CMatlibVector<CComplex>& vecI)
232 {_VECOPCL(-= vecI[i]);}
233
234 protected:
235 EVecTy eVType;
236 int iVectorLength;
237 T* pData;
238 };
239
240 /* operator* ---------------------------------------------------------------- */
241 inline CMatlibVector<CComplex> // cv, cv
242 operator*(const CMatlibVector<CComplex>& cvA, const CMatlibVector<CComplex>& cvB)
243 {_VECOP(CComplex, cvA.GetSize(), cvA[i] * cvB[i]);}
244 inline CMatlibVector<CReal> // rv, rv
245 operator*(const CMatlibVector<CReal>& rvA, const CMatlibVector<CReal>& rvB)
246 {_VECOP(CReal, rvA.GetSize(), rvA[i] * rvB[i]);}
247
248 inline CMatlibVector<CComplex> // cv, rv
249 operator*(const CMatlibVector<CComplex>& cvA, const CMatlibVector<CReal>& rvB)
250 {_VECOP(CComplex, cvA.GetSize(), cvA[i] * rvB[i]);}
251 inline CMatlibVector<CComplex> // rv, cv
252 operator*(const CMatlibVector<CReal>& rvB, const CMatlibVector<CComplex>& cvA)
253 {_VECOP(CComplex, cvA.GetSize(), cvA[i] * rvB[i]);}
254
255 template<class T> inline
256 CMatlibVector<T> // Tv, r
257 operator*(const CMatlibVector<T>& vecA, const CReal& rB)
258 {_VECOP(T, vecA.GetSize(), vecA[i] * rB);}
259 template<class T> inline
260 CMatlibVector<T> // r, Tv
261 operator*(const CReal& rA, const CMatlibVector<T>& vecB)
262 {_VECOP(T, vecB.GetSize(), rA * vecB[i]);}
263
264 template<class T> inline
265 CMatlibVector<CComplex> // Tv, c
266 operator*(const CMatlibVector<T>& vecA, const CComplex& cB)
267 {_VECOP(CComplex, vecA.GetSize(), vecA[i] * cB);}
268 template<class T> inline
269 CMatlibVector<CComplex> // c, Tv
270 operator*(const CComplex& cA, const CMatlibVector<T>& vecB)
271 {_VECOP(CComplex, vecB.GetSize(), cA * vecB[i]);}
272
273
274 /* operator/ ---------------------------------------------------------------- */
275 inline CMatlibVector<CComplex> // cv, cv
276 operator/(const CMatlibVector<CComplex>& cvA, const CMatlibVector<CComplex>& cvB)
277 {_VECOP(CComplex, cvA.GetSize(), cvA[i] / cvB[i]);}
278 inline CMatlibVector<CReal> // rv, rv
279 operator/(const CMatlibVector<CReal>& rvA, const CMatlibVector<CReal>& rvB)
280 {_VECOP(CReal, rvA.GetSize(), rvA[i] / rvB[i]);}
281
282 inline CMatlibVector<CComplex> // cv, rv
283 operator/(const CMatlibVector<CComplex>& cvA, const CMatlibVector<CReal>& rvB)
284 {_VECOP(CComplex, cvA.GetSize(), cvA[i] / rvB[i]);}
285 inline CMatlibVector<CComplex> // rv, cv
286 operator/(const CMatlibVector<CReal>& rvA, const CMatlibVector<CComplex>& cvB)
287 {_VECOP(CComplex, rvA.GetSize(), rvA[i] / cvB[i]);}
288
289 template<class T> inline
290 CMatlibVector<T> // Tv, r
291 operator/(const CMatlibVector<T>& vecA, const CReal& rB)
292 {_VECOP(T, vecA.GetSize(), vecA[i] / rB);}
293 template<class T> inline
294 CMatlibVector<T> // r, Tv
295 operator/(const CReal& rA, const CMatlibVector<T>& vecB)
296 {_VECOP(T, vecB.GetSize(), rA / vecB[i]);}
297
298 template<class T> inline
299 CMatlibVector<CComplex> // Tv, c
300 operator/(const CMatlibVector<T>& vecA, const CComplex& cB)
301 {_VECOP(CComplex, vecA.GetSize(), vecA[i] / cB);}
302 template<class T> inline
303 CMatlibVector<CComplex> // c, Tv
304 operator/(const CComplex& cA, const CMatlibVector<T>& vecB)
305 {_VECOP(CComplex, vecB.GetSize(), cA / vecB[i]);}
306
307
308 /* operator+ ---------------------------------------------------------------- */
309 inline CMatlibVector<CComplex> // cv, cv
310 operator+(const CMatlibVector<CComplex>& cvA, const CMatlibVector<CComplex>& cvB)
311 {_VECOP(CComplex, cvA.GetSize(), cvA[i] + cvB[i]);}
312 inline CMatlibVector<CReal> // rv, rv
313 operator+(const CMatlibVector<CReal>& rvA, const CMatlibVector<CReal>& rvB)
314 {_VECOP(CReal, rvA.GetSize(), rvA[i] + rvB[i]);}
315
316 inline CMatlibVector<CComplex> // cv, rv
317 operator+(const CMatlibVector<CComplex>& cvA, const CMatlibVector<CReal>& rvB)
318 {_VECOP(CComplex, cvA.GetSize(), cvA[i] + rvB[i]);}
319 inline CMatlibVector<CComplex> // rv, cv
320 operator+(const CMatlibVector<CReal>& rvA, const CMatlibVector<CComplex>& cvB)
321 {_VECOP(CComplex, rvA.GetSize(), rvA[i] + cvB[i]);}
322
323 template<class T> inline
324 CMatlibVector<T> // Tv, r
325 operator+(const CMatlibVector<T>& vecA, const CReal& rB)
326 {_VECOP(T, vecA.GetSize(), vecA[i] + rB);}
327 template<class T> inline
328 CMatlibVector<T> // r, Tv
329 operator+(const CReal& rA, const CMatlibVector<T>& vecB)
330 {_VECOP(T, vecB.GetSize(), rA + vecB[i]);}
331
332 template<class T> inline
333 CMatlibVector<CComplex> // Tv, c
334 operator+(const CMatlibVector<T>& vecA, const CComplex& cB)
335 {_VECOP(CComplex, vecA.GetSize(), vecA[i] + cB);}
336 template<class T> inline
337 CMatlibVector<CComplex> // c, Tv
338 operator+(const CComplex& cA, const CMatlibVector<T>& vecB)
339 {_VECOP(CComplex, vecB.GetSize(), cA + vecB[i]);}
340
341
342 /* operator- ---------------------------------------------------------------- */
343 inline CMatlibVector<CComplex> // cv, cv
344 operator-(const CMatlibVector<CComplex>& cvA, const CMatlibVector<CComplex>& cvB)
345 {_VECOP(CComplex, cvA.GetSize(), cvA[i] - cvB[i]);}
346 inline CMatlibVector<CReal> // rv, rv
347 operator-(const CMatlibVector<CReal>& rvA, const CMatlibVector<CReal>& rvB)
348 {_VECOP(CReal, rvA.GetSize(), rvA[i] - rvB[i]);}
349
350 inline CMatlibVector<CComplex> // cv, rv
351 operator-(const CMatlibVector<CComplex>& cvA, const CMatlibVector<CReal>& rvB)
352 {_VECOP(CComplex, cvA.GetSize(), cvA[i] - rvB[i]);}
353 inline CMatlibVector<CComplex> // rv, cv
354 operator-(const CMatlibVector<CReal>& rvA, const CMatlibVector<CComplex>& cvB)
355 {_VECOP(CComplex, rvA.GetSize(), rvA[i] - cvB[i]);}
356
357 template<class T> inline
358 CMatlibVector<T> // Tv, r
359 operator-(const CMatlibVector<T>& vecA, const CReal& rB)
360 {_VECOP(T, vecA.GetSize(), vecA[i] - rB);}
361 template<class T> inline
362 CMatlibVector<T> // r, Tv
363 operator-(const CReal& rA, const CMatlibVector<T>& vecB)
364 {_VECOP(T, vecB.GetSize(), rA - vecB[i]);}
365
366 template<class T> inline
367 CMatlibVector<CComplex> // Tv, c
368 operator-(const CMatlibVector<T>& vecA, const CComplex& cB)
369 {_VECOP(CComplex, vecA.GetSize(), vecA[i] - cB);}
370 template<class T> inline
371 CMatlibVector<CComplex> // c, Tv
372 operator-(const CComplex& cA, const CMatlibVector<T>& vecB)
373 {_VECOP(CComplex, vecB.GetSize(), cA - vecB[i]);}
374
375
376 /* Implementation **************************************************************
377 (the implementation of template classes must be in the header file!) */
378 template<class T>
CMatlibVector(CMatlibVector<T> & vecI)379 CMatlibVector<T>::CMatlibVector(CMatlibVector<T>& vecI) :
380 eVType(VTY_CONST/*VTY_TEMP*/), iVectorLength(vecI.GetSize()), pData(NULL)
381 {
382 /* The copy constructor for the constant vector is a real copying
383 task. But in the case of a temporary buffer only the pointer
384 of the temporary buffer is used. The buffer of the temporary
385 vector is then destroyed!!! Therefore the usage of "VTY_TEMP"
386 should be done if the vector IS NOT USED IN A FUNCTION CALL,
387 otherwise this vector will be destroyed afterwards (if the
388 function argument is not declared with "&") */
389 if (iVectorLength > 0)
390 {
391 if (vecI.eVType == VTY_CONST)
392 {
393 /* Allocate data block for vector */
394 pData = new T[iVectorLength];
395
396 /* Copy vector */
397 for (int i = 0; i < iVectorLength; i++)
398 pData[i] = vecI[i];
399 }
400 else
401 {
402 /* We can define the copy constructor as a destroying operator of
403 the input vector for performance reasons. This
404 saves us from always copy the entire vector */
405 /* Take data pointer from input vector (steal it) */
406 pData = vecI.pData;
407
408 /* Destroy other vector (temporary vectors only) */
409 vecI.pData = NULL;
410 }
411 }
412 }
413
414 /* Copy constructor for constant Matlib vectors */
415 template<class T>
CMatlibVector(const CMatlibVector<T> & vecI)416 CMatlibVector<T>::CMatlibVector(const CMatlibVector<T>& vecI) :
417 eVType(VTY_CONST), iVectorLength(vecI.GetSize()), pData(NULL)
418 {
419 if (iVectorLength > 0)
420 {
421 /* Allocate data block for vector */
422 pData = new T[iVectorLength];
423
424 /* Copy vector */
425 for (int i = 0; i < iVectorLength; i++)
426 pData[i] = vecI[i];
427 }
428 }
429
430 template<class T>
Init(const int iIniLen,const T tIniVal)431 void CMatlibVector<T>::Init(const int iIniLen, const T tIniVal)
432 {
433 iVectorLength = iIniLen;
434
435 /* Allocate data block for vector */
436 if (iVectorLength > 0)
437 {
438 if (pData != NULL)
439 delete[] pData;
440
441 pData = new T[iVectorLength];
442
443 /* Init with init value */
444 for (int i = 0; i < iVectorLength; i++)
445 pData[i] = tIniVal;
446 }
447 }
448
449 template<class T> inline
operator()450 CMatlibVector<T> CMatlibVector<T>::operator()(const int iFrom,
451 const int iTo) const
452 {
453 /* This is also capable of "wrap around" blocks (if the value in "iFrom" is
454 larger then the "iTo" value) */
455 int i;
456 const int iStartVal = iFrom - 1;
457
458 if (iFrom > iTo)
459 {
460 /* Wrap around case */
461 CMatlibVector<T> vecRet(iVectorLength - iStartVal + iTo, VTY_TEMP);
462
463 int iCurPos = 0;
464
465 for (i = iStartVal; i < iVectorLength; i++)
466 vecRet[iCurPos++] = operator[](i);
467
468 for (i = 0; i < iTo; i++)
469 vecRet[iCurPos++] = operator[](i);
470
471 return vecRet;
472 }
473 else
474 {
475 CMatlibVector<T> vecRet(iTo - iStartVal, VTY_TEMP);
476
477 for (i = iStartVal; i < iTo; i++)
478 vecRet[i - iStartVal] = operator[](i);
479
480 return vecRet;
481 }
482 }
483
484 template<class T> inline
operator()485 CMatlibVector<T> CMatlibVector<T>::operator()(const int iFrom,
486 const int iStep,
487 const int iTo) const
488 {
489 CMatlibVector<T> vecRet(size_t(abs(float(iTo - iFrom)) / abs(float(iStep))) + 1, VTY_TEMP);
490 int iOutPos = 0;
491 int i;
492
493 if (iFrom > iTo)
494 {
495 const int iEnd = iTo - 2;
496 for (i = iFrom - 1; i > iEnd; i += iStep)
497 {
498 vecRet[iOutPos] = operator[](i);
499 iOutPos++;
500 }
501 }
502 else
503 {
504 for (i = iFrom - 1; i < iTo; i += iStep)
505 {
506 vecRet[iOutPos] = operator[](i);
507 iOutPos++;
508 }
509 }
510
511 return vecRet;
512 }
513
514 template<class T> inline
PutIn(const int iFrom,const int iTo,CMatlibVector<T> & vecI)515 CMatlibVector<T>& CMatlibVector<T>::PutIn(const int iFrom,
516 const int iTo,
517 CMatlibVector<T>& vecI)
518 {
519 const int iStart = iFrom - 1;
520 const int iEnd = iTo - iStart;
521
522 for (int i = 0; i < iEnd; i++)
523 operator[](i + iStart) = vecI[i];
524
525 return *this;
526 }
527
528 template<class T> inline
Merge(const CMatlibVector<T> & vecA,T & tB)529 CMatlibVector<T>& CMatlibVector<T>::Merge(const CMatlibVector<T>& vecA, T& tB)
530 {
531 const int iSizeA = vecA.GetSize();
532
533 for (int i = 0; i < iSizeA; i++)
534 operator[](i) = vecA[i];
535
536 operator[](iSizeA) = tB;
537
538 return *this;
539 }
540
541 template<class T> inline
Merge(const CMatlibVector<T> & vecA,const CMatlibVector<T> & vecB)542 CMatlibVector<T>& CMatlibVector<T>::Merge(const CMatlibVector<T>& vecA,
543 const CMatlibVector<T>& vecB)
544 {
545 int i;
546 const int iSizeA = vecA.GetSize();
547 const int iSizeB = vecB.GetSize();
548
549 /* Put first vector */
550 for (i = 0; i < iSizeA; i++)
551 operator[](i) = vecA[i];
552
553 /* Put second vector behind the first one, both
554 together must have length of *this */
555 for (i = 0; i < iSizeB; i++)
556 operator[](i + iSizeA) = vecB[i];
557
558 return *this;
559 }
560
561 template<class T> inline
Merge(const CMatlibVector<T> & vecA,const CMatlibVector<T> & vecB,const CMatlibVector<T> & vecC)562 CMatlibVector<T>& CMatlibVector<T>::Merge(const CMatlibVector<T>& vecA,
563 const CMatlibVector<T>& vecB,
564 const CMatlibVector<T>& vecC)
565 {
566 int i;
567 const int iSizeA = vecA.GetSize();
568 const int iSizeB = vecB.GetSize();
569 const int iSizeC = vecC.GetSize();
570 const int iSizeAB = iSizeA + iSizeB;
571
572 /* Put first vector */
573 for (i = 0; i < iSizeA; i++)
574 operator[](i) = vecA[i];
575
576 /* Put second vector behind the first one */
577 for (i = 0; i < iSizeB; i++)
578 operator[](i + iSizeA) = vecB[i];
579
580 /* Put third vector behind previous put vectors */
581 for (i = 0; i < iSizeC; i++)
582 operator[](i + iSizeAB) = vecC[i];
583
584 return *this;
585 }
586
587
588 /******************************************************************************/
589 /* CMatlibMatrix class ********************************************************/
590 /******************************************************************************/
591 /*
592 We define: Matrix[row][column]
593 */
594 template<class T>
595 class CMatlibMatrix
596 {
597 public:
598 /* Construction, Destruction -------------------------------------------- */
CMatlibMatrix()599 CMatlibMatrix() : eVType(VTY_CONST), iRowSize(0), ppData(NULL) {}
600 CMatlibMatrix(const int iNRowLen, const int iNColLen,
eVType(eNTy)601 const EVecTy eNTy = VTY_CONST) : eVType(eNTy),
602 iRowSize(0), ppData(NULL) {Init(iNRowLen, iNColLen);}
CMatlibMatrix(const int iNRowLen,const int iNColLen,const T tIniVal)603 CMatlibMatrix(const int iNRowLen, const int iNColLen, const T tIniVal) :
604 eVType(VTY_CONST), iRowSize(0), ppData(NULL)
605 {Init(iNRowLen, iNColLen, tIniVal);}
606 CMatlibMatrix(const CMatlibMatrix<T>& matI);
607
~CMatlibMatrix()608 virtual ~CMatlibMatrix() {if (ppData != NULL) delete[] ppData;}
609
610 void Init(const int iNRowLen, const int iNColLen, const T tIniVal = 0);
GetRowSize()611 inline int GetRowSize() const {return iRowSize;}
GetColSize()612 inline int GetColSize() const
613 {if (iRowSize > 0) return ppData[0].GetSize(); else return 0;}
614
615 /* Operator[] (Regular indices!!!) */
616 inline CMatlibVector<T>& operator[](int const iPos) const
617 {_TESTRNGRM(iPos); return ppData[iPos];}
618 inline CMatlibVector<T>& operator[](int const iPos)
619 {_TESTRNGWM(iPos); return ppData[iPos];} // For use as l value
620
621 /* Operator() */
operator()622 inline CMatlibVector<T>& operator()(int const iPos) const
623 {_TESTRNGRM(iPos - 1); return ppData[iPos - 1];}
operator()624 inline CMatlibVector<T>& operator()(int const iPos)
625 {_TESTRNGWM(iPos - 1); return ppData[iPos - 1];} // For use as l value
626
627 CMatlibMatrix<T> operator()(const int iRowFrom, const int iRowTo,
628 const int iColFrom, const int iColTo) const;
629
630 /* operator= */
631 inline CMatlibMatrix<T>& operator=(const CMatlibMatrix<CReal>& matI)
632 {_TESTSIZEM(matI.GetRowSize()); _MATOPCL(= matI[i]);}
633 inline CMatlibMatrix<CComplex>& operator=(const CMatlibMatrix<CComplex>& matI)
634 {_TESTSIZEM(matI.GetRowSize()); _MATOPCL(= matI[i]);}
635
636 /* operator+= */
637 inline CMatlibMatrix<T>& operator+=(const CMatlibMatrix<CReal>& matI)
638 {_MATOPCL(+= matI[i]);}
639 inline CMatlibMatrix<CComplex>& operator+=(const CMatlibMatrix<CComplex>& matI)
640 {_MATOPCL(+= matI[i]);}
641
642 /* operator-= */
643 inline CMatlibMatrix<T>& operator-=(const CMatlibMatrix<CReal>& matI)
644 {_MATOPCL(-= matI[i]);}
645 inline CMatlibMatrix<CComplex>& operator-=(const CMatlibMatrix<CComplex>& matI)
646 {_MATOPCL(-= matI[i]);}
647
648 /* operator*= */
649 inline CMatlibMatrix<T>& operator*=(const CReal& rI)
650 {_MATOPCL(*= rI);}
651 inline CMatlibMatrix<CComplex>& operator*=(const CComplex& cI)
652 {_MATOPCL(*= cI);}
653
654 /* operator/= */
655 inline CMatlibMatrix<T>& operator/=(const CReal& rI)
656 {_MATOPCL(/= rI);}
657 inline CMatlibMatrix<CComplex>& operator/=(const CComplex& cI)
658 {_MATOPCL(/= cI);}
659
660 protected:
661 EVecTy eVType;
662 int iRowSize;
663 CMatlibVector<T>* ppData;
664 };
665
666
667 /* Help functions *************************************************************/
668 /* operator+ */
669 inline CMatlibMatrix<CComplex> // cm, cm
670 operator+(const CMatlibMatrix<CComplex>& cmA, const CMatlibMatrix<CComplex>& cmB)
671 {
672 const int iRowSizeA = cmA.GetRowSize();
673 const int iColSizeA = cmA.GetColSize();
674 CMatlibMatrix<CComplex> matRet(iRowSizeA, iColSizeA, VTY_TEMP);
675
676 for (int j = 0; j < iRowSizeA; j++)
677 {
678 for (int i = 0; i < iColSizeA; i++)
679 matRet[j][i] = cmA[j][i] + cmB[j][i];
680 }
681
682 return matRet;
683 }
684
685 /* operator- */
686 inline CMatlibMatrix<CComplex> // cm, cm
687 operator-(const CMatlibMatrix<CComplex>& cmA, const CMatlibMatrix<CComplex>& cmB)
688 {
689 const int iRowSizeA = cmA.GetRowSize();
690 const int iColSizeA = cmA.GetColSize();
691 CMatlibMatrix<CComplex> matRet(iRowSizeA, iColSizeA, VTY_TEMP);
692
693 for (int j = 0; j < iRowSizeA; j++)
694 {
695 for (int i = 0; i < iColSizeA; i++)
696 matRet[j][i] = cmA[j][i] - cmB[j][i];
697 }
698
699 return matRet;
700 }
701
702 /* operator* */
703 inline CMatlibVector<CComplex> // cm, cv
704 operator*(const CMatlibMatrix<CComplex>& cmA, const CMatlibVector<CComplex>& cvB)
705 {
706 const int iRowSizeA = cmA.GetRowSize();
707 const int iSizeB = cvB.GetSize();
708 CMatlibVector<CComplex> vecRet(iSizeB, VTY_TEMP);
709
710 for (int j = 0; j < iRowSizeA; j++)
711 {
712 vecRet[j] = (CReal) 0.0;
713
714 for (int i = 0; i < iSizeB; i++)
715 vecRet[j] += cmA[j][i] * cvB[i];
716 }
717
718 return vecRet;
719 }
720
721 /* operator* */
722 inline CMatlibVector<CReal> // rm, rv
723 operator*(const CMatlibMatrix<CReal>& rmA, const CMatlibVector<CReal>& rvB)
724 {
725 const int iRowSizeA = rmA.GetRowSize();
726 const int iSizeB = rvB.GetSize();
727 CMatlibVector<CReal> vecRet(iSizeB, VTY_TEMP);
728
729 for (int j = 0; j < iRowSizeA; j++)
730 {
731 vecRet[j] = (CReal) 0.0;
732
733 for (int i = 0; i < iSizeB; i++)
734 vecRet[j] += rmA[j][i] * rvB[i];
735 }
736
737 return vecRet;
738 }
739
740 /* operator* */
741 inline CMatlibMatrix<CComplex> // cm, cm
742 operator*(const CMatlibMatrix<CComplex>& cmA, const CMatlibMatrix<CComplex>& cmB)
743 {
744 const int iRowSizeA = cmA.GetRowSize();
745 const int iRowSizeB = cmB.GetRowSize();
746 const int iColSizeB = cmB.GetColSize();
747 CMatlibMatrix<CComplex> matRet(iRowSizeA, iColSizeB, VTY_TEMP);
748
749 for (int k = 0; k < iColSizeB; k++)
750 {
751 for (int j = 0; j < iRowSizeA; j++)
752 {
753 matRet[j][k] = (CReal) 0.0;
754
755 for (int i = 0; i < iRowSizeB; i++)
756 matRet[j][k] += cmA[j][i] * cmB[i][k];
757 }
758 }
759
760 return matRet;
761 }
762
763 /* operator* */
764 inline CMatlibMatrix<CComplex> // c, cm
765 operator*(const CComplex& cA, const CMatlibMatrix<CComplex>& cmB)
766 {
767 const int iRowSizeB = cmB.GetRowSize();
768 const int iColSizeB = cmB.GetColSize();
769 CMatlibMatrix<CComplex> matRet(iRowSizeB, iColSizeB, VTY_TEMP);
770
771 for (int k = 0; k < iColSizeB; k++)
772 {
773 for (int j = 0; j < iRowSizeB; j++)
774 matRet[j][k] = cA * cmB[j][k];
775 }
776
777 return matRet;
778 }
779
780 /* operator* */
781 inline CMatlibMatrix<CReal> // r, rm
782 operator*(const CReal& rA, const CMatlibMatrix<CReal>& rmB)
783 {
784 const int iRowSizeB = rmB.GetRowSize();
785 const int iColSizeB = rmB.GetColSize();
786 CMatlibMatrix<CReal> matRet(iRowSizeB, iColSizeB, VTY_TEMP);
787
788 for (int k = 0; k < iColSizeB; k++)
789 {
790 for (int j = 0; j < iRowSizeB; j++)
791 matRet[j][k] = rA * rmB[j][k];
792 }
793
794 return matRet;
795 }
796
797
798 /* Implementation **************************************************************
799 (the implementation of template classes must be in the header file!) */
800 template<class T>
CMatlibMatrix(const CMatlibMatrix<T> & matI)801 CMatlibMatrix<T>::CMatlibMatrix(const CMatlibMatrix<T>& matI) :
802 eVType(VTY_CONST), iRowSize(matI.GetRowSize()), ppData(NULL)
803 {
804 if (iRowSize > 0)
805 {
806 /* Allocate data block for vector */
807 ppData = new CMatlibVector<T>[iRowSize];
808
809 /* Init column vectors and copy */
810 for (int i = 0; i < iRowSize; i++)
811 {
812 ppData[i].Init(matI.GetColSize());
813
814 /* Copy entire vector */
815 ppData[i] = matI[i];
816 }
817 }
818 }
819
820 template<class T>
Init(const int iNRowLen,const int iNColLen,const T tIniVal)821 void CMatlibMatrix<T>::Init(const int iNRowLen, const int iNColLen, const T tIniVal)
822 {
823 iRowSize = iNRowLen;
824
825 /* Allocate data block for vector */
826 if (iRowSize > 0)
827 {
828 if (ppData != NULL)
829 delete[] ppData;
830
831 ppData = new CMatlibVector<T>[iRowSize];
832
833 /* Init column vectors and set to init value */
834 for (int i = 0; i < iRowSize; i++)
835 ppData[i].Init(iNColLen, tIniVal);
836 }
837 }
838
839 template<class T> inline
operator()840 CMatlibMatrix<T> CMatlibMatrix<T>::operator()(const int iRowFrom, const int iRowTo,
841 const int iColFrom, const int iColTo) const
842 {
843 const int iStartRow = iRowFrom - 1;
844 const int iStartCol = iColFrom - 1;
845 CMatlibMatrix<T> matRet(iRowTo - iStartRow, iColTo - iStartCol, VTY_TEMP);
846
847 for (int j = iStartRow; j < iRowTo; j++)
848 {
849 for (int i = iStartCol; i < iColTo; i++)
850 matRet[j - iStartRow][i - iStartCol] = operator[](j)[i];
851 }
852
853 return matRet;
854 }
855
856
857 /* Include toolboxes after all type definitions */
858 #include "MatlibStdToolbox.h"
859 #include "MatlibSigProToolbox.h"
860
861
862 #endif /* _MATLIB_H_ */
863