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>
156 class MatrixFixWidth
157 {
158 protected:
159 int height;
160 double * data;
161
162 public:
163 ///
MatrixFixWidth()164 MatrixFixWidth ()
165 { height = 0; data = 0; }
166 ///
MatrixFixWidth(int h)167 MatrixFixWidth (int h)
168 { height = h; data = new double[WIDTH*height]; }
169 ///
~MatrixFixWidth()170 ~MatrixFixWidth ()
171 { delete [] data; }
172
SetSize(int h)173 void SetSize (int h)
174 {
175 if (h != height)
176 {
177 delete data;
178 height = h;
179 data = new double[WIDTH*height];
180 }
181 }
182
183 ///
Height() const184 int Height() const { return height; }
185
186 ///
Width() const187 int Width() const { return WIDTH; }
188
189 ///
operator =(double v)190 MatrixFixWidth & operator= (double v)
191 {
192 for (int i = 0; i < height*WIDTH; i++)
193 data[i] = v;
194 return *this;
195 }
196
197 ///
Mult(const FlatVector & v,FlatVector & prod) const198 void Mult (const FlatVector & v, FlatVector & prod) const
199 {
200 double sum;
201 const double * mp, * sp;
202 double * dp;
203
204 /*
205 if (prod.Size() != height)
206 {
207 cerr << "MatrixFixWidth::Mult: wrong vector size " << endl;
208 assert (1);
209 }
210 */
211
212 mp = data;
213 dp = &prod[0];
214 for (int i = 0; i < height; i++)
215 {
216 sum = 0;
217 sp = &v[0];
218
219 for (int j = 0; j < WIDTH; j++)
220 {
221 sum += *mp * *sp;
222 mp++;
223 sp++;
224 }
225
226 *dp = sum;
227 dp++;
228 }
229 }
230
operator ()(int i,int j)231 double & operator() (int i, int j)
232 { return data[i*WIDTH+j]; }
233
operator ()(int i,int j) const234 const double & operator() (int i, int j) const
235 { return data[i*WIDTH+j]; }
236
237
operator *=(double v)238 MatrixFixWidth & operator*= (double v)
239 {
240 if (data)
241 for (int i = 0; i < height*WIDTH; i++)
242 data[i] *= v;
243 return *this;
244 }
245
246
247
Get(int i,int j) const248 const double & Get(int i, int j) const { return data[(i-1)*WIDTH+j-1]; }
249 ///
Get(int i) const250 const double & Get(int i) const { return data[i-1]; }
251 ///
Set(int i,int j,double v)252 void Set(int i, int j, double v) { data[(i-1)*WIDTH+j-1] = v; }
253 ///
Elem(int i,int j)254 double & Elem(int i, int j) { return data[(i-1)*WIDTH+j-1]; }
255 ///
ConstElem(int i,int j) const256 const double & ConstElem(int i, int j) const { return data[(i-1)*WIDTH+j-1]; }
257 };
258
259
260 template <int WIDTH>
operator <<(ostream & ost,const MatrixFixWidth<WIDTH> & m)261 extern ostream & operator<< (ostream & ost, const MatrixFixWidth<WIDTH> & m)
262 {
263 for (int i = 0; i < m.Height(); i++)
264 {
265 for (int j = 0; j < m.Width(); j++)
266 ost << m.Get(i+1,j+1) << " ";
267 ost << endl;
268 }
269 return ost;
270 };
271
272
273 extern DLL_HEADER void CalcAtA (const DenseMatrix & a, DenseMatrix & m2);
274 extern DLL_HEADER void CalcInverse (const DenseMatrix & m1, DenseMatrix & m2);
275
276
277 #endif
278