Subversion Repositories gelsvn

Rev

Rev 89 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 89 Rev 119
1
#if !defined(MATRIX_H_HAA_AGUST_2001)
1
#if !defined(MATRIX_H_HAA_AGUST_2001)
2
#define MATRIX_H_HAA_AGUST_2001
2
#define MATRIX_H_HAA_AGUST_2001
3
 
3
 
4
#include <cassert>
4
#include <cassert>
5
#include <iostream>
5
#include <iostream>
6
#include <cmath>
6
#include <cmath>
7
#include "Vector.h"
7
#include "Vector.h"
8
#include "CGLA/ArithMatFloat.h"
8
#include "CGLA/ArithMatFloat.h"
9
 
9
 
10
namespace LinAlg
10
namespace LinAlg
11
{
11
{
12
 
12
 
13
	/*!
13
	/*!
14
		\file Matrix.h
14
		\file Matrix.h
15
		\brief The Matrix type
15
		\brief The Matrix type
16
	*/
16
	*/
17
 
17
 
18
 
18
 
19
	/*!
19
	/*!
20
		\brief The Matrix type.
20
		\brief The Matrix type.
21
 
21
 
22
		This matrix teplate is one of the basic linear algebra types. In 
22
		This matrix teplate is one of the basic linear algebra types. In 
23
		principle it an overloaded Array of type T. The typedef Matrix with 
23
		principle it an overloaded Array of type T. The typedef Matrix with 
24
		T set to double
24
		T set to double
25
 
25
 
26
		typedef CMatrixType<double> CMatrix;
26
		typedef CMatrixType<double> CMatrix;
27
 
27
 
28
		is the delclaration intended to be used as the matrix type outside 
28
		is the delclaration intended to be used as the matrix type outside 
29
		the this and related files. The reson for making a template is 
29
		the this and related files. The reson for making a template is 
30
		twofold, first it allows for an esay change of precision in this
30
		twofold, first it allows for an esay change of precision in this
31
		linear algebre package. Secondly it allows for other MAtrix types 
31
		linear algebre package. Secondly it allows for other MAtrix types 
32
		in special cases.
32
		in special cases.
33
 
33
 
34
		The data structure consists of the array of type T, T* Data, 
34
		The data structure consists of the array of type T, T* Data, 
35
		with the Ilife vector, T** Ilife, refering to the begining 
35
		with the Ilife vector, T** Ilife, refering to the begining 
36
		of each of the columns. Hereby a coherent data area is avalable and 
36
		of each of the columns. Hereby a coherent data area is avalable and 
37
		the concept of 2-dimensionality is acheived by the Ilife vector. The 
37
		the concept of 2-dimensionality is acheived by the Ilife vector. The 
38
		pointer to the data area can be obtained by operator [] (int i) 
38
		pointer to the data area can be obtained by operator [] (int i) 
39
		with i=0, e.g. 
39
		with i=0, e.g. 
40
 
40
 
41
		Matrix A;
41
		Matrix A;
42
 
42
 
43
		T* pData=A[0];  
43
		T* pData=A[0];  
44
 
44
 
45
		will  make pData point to the data.
45
		will  make pData point to the data.
46
 
46
 
47
		It should be noted that much of the functionality associated with this
47
		It should be noted that much of the functionality associated with this
48
		Matrix type and the associated Vector type is located in the 
48
		Matrix type and the associated Vector type is located in the 
49
		LapackFunc.h file.
49
		LapackFunc.h file.
50
 
50
 
51
		\see LapackFunc.h
51
		\see LapackFunc.h
52
		\see CVectorType  
52
		\see CVectorType  
53
		\see Additional matrix operators (Located in this file Matrix.h).
53
		\see Additional matrix operators (Located in this file Matrix.h).
54
		\see CMatrix
54
		\see CMatrix
55
		\see CVector
55
		\see CVector
56
		\see CVec2Type
56
		\see CVec2Type
57
		\see CVec3Type
57
		\see CVec3Type
58
		\see CVec2
58
		\see CVec2
59
		\see CVec3
59
		\see CVec3
60
 
60
 
61
		\author Henrik Aanæs
61
		\author Henrik Aanæs
62
		\version Aug 2001
62
		\version Aug 2001
63
	*/
63
	*/
64
	template <class T>
64
	template <class T>
65
	class CMatrixType 
65
	class CMatrixType 
66
	{
66
	{
67
		friend class CVectorType<T>;
67
		friend class CVectorType<T>;
68
 
68
 
69
	private:
69
	private:
70
		///The number of rows in the matrix
70
		///The number of rows in the matrix
71
		int nRows;
71
		int nRows;
72
		///The number of columns in the matrix
72
		///The number of columns in the matrix
73
		int nCols;
73
		int nCols;
74
		/// The number of elements in the matrix, i.e. nElems=nRows*nCols.
74
		/// The number of elements in the matrix, i.e. nElems=nRows*nCols.
75
		int nElems;
75
		int nElems;
76
		///The pointer to the data containing the elements of the matrix.
76
		///The pointer to the data containing the elements of the matrix.
77
		T* Data;
77
		T* Data;
78
		///The Ilife vector containing pointers to the start of the columns. That is pointers to elements in Data.
78
		///The Ilife vector containing pointers to the start of the columns. That is pointers to elements in Data.
79
		T** Ilife;
79
		T** Ilife;
80
	
80
	
81
		///Sets all the elements in the matrix to Scalar.
81
		///Sets all the elements in the matrix to Scalar.
82
		void SetScalar(const T& Scalar) 
82
		void SetScalar(const T& Scalar) 
83
		{
83
		{
84
			T* Itt=Data;
84
			T* Itt=Data;
85
			T* Stop=&(Data[nElems]);
85
			T* Stop=&(Data[nElems]);
86
			while(Itt!=Stop)
86
			while(Itt!=Stop)
87
				{
87
				{
88
					*Itt=Scalar;
88
					*Itt=Scalar;
89
					++Itt;
89
					++Itt;
90
				}
90
				}
91
		} 
91
		} 
92
	
92
	
93
		///Clears the allocated memory.
93
		///Clears the allocated memory.
94
		void CleanUp(void)
94
		void CleanUp(void)
95
		{
95
		{
96
			if(Data!=NULL) 
96
			if(Data!=NULL) 
97
				delete [] Data; 
97
				delete [] Data; 
98
			if(Ilife!=NULL) 
98
			if(Ilife!=NULL) 
99
				delete [] Ilife;
99
				delete [] Ilife;
100
		};
100
		};
101
	
101
	
102
		///Allocates memory. Note that nRows,nCols and nElems should be set first.
102
		///Allocates memory. Note that nRows,nCols and nElems should be set first.
103
		void Allocate(void)
103
		void Allocate(void)
104
		{ 
104
		{ 
105
			Data= new T[nElems]; 
105
			Data= new T[nElems]; 
106
			Ilife=new T*[nRows];	
106
			Ilife=new T*[nRows];	
107
			for(int cRow=0;cRow<nRows;cRow++) 
107
			for(int cRow=0;cRow<nRows;cRow++) 
108
				Ilife[cRow]=&Data[cRow*nCols];
108
				Ilife[cRow]=&Data[cRow*nCols];
109
		};
109
		};
110
	
110
	
111
	public:
111
	public:
112
		///	\name The Constructors and Destructors of the CMatrixType class.
112
		///	\name The Constructors and Destructors of the CMatrixType class.
113
		//@{
113
		//@{
114
		CMatrixType<T>() : nRows(0),nCols(0),nElems(0),Data(NULL),Ilife(NULL) {};
114
		CMatrixType<T>() : nRows(0),nCols(0),nElems(0),Data(NULL),Ilife(NULL) {};
115
 
115
 
116
		CMatrixType<T>(const int Row,const int Col) : nRows(Row),nCols(Col),nElems(Row*Col) {Allocate();}
116
		CMatrixType<T>(const int Row,const int Col) : nRows(Row),nCols(Col),nElems(Row*Col) {Allocate();}
117
 
117
 
118
		CMatrixType<T>(const int Row,const int Col,const T& Scalar) : nRows(Row),nCols(Col),nElems(Row*Col) {Allocate(); SetScalar(Scalar);};
118
		CMatrixType<T>(const int Row,const int Col,const T& Scalar) : nRows(Row),nCols(Col),nElems(Row*Col) {Allocate(); SetScalar(Scalar);};
119
 
119
 
120
		CMatrixType<T>(const CMatrixType<T>& Rhs) : nRows(Rhs.nRows),nCols(Rhs.nCols),nElems(Rhs.nRows*Rhs.nCols) 
120
		CMatrixType<T>(const CMatrixType<T>& Rhs) : nRows(Rhs.nRows),nCols(Rhs.nCols),nElems(Rhs.nRows*Rhs.nCols) 
121
		{ 
121
		{ 
122
			Allocate();
122
			Allocate();
123
			memcpy(Data,Rhs.Data,nElems*sizeof(T));
123
			memcpy(Data,Rhs.Data,nElems*sizeof(T));
124
		}
124
		}
125
 
125
 
126
		template <class VVT, class HVT, class MT, unsigned int ROWS>
126
		template <class VVT, class HVT, class MT, unsigned int ROWS>
127
		CMatrixType<T>(const CGLA::ArithMatFloat<VVT,HVT,MT,ROWS>& Rhs): 
127
		CMatrixType<T>(const CGLA::ArithMatFloat<VVT,HVT,MT,ROWS>& Rhs): 
128
			nRows(Rhs.get_v_dim()),nCols(Rhs.get_h_dim()),
128
			nRows(Rhs.get_v_dim()),nCols(Rhs.get_h_dim()),
129
			nElems(nRows*nCols) 
129
			nElems(nRows*nCols) 
130
		{ 
130
		{ 
131
			Allocate();
131
			Allocate();
132
			for(int i=0;i<nElems;++i)
132
			for(int i=0;i<nElems;++i)
133
				Data[i] = Rhs.get()[i];
133
				Data[i] = Rhs.get()[i];
134
		}
134
		}
135
 
135
 
136
		CMatrixType<T>(const CVectorType<T>& Rhs) : nRows(Rhs.Length()),nCols(1),nElems(Rhs.Length()) 
136
		CMatrixType<T>(const CVectorType<T>& Rhs) : nRows(Rhs.Length()),nCols(1),nElems(Rhs.Length()) 
137
		{ 
137
		{ 
138
			Allocate();
138
			Allocate();
139
			memcpy(Data,Rhs.Data,nElems*sizeof(T));
139
			memcpy(Data,Rhs.Data,nElems*sizeof(T));
140
		}
140
		}
141
		virtual ~CMatrixType<T>(){CleanUp();};
141
		virtual ~CMatrixType<T>(){CleanUp();};
142
		//@}
142
		//@}
143
 
143
 
144
		/*!
144
		/*!
145
			\name	Asignment opertors.
145
			\name	Asignment opertors.
146
			The assignment operators of the CMatrixType class.
146
			The assignment operators of the CMatrixType class.
147
			Note that the matrix is resized if it does not have the correct size.
147
			Note that the matrix is resized if it does not have the correct size.
148
		*/
148
		*/
149
		//@{
149
		//@{
150
		CMatrixType<T>& operator=(const T& Rhs){ SetScalar(Rhs);return *this;}
150
		CMatrixType<T>& operator=(const T& Rhs){ SetScalar(Rhs);return *this;}
151
 
151
 
152
    CMatrixType<T>& operator=(const CMatrixType<T>& Rhs)
152
    CMatrixType<T>& operator=(const CMatrixType<T>& Rhs)
153
		{ 
153
		{ 
154
			Resize(Rhs.nRows,Rhs.nCols); 
154
			Resize(Rhs.nRows,Rhs.nCols); 
155
			memcpy(Data,Rhs.Data,nElems*sizeof(T));
155
			memcpy(Data,Rhs.Data,nElems*sizeof(T));
156
			return *this;
156
			return *this;
157
		}
157
		}
158
 
158
 
159
		template <class VVT, class HVT, class MT, unsigned int ROWS>
159
		template <class VVT, class HVT, class MT, unsigned int ROWS>
160
		CMatrixType<T>& operator=(const CGLA::ArithMatFloat<VVT,HVT,MT,ROWS>& Rhs)
160
		CMatrixType<T>& operator=(const CGLA::ArithMatFloat<VVT,HVT,MT,ROWS>& Rhs)
161
		{ 
161
		{ 
162
			Resize(Rhs.get_v_dim(),Rhs.get_h_dim()); 
162
			Resize(Rhs.get_v_dim(),Rhs.get_h_dim()); 
163
			for(int i=0;i<nElems;++i)
163
			for(int i=0;i<nElems;++i)
164
				Data[i] = Rhs.get()[i];
164
				Data[i] = Rhs.get()[i];
165
			return *this;
165
			return *this;
166
		}
166
		}
167
 
167
 
168
		CMatrixType<T>& operator=(const CVectorType<T>&Rhs)
168
		CMatrixType<T>& operator=(const CVectorType<T>&Rhs)
169
		{
169
		{
170
			Resize(Rhs.Lenght(),1);
170
			Resize(Rhs.Lenght(),1);
171
			memcpy(Data,Rhs.Dat,nElems*sizeof(T));
171
			memcpy(Data,Rhs.Dat,nElems*sizeof(T));
172
			return *this;
172
			return *this;
173
		}
173
		}
174
		//@}
174
		//@}
175
 
175
 
176
 
176
 
177
		/*!
177
		/*!
178
			\name Dimension functions.
178
			\name Dimension functions.
179
			Functions dealing with the dimension of the matrix.
179
			Functions dealing with the dimension of the matrix.
180
		*/
180
		*/
181
		//@{
181
		//@{
182
		/// Returns the number of rows in the matrix
182
		/// Returns the number of rows in the matrix
183
		const int Rows(void) const {return nRows;};
183
		const int Rows(void) const {return nRows;};
184
		///Retruns the number of columns in the matrix. 
184
		///Retruns the number of columns in the matrix. 
185
		const int Cols(void) const {return nCols;};
185
		const int Cols(void) const {return nCols;};
186
		/// Resizes the matrix IF it does not already have the desired dimensions.
186
		/// Resizes the matrix IF it does not already have the desired dimensions.
187
		void Resize(const int Row,const int Col)
187
		void Resize(const int Row,const int Col)
188
		{
188
		{
189
			if(Row!=nRows || Col!=nCols)
189
			if(Row!=nRows || Col!=nCols)
190
				{
190
				{
191
					CleanUp();
191
					CleanUp();
192
					nRows=Row;
192
					nRows=Row;
193
					nCols=Col;
193
					nCols=Col;
194
					nElems=Row*Col;
194
					nElems=Row*Col;
195
					Data= new T[nElems];
195
					Data= new T[nElems];
196
					Ilife=new T*[nRows];
196
					Ilife=new T*[nRows];
197
					for(int cRow=0;cRow<nRows;cRow++)
197
					for(int cRow=0;cRow<nRows;cRow++)
198
						{
198
						{
199
							Ilife[cRow]=&Data[cRow*nCols];
199
							Ilife[cRow]=&Data[cRow*nCols];
200
						}
200
						}
201
				}
201
				}
202
		}
202
		}
203
		//@}
203
		//@}
204
	
204
	
205
		/*!
205
		/*!
206
			\name Acces functions and operators.
206
			\name Acces functions and operators.
207
			The [] operators are the ones intaended for usual use. The get 
207
			The [] operators are the ones intaended for usual use. The get 
208
			and set functions are added to allow for genral Matrix function 
208
			and set functions are added to allow for genral Matrix function 
209
			templates.
209
			templates.
210
		*/
210
		*/
211
		//@{
211
		//@{
212
		T* operator[](const int Row) 
212
		T* operator[](const int Row) 
213
		{
213
		{
214
			assert(Row>=0);
214
			assert(Row>=0);
215
			assert(Row<nRows);
215
			assert(Row<nRows);
216
		
216
		
217
			return Ilife[Row];
217
			return Ilife[Row];
218
		}
218
		}
219
	
219
	
220
		const T* operator[](const int Row) const 
220
		const T* operator[](const int Row) const 
221
		{
221
		{
222
			assert(Row>=0);
222
			assert(Row>=0);
223
			assert(Row<nRows);
223
			assert(Row<nRows);
224
		
224
		
225
			return Ilife[Row];
225
			return Ilife[Row];
226
		}
226
		}
227
 
227
 
228
		const T& get(const int Row,const int Col) const
228
		const T& get(const int Row,const int Col) const
229
		{
229
		{
230
			assert(Row>=0 && Row<nRows);
230
			assert(Row>=0 && Row<nRows);
231
			assert(Col>=0 && Col<nCols);
231
			assert(Col>=0 && Col<nCols);
232
			return Ilife[Row][Col];
232
			return Ilife[Row][Col];
233
		}
233
		}
234
 
234
 
235
		void set(const int Row,const int Col,const T val) 
235
		void set(const int Row,const int Col,const T val) 
236
		{
236
		{
237
			assert(Row>=0 && Row<nRows);
237
			assert(Row>=0 && Row<nRows);
238
			assert(Col>=0 && Col<nCols);
238
			assert(Col>=0 && Col<nCols);
239
			Ilife[Row][Col]=val;
239
			Ilife[Row][Col]=val;
240
		}
240
		}
241
		//@}
241
		//@}
242
	
242
	
243
		/*! 
243
		/*! 
244
			\name Arithmetic operators.
244
			\name Arithmetic operators.
245
			The aritmic operators of the CMatrixType class.
245
			The aritmic operators of the CMatrixType class.
246
		*/
246
		*/
247
		//@{
247
		//@{
248
		CMatrixType<T> operator+(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret+=Rhs;};
248
		CMatrixType<T> operator+(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret+=Rhs;};
249
		CMatrixType<T> operator-(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret-=Rhs;};
249
		CMatrixType<T> operator-(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret-=Rhs;};
250
		CMatrixType<T> operator*(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret*=Rhs;};
250
		CMatrixType<T> operator*(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret*=Rhs;};
251
		CMatrixType<T> operator/(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret/=Rhs;};
251
		CMatrixType<T> operator/(const T& Rhs) const {CMatrixType<T>Ret(*this); return Ret/=Rhs;};
252
		CMatrixType<T>& operator+=(const T& Rhs)
252
		CMatrixType<T>& operator+=(const T& Rhs)
253
		{
253
		{
254
 
254
 
255
			T* Itt=Data;
255
			T* Itt=Data;
256
			T* Stop=&(Data[nElems]);
256
			T* Stop=&(Data[nElems]);
257
			while(Itt!=Stop)
257
			while(Itt!=Stop)
258
				{
258
				{
259
					(*Itt)+=Rhs;
259
					(*Itt)+=Rhs;
260
					++Itt;
260
					++Itt;
261
				}
261
				}
262
		
262
		
263
			return *this; 
263
			return *this; 
264
		}
264
		}
265
	
265
	
266
		CMatrixType<T>& operator-=(const T& Rhs)
266
		CMatrixType<T>& operator-=(const T& Rhs)
267
		{
267
		{
268
			T* Itt=Data;
268
			T* Itt=Data;
269
			T* Stop=&(Data[nElems]);
269
			T* Stop=&(Data[nElems]);
270
			while(Itt!=Stop)
270
			while(Itt!=Stop)
271
				{
271
				{
272
					(*Itt)-=Rhs;
272
					(*Itt)-=Rhs;
273
					++Itt;
273
					++Itt;
274
				}
274
				}
275
		
275
		
276
			return *this; 
276
			return *this; 
277
		}
277
		}
278
	
278
	
279
		CMatrixType<T>& operator*=(const T& Rhs)
279
		CMatrixType<T>& operator*=(const T& Rhs)
280
		{
280
		{
281
			T* Itt=Data;
281
			T* Itt=Data;
282
			T* Stop=&(Data[nElems]);
282
			T* Stop=&(Data[nElems]);
283
			while(Itt!=Stop)
283
			while(Itt!=Stop)
284
				{
284
				{
285
					(*Itt)*=Rhs;
285
					(*Itt)*=Rhs;
286
					++Itt;
286
					++Itt;
287
				}
287
				}
288
		
288
		
289
			return *this; 
289
			return *this; 
290
		}
290
		}
291
	
291
	
292
		CMatrixType<T>& operator/=(const T& Rhs)
292
		CMatrixType<T>& operator/=(const T& Rhs)
293
		{
293
		{
294
			T* Itt=Data;
294
			T* Itt=Data;
295
			T* Stop=&(Data[nElems]);
295
			T* Stop=&(Data[nElems]);
296
			while(Itt!=Stop)
296
			while(Itt!=Stop)
297
				{
297
				{
298
					(*Itt)/=Rhs;
298
					(*Itt)/=Rhs;
299
					++Itt;
299
					++Itt;
300
				}
300
				}
301
		
301
		
302
			return *this; 
302
			return *this; 
303
		}
303
		}
304
 
304
 
305
		CMatrixType<T> operator+(const CMatrixType<T>& Rhs) const {CMatrixType<T>Ret(*this); return Ret+=Rhs;};
305
		CMatrixType<T> operator+(const CMatrixType<T>& Rhs) const {CMatrixType<T>Ret(*this); return Ret+=Rhs;};
306
		CMatrixType<T> operator-(const CMatrixType<T>& Rhs) const {CMatrixType<T>Ret(*this); return Ret-=Rhs;};
306
		CMatrixType<T> operator-(const CMatrixType<T>& Rhs) const {CMatrixType<T>Ret(*this); return Ret-=Rhs;};
307
 
307
 
308
		CMatrixType<T>& operator+=(const CMatrixType<T>& Rhs)
308
		CMatrixType<T>& operator+=(const CMatrixType<T>& Rhs)
309
		{
309
		{
310
			assert(nRows==Rhs.nRows);
310
			assert(nRows==Rhs.nRows);
311
			assert(nCols==Rhs.nCols);
311
			assert(nCols==Rhs.nCols);
312
 
312
 
313
			T* IttRhs=Rhs.Data;
313
			T* IttRhs=Rhs.Data;
314
			T* Itt=Data;
314
			T* Itt=Data;
315
			T* Stop=&(Data[nElems]);
315
			T* Stop=&(Data[nElems]);
316
			while(Itt!=Stop)
316
			while(Itt!=Stop)
317
				{
317
				{
318
					(*Itt)+=(*IttRhs);
318
					(*Itt)+=(*IttRhs);
319
					++Itt;
319
					++Itt;
320
					++IttRhs;
320
					++IttRhs;
321
				}
321
				}
322
 
322
 
323
			return *this; 
323
			return *this; 
324
		}
324
		}
325
	
325
	
326
		CMatrixType<T>& operator-=(const CMatrixType<T>& Rhs)
326
		CMatrixType<T>& operator-=(const CMatrixType<T>& Rhs)
327
		{
327
		{
328
			assert(nRows==Rhs.nRows);
328
			assert(nRows==Rhs.nRows);
329
			assert(nCols==Rhs.nCols);
329
			assert(nCols==Rhs.nCols);
330
			T* IttRhs=Rhs.Data;
330
			T* IttRhs=Rhs.Data;
331
			T* Itt=Data;
331
			T* Itt=Data;
332
			T* Stop=&(Data[nElems]);
332
			T* Stop=&(Data[nElems]);
333
			while(Itt!=Stop)
333
			while(Itt!=Stop)
334
				{
334
				{
335
					(*Itt)-=(*IttRhs);
335
					(*Itt)-=(*IttRhs);
336
					++Itt;
336
					++Itt;
337
					++IttRhs;
337
					++IttRhs;
338
				}
338
				}
339
			return *this; 
339
			return *this; 
340
		}
340
		}
341
 
341
 
342
		CMatrixType<T> operator* (const CMatrixType<T>& Rhs) const
342
		CMatrixType<T> operator* (const CMatrixType<T>& Rhs) const
343
		{
343
		{
344
			assert(nCols==Rhs.nRows);
344
			assert(nCols==Rhs.nRows);
345
			CMatrixType<T> Ret(nRows,Rhs.nCols);
345
			CMatrixType<T> Ret(nRows,Rhs.nCols);
346
 
346
 
347
			for(int cRow=0;cRow<nRows;cRow++)
347
			for(int cRow=0;cRow<nRows;cRow++)
348
				{
348
				{
349
					for(int cCol=0;cCol<Rhs.nCols;cCol++)
349
					for(int cCol=0;cCol<Rhs.nCols;cCol++)
350
						{
350
						{
351
							T* Ret_ij=&(Ret[cRow][cCol]);
351
							T* Ret_ij=&(Ret[cRow][cCol]);
352
							*Ret_ij=0;
352
							*Ret_ij=0;
353
							for(int cI=0;cI<nCols;cI++)
353
							for(int cI=0;cI<nCols;cI++)
354
								{
354
								{
355
									(*Ret_ij)+=(Ilife[cRow])[cI]*Rhs[cI][cCol];
355
									(*Ret_ij)+=(Ilife[cRow])[cI]*Rhs[cI][cCol];
356
								}
356
								}
357
						}
357
						}
358
				}
358
				}
359
			return Ret;
359
			return Ret;
360
		}
360
		}
361
 
361
 
362
		CVectorType<T> operator* (const CVectorType<T>& Rhs) const 
362
		CVectorType<T> operator* (const CVectorType<T>& Rhs) const 
363
		{
363
		{
364
			assert(nCols==Rhs.nElems);
364
			assert(nCols==Rhs.nElems);
365
			CVectorType<T> Ret(nRows,0);
365
			CVectorType<T> Ret(nRows,0);
366
 
366
 
367
			T* ThisItt=Data;
367
			T* ThisItt=Data;
368
			T* Stop=&(Data[nElems]);
368
			T* Stop=&(Data[nElems]);
369
			T* RhsItt;
369
			T* RhsItt;
370
			T* RetItt=Ret.Data;
370
			T* RetItt=Ret.Data;
371
			while(ThisItt!=Stop)
371
			while(ThisItt!=Stop)
372
				{
372
				{
373
					RhsItt=Rhs.Data;
373
					RhsItt=Rhs.Data;
374
					for(int cCol=nCols;cCol>0;--cCol)
374
					for(int cCol=nCols;cCol>0;--cCol)
375
						{
375
						{
376
							(*RetItt)+=(*RhsItt)*(*ThisItt);
376
							(*RetItt)+=(*RhsItt)*(*ThisItt);
377
							++ThisItt;
377
							++ThisItt;
378
							++RhsItt;
378
							++RhsItt;
379
						}
379
						}
380
					++RetItt;
380
					++RetItt;
381
				}
381
				}
382
			return Ret;
382
			return Ret;
383
		}
383
		}
384
		//@}
384
		//@}
385
	
385
	
386
		/*!
386
		/*!
387
			\name Matrix Transpose.
387
			\name Matrix Transpose.
388
			Functions to transpose the matrix
388
			Functions to transpose the matrix
389
		*/
389
		*/
390
		//@{
390
		//@{
391
		///Transposes the matrix itself.
391
		///Transposes the matrix itself.
392
		CMatrixType<T>& Transpose(void)
392
		CMatrixType<T>& Transpose(void)
393
		{
393
		{
394
			T* NewData =new T[nElems];
394
			T* NewData =new T[nElems];
395
			T*RowP;
395
			T*RowP;
396
			int cRow;
396
			int cRow;
397
			for(cRow=0;cRow<nRows;cRow++)
397
			for(cRow=0;cRow<nRows;cRow++)
398
				{
398
				{
399
					RowP=Ilife[cRow];
399
					RowP=Ilife[cRow];
400
					for(int cCol=0;cCol<nCols;cCol++)
400
					for(int cCol=0;cCol<nCols;cCol++)
401
						{
401
						{
402
							NewData[cCol*nRows+cRow]=RowP[cCol];
402
							NewData[cCol*nRows+cRow]=RowP[cCol];
403
						}
403
						}
404
				}
404
				}
405
		
405
		
406
			CleanUp();
406
			CleanUp();
407
		
407
		
408
			Data=NewData;
408
			Data=NewData;
409
			int temp=nRows;
409
			int temp=nRows;
410
			nRows=nCols;
410
			nRows=nCols;
411
			nCols=temp;
411
			nCols=temp;
412
			Ilife=new T*[nRows];	
412
			Ilife=new T*[nRows];	
413
			for(cRow=0;cRow<nRows;cRow++) 
413
			for(cRow=0;cRow<nRows;cRow++) 
414
				{
414
				{
415
					Ilife[cRow]=&Data[cRow*nCols];
415
					Ilife[cRow]=&Data[cRow*nCols];
416
				}
416
				}
417
		
417
		
418
			return *this;
418
			return *this;
419
		}
419
		}
420
 
420
 
421
		///Returns the Transpose of the matrix is returned in AT.
421
		///Returns the Transpose of the matrix is returned in AT.
422
		void Transposed(CMatrixType<T>& AT) const
422
		void Transposed(CMatrixType<T>& AT) const
423
		{
423
		{
424
			AT.Resize(nCols,nRows);
424
			AT.Resize(nCols,nRows);
425
			for(int cRow=nRows-1;cRow>=0;--cRow)
425
			for(int cRow=nRows-1;cRow>=0;--cRow)
426
				{
426
				{
427
					for(int cCol=nCols-1;cCol>=0;--cCol)
427
					for(int cCol=nCols-1;cCol>=0;--cCol)
428
						{
428
						{
429
							AT[cCol][cRow]=Ilife[cRow][cCol];
429
							AT[cCol][cRow]=Ilife[cRow][cCol];
430
						}
430
						}
431
				}
431
				}
432
		}
432
		}
433
 
433
 
434
		///Returns the Transpose of the matrix is returned.
434
		///Returns the Transpose of the matrix is returned.
435
		CMatrixType<T> Transposed(void) const
435
		CMatrixType<T> Transposed(void) const
436
		{
436
		{
437
			CMatrixType<T> Ret(nCols,nRows);
437
			CMatrixType<T> Ret(nCols,nRows);
438
			Transposed(Ret);
438
			Transposed(Ret);
439
			return Ret;
439
			return Ret;
440
		}
440
		}
441
		//@}
441
		//@}
442
 
442
 
443
 
443
 
444
		/*!
444
		/*!
445
			\name Elementwise functions.
445
			\name Elementwise functions.
446
			Functions the perform some elementary operation on each 
446
			Functions the perform some elementary operation on each 
447
			of the elements of the matrix.
447
			of the elements of the matrix.
448
		*/
448
		*/
449
		//@{
449
		//@{
450
		///Pairwise multipling the elemtns of the two matrices. Equvivalent to MatLAb's .*
450
		///Pairwise multipling the elemtns of the two matrices. Equvivalent to MatLAb's .*
451
		void ElemMult(const CMatrixType<T>& Rhs)
451
		void ElemMult(const CMatrixType<T>& Rhs)
452
		{
452
		{
453
			T*Itt=Data;
453
			T*Itt=Data;
454
			T*Stop=&(Data[nElems]);
454
			T*Stop=&(Data[nElems]);
455
			T*RhsItt=Rhs.Data;
455
			T*RhsItt=Rhs.Data;
456
			while(Itt!=Stop)
456
			while(Itt!=Stop)
457
				{
457
				{
458
					(*Itt)*=(*RhsItt);
458
					(*Itt)*=(*RhsItt);
459
					++Itt;
459
					++Itt;
460
					++RhsItt;
460
					++RhsItt;
461
				}
461
				}
462
		}
462
		}
463
 
463
 
464
		///Pairwise dividnig the the elemtns of the two matrices. Equvivalent to MatLAb's ./
464
		///Pairwise dividnig the the elemtns of the two matrices. Equvivalent to MatLAb's ./
465
		void ElemDiv(const CMatrixType<T>& Rhs)
465
		void ElemDiv(const CMatrixType<T>& Rhs)
466
		{
466
		{
467
			T*Itt=Data;
467
			T*Itt=Data;
468
			T*Stop=&(Data[nElems]);
468
			T*Stop=&(Data[nElems]);
469
			T*RhsItt=Rhs.Data;
469
			T*RhsItt=Rhs.Data;
470
			while(Itt!=Stop)
470
			while(Itt!=Stop)
471
				{
471
				{
472
					(*Itt)/=(*RhsItt);
472
					(*Itt)/=(*RhsItt);
473
					++Itt;
473
					++Itt;
474
					++RhsItt;
474
					++RhsItt;
475
				}
475
				}
476
		}
476
		}
477
 
477
 
478
		///The elements of the matrix squared. Equvivalent to MatLAb's .^2
478
		///The elements of the matrix squared. Equvivalent to MatLAb's .^2
479
		void ElemSqr(void)
479
		void ElemSqr(void)
480
		{
480
		{
481
			T*Itt=Data;
481
			T*Itt=Data;
482
			T*Stop=&(Data[nElems]);
482
			T*Stop=&(Data[nElems]);
483
			while(Itt!=Stop)
483
			while(Itt!=Stop)
484
				{
484
				{
485
					(*Itt)*=(*Itt);
485
					(*Itt)*=(*Itt);
486
					++Itt;
486
					++Itt;
487
				}
487
				}
488
		}
488
		}
489
 
489
 
490
		///The square root of the elements. Equvivalent to MatLAb's .^0.5
490
		///The square root of the elements. Equvivalent to MatLAb's .^0.5
491
		void ElemSqrt(void)
491
		void ElemSqrt(void)
492
		{
492
		{
493
			T*Itt=Data;
493
			T*Itt=Data;
494
			T*Stop=&(Data[nElems]);
494
			T*Stop=&(Data[nElems]);
495
			while(Itt!=Stop)
495
			while(Itt!=Stop)
496
				{
496
				{
497
					(*Itt)=sqrt(*Itt);
497
					(*Itt)=sqrt(*Itt);
498
					++Itt;
498
					++Itt;
499
				}
499
				}
500
		}
500
		}
501
		//@}
501
		//@}
502
 
502
 
503
		/*!
503
		/*!
504
			\name Norms.
504
			\name Norms.
505
			Computes various norms of the matrix.
505
			Computes various norms of the matrix.
506
		*/
506
		*/
507
		//@{
507
		//@{
508
		///The Max norm returns the maximum absolute value
508
		///The Max norm returns the maximum absolute value
509
		const T NormMax(void) const
509
		const T NormMax(void) const
510
		{
510
		{
511
			T Ret=0;
511
			T Ret=0;
512
			T Abs;
512
			T Abs;
513
			const T*Itt=Data;
513
			const T*Itt=Data;
514
			const T*Stop=&(Data[nElems]);
514
			const T*Stop=&(Data[nElems]);
515
			while(Itt!=Stop)
515
			while(Itt!=Stop)
516
				{
516
				{
517
					Abs=(*Itt)>0?(*Itt):-(*Itt);
517
					Abs=(*Itt)>0?(*Itt):-(*Itt);
518
					if(Ret<Abs)
518
					if(Ret<Abs)
519
						Ret=Abs;
519
						Ret=Abs;
520
 
520
 
521
					++Itt;
521
					++Itt;
522
				}
522
				}
523
			return Ret;
523
			return Ret;
524
		}
524
		}
525
 
525
 
526
		/// The Frobenius norm returns the square root of the element sum of squares.
526
		/// The Frobenius norm returns the square root of the element sum of squares.
527
		const T NormFrobenius(void) const
527
		const T NormFrobenius(void) const
528
		{
528
		{
529
			T Ret=0;
529
			T Ret=0;
530
			const T*Itt=Data;
530
			const T*Itt=Data;
531
			const T*Stop=&(Data[nElems]);
531
			const T*Stop=&(Data[nElems]);
532
			while(Itt!=Stop)
532
			while(Itt!=Stop)
533
				{
533
				{
534
					Ret+=(*Itt)*(*Itt);
534
					Ret+=(*Itt)*(*Itt);
535
					++Itt;
535
					++Itt;
536
				}
536
				}
537
			return sqrt(Ret);
537
			return sqrt(Ret);
538
		}
538
		}
539
		//@}
539
		//@}
540
 
540
 
541
 
541
 
542
		/*!
542
		/*!
543
			\name Sub-part acces functions
543
			\name Sub-part acces functions
544
			These functions are used to acces sub parts of the matrix.
544
			These functions are used to acces sub parts of the matrix.
545
		*/
545
		*/
546
		//@{
546
		//@{
547
		///Get the cRow'th row and set it equal to Row. Equvivalent to Row=this(cRow,:) in Matlab.
547
		///Get the cRow'th row and set it equal to Row. Equvivalent to Row=this(cRow,:) in Matlab.
548
		void GetRow(CVectorType<T>& Row,const int cRow) const
548
		void GetRow(CVectorType<T>& Row,const int cRow) const
549
		{
549
		{
550
			Row.Resize(nCols);
550
			Row.Resize(nCols);
551
			memcpy(Row.Data,Ilife[cRow],nCols*sizeof(T));
551
			memcpy(Row.Data,Ilife[cRow],nCols*sizeof(T));
552
		}
552
		}
553
 
553
 
554
		///Sets the cRow'th row of the matrix equal to Row. Equvivalent to this(cRow,:)=Row in Matlab.
554
		///Sets the cRow'th row of the matrix equal to Row. Equvivalent to this(cRow,:)=Row in Matlab.
555
		void SetRow(CVectorType<T>& Row,const int cRow)
555
		void SetRow(CVectorType<T>& Row,const int cRow)
556
		{
556
		{
557
			assert(Row.Length()==nCols);
557
			assert(Row.Length()==nCols);
558
			memcpy(Ilife[cRow],Row.Data,nCols*sizeof(T));
558
			memcpy(Ilife[cRow],Row.Data,nCols*sizeof(T));
559
		}
559
		}
560
 
560
 
561
		///Copies row 'from' of this matrix to row 'to' of matrix A. Equvivalent to A(to,:)=this(from,:) in Matlab.
561
		///Copies row 'from' of this matrix to row 'to' of matrix A. Equvivalent to A(to,:)=this(from,:) in Matlab.
562
		void RowTransfere(CMatrixType<T>& A, const int from,const int to) const
562
		void RowTransfere(CMatrixType<T>& A, const int from,const int to) const
563
		{
563
		{
564
			assert(A.Cols()==nCols);
564
			assert(A.Cols()==nCols);
565
			memcpy(A.Ilife[to],Ilife[from],nCols*sizeof(T));
565
			memcpy(A.Ilife[to],Ilife[from],nCols*sizeof(T));
566
		}
566
		}
567
 
567
 
568
		void GetCol(CVectorType<T>&Col, const int nCol) const
568
		void GetCol(CVectorType<T>&Col, const int nCol) const
569
		{
569
		{
570
			Col.Resize(nRows);
570
			Col.Resize(nRows);
571
			const T* Itt=&Data[nCol];
571
			const T* Itt=&Data[nCol];
572
			T* VecItt=&Col[0];
572
			T* VecItt=&Col[0];
573
			for(int cRow=0;cRow<nRows;cRow++)
573
			for(int cRow=0;cRow<nRows;cRow++)
574
				{
574
				{
575
					*VecItt=*Itt;
575
					*VecItt=*Itt;
576
					++VecItt;
576
					++VecItt;
577
					Itt+=nCols;
577
					Itt+=nCols;
578
				}
578
				}
579
		}
579
		}
580
 
580
 
581
		void SetCol(CVectorType<T>&Col, const int nCol)
581
		void SetCol(CVectorType<T>&Col, const int nCol)
582
		{
582
		{
583
			assert(nRows==Col.Length());
583
			assert(nRows==Col.Length());
584
			assert(nCol<nCols);
584
			assert(nCol<nCols);
585
 
585
 
586
			T* Itt=&Data[nCol];
586
			T* Itt=&Data[nCol];
587
			T* VecItt=&Col[0];
587
			T* VecItt=&Col[0];
588
			for(int cRow=0;cRow<nRows;cRow++)
588
			for(int cRow=0;cRow<nRows;cRow++)
589
				{
589
				{
590
					*Itt=*VecItt;
590
					*Itt=*VecItt;
591
					++VecItt;
591
					++VecItt;
592
					Itt+=nCols;
592
					Itt+=nCols;
593
				}
593
				}
594
		}
594
		}
595
 
595
 
596
		void GetSubMatrix(CMatrixType<T>& A, const int StartRow, const int StartCol, const int RowDim, const int ColDim) const
596
		void GetSubMatrix(CMatrixType<T>& A, const int StartRow, const int StartCol, const int RowDim, const int ColDim) const
597
		{
597
		{
598
			assert(StartRow+RowDim<=nRows);
598
			assert(StartRow+RowDim<=nRows);
599
			assert(StartCol+ColDim<=nCols);
599
			assert(StartCol+ColDim<=nCols);
600
			A.Resize(RowDim,ColDim);
600
			A.Resize(RowDim,ColDim);
601
			const T* Itt=&(Ilife[StartRow][StartCol]);
601
			const T* Itt=&(Ilife[StartRow][StartCol]);
602
 
602
 
603
			for(int cRow=0;cRow<A.Rows();cRow++)
603
			for(int cRow=0;cRow<A.Rows();cRow++)
604
				{
604
				{
605
					memcpy(A[cRow],Itt,A.Cols()*sizeof(T));
605
					memcpy(A[cRow],Itt,A.Cols()*sizeof(T));
606
					Itt+=nCols;
606
					Itt+=nCols;
607
				}
607
				}
608
		}
608
		}
609
 
609
 
610
		void SetSubMatrix(CMatrixType<T>& A, const int StartRow, const int StartCol)
610
		void SetSubMatrix(CMatrixType<T>& A, const int StartRow, const int StartCol)
611
		{
611
		{
612
			assert(A.Rows()+StartRow<=nRows);
612
			assert(A.Rows()+StartRow<=nRows);
613
			assert(A.Cols()+StartCol<=nCols);
613
			assert(A.Cols()+StartCol<=nCols);
614
			T* Itt=&(Ilife[StartRow][StartCol]);
614
			T* Itt=&(Ilife[StartRow][StartCol]);
615
 
615
 
616
			for(int cRow=0;cRow<A.Rows();cRow++)
616
			for(int cRow=0;cRow<A.Rows();cRow++)
617
				{
617
				{
618
					memcpy(Itt,A[cRow],A.Cols()*sizeof(T));
618
					memcpy(Itt,A[cRow],A.Cols()*sizeof(T));
619
					Itt+=nCols;
619
					Itt+=nCols;
620
				}
620
				}
621
		}
621
		}
622
 
622
 
623
		//@}
623
		//@}
624
 
624
 
625
 
625
 
626
 
626
 
627
	};
627
	};
628
 
628
 
629
	/*!
629
	/*!
630
		\name Additional matrix operators
630
		\name Additional matrix operators
631
		These are operators heavily associated with CMAtrixType, but not included
631
		These are operators heavily associated with CMAtrixType, but not included
632
		in the class definition it self.
632
		in the class definition it self.
633
	*/
633
	*/
634
	//@{
634
	//@{
635
	template<class T>
635
	template<class T>
636
	inline CMatrixType<T> operator+(const T& Lhs,const CMatrixType<T>&Rhs ) 
636
	inline CMatrixType<T> operator+(const T& Lhs,const CMatrixType<T>&Rhs ) 
637
	{
637
	{
638
		return Rhs+Lhs;
638
		return Rhs+Lhs;
639
	}
639
	}
640
 
640
 
641
	template<class T>
641
	template<class T>
642
	inline CMatrixType<T> operator*(const T& Lhs,const CMatrixType<T>&Rhs ) 
642
	inline CMatrixType<T> operator*(const T& Lhs,const CMatrixType<T>&Rhs ) 
643
	{
643
	{
644
		return Rhs*Lhs;
644
		return Rhs*Lhs;
645
	}
645
	}
646
 
646
 
647
	template <class T>
647
	template <class T>
648
	std::ostream& operator<<(std::ostream &s, const CMatrixType<T> &A)
648
	std::ostream& operator<<(std::ostream &s, const CMatrixType<T> &A)
649
	{
649
	{
650
    int nRows=A.Rows();
650
    int nRows=A.Rows();
651
    int nCols=A.Cols();
651
    int nCols=A.Cols();
652
	
652
	
653
    for (int cRow=0; cRow<nRows; cRow++)
653
    for (int cRow=0; cRow<nRows; cRow++)
654
			{
654
			{
655
        for (int cCol=0; cCol<nCols; cCol++)
655
        for (int cCol=0; cCol<nCols; cCol++)
656
					{
656
					{
657
            s << A[cRow][cCol] << " ";
657
            s << A[cRow][cCol] << " ";
658
					}
658
					}
659
        s << "\n";
659
        s << "\n";
660
			}
660
			}
661
    return s;
661
    return s;
662
	}
662
	}
663
 
663
 
664
	template <class T>
664
	template <class T>
665
	std::istream& operator>>(std::istream &s, CMatrixType<T> &A)
665
	std::istream& operator>>(std::istream &s, CMatrixType<T> &A)
666
	{
666
	{
667
    int nRows;
667
    int nRows;
668
		int nCols;
668
		int nCols;
669
 
669
 
670
    s >> nRows >> nCols;
670
    s >> nRows >> nCols;
671
 
671
 
672
    if ( nRows!=A.Rows() || nCols!=A.Cols() )
672
    if ( nRows!=A.Rows() || nCols!=A.Cols() )
673
			A.Resize(nRows,nCols);
673
			A.Resize(nRows,nCols);
674
 
674
 
675
		T* RowP;
675
		T* RowP;
676
 
676
 
677
 
677
 
678
    for (int cRow=0; cRow<nRows;cRow++)
678
    for (int cRow=0; cRow<nRows;cRow++)
679
			{
679
			{
680
				RowP=A[cRow];
680
				RowP=A[cRow];
681
        for (int cCol=0;cCol<nCols;cCol++)
681
        for (int cCol=0;cCol<nCols;cCol++)
682
					{
682
					{
683
            s >>  RowP[cCol];
683
            s >>  RowP[cCol];
684
					}
684
					}
685
			}
685
			}
686
 
686
 
687
    return s;
687
    return s;
688
	}
688
	}
689
	//@}
689
	//@}
690
 
690
 
691
	/// The Matrix annotation intended for use.
691
	/// The Matrix annotation intended for use.
692
	typedef CMatrixType<double> CMatrix;
692
	typedef CMatrixType<double> CMatrix;
693
}
693
}
694
#endif // !defined(MATRIX_H_HAA_AGUST_2001)
694
#endif // !defined(MATRIX_H_HAA_AGUST_2001)
695
 
695
 
696
 
696