1 #ifndef FILE_DENSEMAT
2 #define FILE_DENSEMAT
3
4 /**************************************************************************/
5 /* File: densemat.hpp */
6 /* Author: Joachim Schoeberl */
7 /* Date: 01. Oct. 94 */
8 /**************************************************************************/
9
10 /**
11 Data type dense matrix
12 */
13
14
15 class DenseMatrix
16 {
17 protected:
18 int height;
19 int width;
20 double * data;
21
22 public:
23 ///
24 DLL_HEADER DenseMatrix ();
25 ///
26 DLL_HEADER DenseMatrix (int h, int w = 0);
27 ///
28 DLL_HEADER DenseMatrix (const DenseMatrix & m2);
29 ///
30 DLL_HEADER ~DenseMatrix ();
31
32 ///
33 DLL_HEADER void SetSize (int h, int w = 0);
34
Height() const35 int Height() const { return height; }
Width() const36 int Width() const {return width; }
37
operator ()(int i,int j)38 double & operator() (int i, int j) { return data[i*width+j]; }
operator ()(int i,int j) const39 double operator() (int i, int j) const { return data[i*width+j]; }
operator ()(int i)40 double & operator() (int i) { return data[i]; }
operator ()(int i) const41 double operator() (int i) const { return data[i]; }
42
43 ///
44 DLL_HEADER DenseMatrix & operator= (const DenseMatrix & m2);
45 ///
46 DLL_HEADER DenseMatrix & operator+= (const DenseMatrix & m2);
47 ///
48 DLL_HEADER DenseMatrix & operator-= (const DenseMatrix & m2);
49
50 ///
51 DLL_HEADER DenseMatrix & operator= (double v);
52 ///
53 DLL_HEADER DenseMatrix & operator*= (double v);
54
55 ///
Mult(const FlatVector & v,FlatVector & prod) const56 DLL_HEADER void Mult (const FlatVector & v, FlatVector & prod) const
57 {
58 double sum;
59 const double * mp, * sp;
60 double * dp;
61
62 #ifdef DEBUG
63 if (prod.Size() != height)
64 {
65 (*myerr) << "Mult: wrong vector size " << endl;
66 }
67 if (!height)
68 {
69 cout << "DenseMatrix::Mult height = 0" << endl;
70 }
71 if (!width)
72 {
73 cout << "DenseMatrix::Mult width = 0" << endl;
74 }
75
76 if (width != v.Size())
77 {
78 (*myerr) << "\nMatrix and Vector don't fit" << endl;
79 }
80 else if (Height() != prod.Size())
81 {
82 (*myerr) << "Base_Matrix::operator*(Vector): prod vector not ok" << endl;
83 }
84 else
85 #endif
86 {
87 mp = data;
88 dp = &prod(0);
89 for (int i = 0; i < height; i++)
90 {
91 sum = 0;
92 sp = &v(0);
93
94 for (int j = 0; j < width; j++)
95 {
96 // sum += Get(i,j) * v.Get(j);
97 sum += *mp * *sp;
98 mp++;
99 sp++;
100 }
101
102 *dp = sum;
103 dp++;
104 }
105 }
106 }
107
108 ///
109 DLL_HEADER void MultTrans (const Vector & v, Vector & prod) const;
110 ///
111 DLL_HEADER void Residuum (const Vector & x, const Vector & b, Vector & res) const;
112 ///
113 DLL_HEADER double Det () const;
114
115 ///
116 friend DenseMatrix operator* (const DenseMatrix & m1, const DenseMatrix & m2);
117 ///
118 friend DenseMatrix operator+ (const DenseMatrix & m1, const DenseMatrix & m2);
119
120 ///
121 friend void Transpose (const DenseMatrix & m1, DenseMatrix & m2);
122 ///
123 friend void Mult (const DenseMatrix & m1, const DenseMatrix & m2, DenseMatrix & m3);
124 ///
125 // friend void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2);
126 ///
127 friend void CalcAAt (const DenseMatrix & a, DenseMatrix & m2);
128 ///
129 // friend void CalcAtA (const DenseMatrix & a, DenseMatrix & m2);
130 ///
131 friend void CalcABt (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2);
132 ///
133 friend void CalcAtB (const DenseMatrix & a, const DenseMatrix & b, DenseMatrix & m2);
134 ///
135 DLL_HEADER void Solve (const Vector & b, Vector & x) const;
136 ///
137 void SolveDestroy (const Vector & b, Vector & x);
138 ///
Get(int i,int j) const139 const double & Get(int i, int j) const { return data[(i-1)*width+j-1]; }
140 ///
Get(int i) const141 const double & Get(int i) const { return data[i-1]; }
142 ///
Set(int i,int j,double v)143 void Set(int i, int j, double v) { data[(i-1)*width+j-1] = v; }
144 ///
Elem(int i,int j)145 double & Elem(int i, int j) { return data[(i-1)*width+j-1]; }
146 ///
ConstElem(int i,int j) const147 const double & ConstElem(int i, int j) const { return data[(i-1)*width+j-1]; }
148 };
149
150
151 extern ostream & operator<< (ostream & ost, const DenseMatrix & m);
152
153
154
155 template <int WIDTH, typename T = double>
156 class MatrixFixWidth
157 {
158 protected:
159 int height;
160 T * __restrict data;
161 bool ownmem;
162 public:
163 ///
MatrixFixWidth()164 MatrixFixWidth ()
165 { height = 0; data = 0; ownmem = false; }
166 ///
MatrixFixWidth(int h)167 MatrixFixWidth (int h)
168 { height = h; data = new T[WIDTH*height]; ownmem = true; }
169 ///
MatrixFixWidth(int h,T * adata)170 MatrixFixWidth (int h, T * adata)
171 { height = h; data = adata; ownmem = false; }
172 ///
MatrixFixWidth(const MatrixFixWidth & m2)173 MatrixFixWidth (const MatrixFixWidth & m2)
174 : height(m2.height), ownmem(true)
175 {
176 data = new T[height*WIDTH];
177 for (int i = 0; i < WIDTH*height; i++)
178 data[i] = m2.data[i];
179 }
180 // : height(m2.height), data(m2.data), ownmem(false)
181 //{ ; }
182 ///
~MatrixFixWidth()183 ~MatrixFixWidth ()
184 { if (ownmem) delete [] data; }
185
SetSize(int h)186 void SetSize (int h)
187 {
188 if (h != height)
189 {
190 if (ownmem) delete data;
191 height = h;
192 data = new T[WIDTH*height];
193 ownmem = true;
194 }
195 }
196
197 ///
Height() const198 int Height() const { return height; }
199
200 ///
Width() const201 int Width() const { return WIDTH; }
operator =(const MatrixFixWidth & m2)202 MatrixFixWidth & operator= (const MatrixFixWidth & m2)
203 {
204 for (int i = 0; i < height*WIDTH; i++)
205 data[i] = m2.data[i];
206 }
207 ///
operator =(T v)208 MatrixFixWidth & operator= (T v)
209 {
210 for (int i = 0; i < height*WIDTH; i++)
211 data[i] = v;
212 return *this;
213 }
214
215 /*
216 ///
217 void Mult (const FlatVector & v, FlatVector & prod) const
218 {
219 T sum;
220 const T * mp, * sp;
221 T * dp;
222
223 mp = data;
224 dp = &prod[0];
225 for (int i = 0; i < height; i++)
226 {
227 sum = 0;
228 sp = &v[0];
229
230 for (int j = 0; j < WIDTH; j++)
231 {
232 sum += *mp * *sp;
233 mp++;
234 sp++;
235 }
236
237 *dp = sum;
238 dp++;
239 }
240 }
241 */
242
operator ()(int i,int j)243 T & operator() (int i, int j)
244 { return data[i*WIDTH+j]; }
245
operator ()(int i,int j) const246 const T & operator() (int i, int j) const
247 { return data[i*WIDTH+j]; }
248
249
operator *=(T v)250 MatrixFixWidth & operator*= (T v)
251 {
252 if (data)
253 for (int i = 0; i < height*WIDTH; i++)
254 data[i] *= v;
255 return *this;
256 }
257
258
259
Get(int i,int j) const260 const T & Get(int i, int j) const { return data[(i-1)*WIDTH+j-1]; }
261 ///
Get(int i) const262 const T & Get(int i) const { return data[i-1]; }
263 ///
Set(int i,int j,T v)264 void Set(int i, int j, T v) { data[(i-1)*WIDTH+j-1] = v; }
265 ///
Elem(int i,int j)266 T & Elem(int i, int j) { return data[(i-1)*WIDTH+j-1]; }
267 ///
ConstElem(int i,int j) const268 const T & ConstElem(int i, int j) const { return data[(i-1)*WIDTH+j-1]; }
269 };
270
271
272 template <int WIDTH>
273 class MatrixFixWidth<WIDTH,double>
274 {
275 protected:
276 int height;
277 double * data;
278 bool ownmem;
279 public:
280 ///
MatrixFixWidth()281 MatrixFixWidth ()
282 { height = 0; data = 0; ownmem = false; }
283 ///
MatrixFixWidth(int h)284 MatrixFixWidth (int h)
285 { height = h; data = new double[WIDTH*height]; ownmem = true; }
286
MatrixFixWidth(const MatrixFixWidth & m2)287 MatrixFixWidth (const MatrixFixWidth & m2)
288 : height(m2.height), ownmem(true)
289 {
290 data = new double[height*WIDTH];
291 for (int i = 0; i < WIDTH*height; i++)
292 data[i] = m2.data[i];
293 }
294
295 ///
MatrixFixWidth(int h,double * adata)296 MatrixFixWidth (int h, double * adata)
297 { height = h; data = adata; ownmem = false; }
298 ///
~MatrixFixWidth()299 ~MatrixFixWidth ()
300 { if (ownmem) delete [] data; }
301
SetSize(int h)302 void SetSize (int h)
303 {
304 if (h != height)
305 {
306 if (ownmem) delete data;
307 height = h;
308 data = new double[WIDTH*height];
309 ownmem = true;
310 }
311 }
312
313 ///
Height() const314 int Height() const { return height; }
315
316 ///
Width() const317 int Width() const { return WIDTH; }
318
319 ///
operator =(double v)320 MatrixFixWidth & operator= (double v)
321 {
322 for (int i = 0; i < height*WIDTH; i++)
323 data[i] = v;
324 return *this;
325 }
326
327 ///
Mult(const FlatVector & v,FlatVector & prod) const328 void Mult (const FlatVector & v, FlatVector & prod) const
329 {
330 double sum;
331 const double * mp, * sp;
332 double * dp;
333
334 /*
335 if (prod.Size() != height)
336 {
337 cerr << "MatrixFixWidth::Mult: wrong vector size " << endl;
338 assert (1);
339 }
340 */
341
342 mp = data;
343 dp = &prod[0];
344 for (int i = 0; i < height; i++)
345 {
346 sum = 0;
347 sp = &v[0];
348
349 for (int j = 0; j < WIDTH; j++)
350 {
351 sum += *mp * *sp;
352 mp++;
353 sp++;
354 }
355
356 *dp = sum;
357 dp++;
358 }
359 }
360
operator ()(int i,int j)361 double & operator() (int i, int j)
362 { return data[i*WIDTH+j]; }
363
operator ()(int i,int j) const364 const double & operator() (int i, int j) const
365 { return data[i*WIDTH+j]; }
366
367
operator *=(double v)368 MatrixFixWidth & operator*= (double v)
369 {
370 if (data)
371 for (int i = 0; i < height*WIDTH; i++)
372 data[i] *= v;
373 return *this;
374 }
375
376
377
Get(int i,int j) const378 const double & Get(int i, int j) const { return data[(i-1)*WIDTH+j-1]; }
379 ///
Get(int i) const380 const double & Get(int i) const { return data[i-1]; }
381 ///
Set(int i,int j,double v)382 void Set(int i, int j, double v) { data[(i-1)*WIDTH+j-1] = v; }
383 ///
Elem(int i,int j)384 double & Elem(int i, int j) { return data[(i-1)*WIDTH+j-1]; }
385 ///
ConstElem(int i,int j) const386 const double & ConstElem(int i, int j) const { return data[(i-1)*WIDTH+j-1]; }
387 };
388
389
390
391
392
393 template <int WIDTH>
operator <<(ostream & ost,const MatrixFixWidth<WIDTH> & m)394 extern ostream & operator<< (ostream & ost, const MatrixFixWidth<WIDTH> & m)
395 {
396 for (int i = 0; i < m.Height(); i++)
397 {
398 for (int j = 0; j < m.Width(); j++)
399 ost << m.Get(i+1,j+1) << " ";
400 ost << endl;
401 }
402 return ost;
403 };
404
405
406 extern DLL_HEADER void CalcAtA (const DenseMatrix & a, DenseMatrix & m2);
407 extern DLL_HEADER void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2);
408
409
410 #endif
411