1 #ifndef _global_h
2 # include "global.h"
3 #endif
4
5 #include <stdarg.h>
6
7 #ifndef _matrix_h
8 # include "matrix.h"
9 #endif
10
11 Matrix MatrixZero;
12
Matrix(int n,int m,...)13 Matrix::Matrix( int n, int m, ... )
14 {
15 va_list argptr;
16
17 va_start(argptr,m);
18 dim = n;
19 if (dim) {
20 data = new Vector [n];
21 for (int i=0;i<n;i++) {
22 data[i].Resize(m);
23 for (int j=0;j<m;j++) {
24 data[i][j] = Real( va_arg(argptr,double) );
25 }
26 }
27 }
28 va_end(argptr);
29 }
30
Matrix(const Matrix & m)31 Matrix::Matrix( const Matrix &m ) {
32 dim = m.dim;
33 data = new Vector[dim];
34 for (int i=0;i<dim;i++) data[i] = m.data[i];
35 }
36
show(const char * str)37 void Matrix::show(const char *str) const {
38 putchar( '[' );
39 for (int i=0;i<dim;i++) {
40 if (i) printf( ", " );
41 data[i].show();
42 }
43 putchar(']');
44 if (str) printf(str);
45 }
46
Cols()47 int Matrix::Cols() const {
48 int max=0;
49
50 for (int i=0;i<dim;i++) if (data[i].dim>max) max=data[i].dim;
51 return max;
52 }
53
IsZero()54 int Matrix::IsZero() const {
55 for (int i=0;i<dim;i++) {
56 if (data[i]!=VectorZero) return 0;
57 }
58 return -1;
59 }
60
Resize(int n)61 void Matrix::Resize(int n) {
62 if (n>dim) {
63 Vector *new_data;
64
65 new_data = new Vector[n]; // allocate new array
66 for (int i=0;i<dim;i++) new_data[i] = data[i]; // copy old data
67 // for ( ;i<n;i++) new_data[i] = VectorZero; // reset new data
68 delete [] data; // free old array
69 data = new_data;
70 dim = n;
71 }
72 }
Resize(int n,int m)73 void Matrix::Resize(int n,int m) {
74 Resize(n);
75 for (int i=0;i<dim;i++) data[i].Resize(m);
76 }
77
78
GetRef(int n)79 Vector &Matrix::GetRef(int n) {
80 static Vector error=VectorZero;
81
82 if (n<0) {
83 fprintf(stderr,"illegal index: Matrix::GetRef(%d)\n",n);
84 return error; // return at least something to play with ...
85 }
86 if (n>=dim) Resize(n+1);
87 return data[n];
88 }
89
90
91 const Matrix &Matrix::operator=(const Matrix &m) {
92 if (dim<m.dim) {
93 if (dim) delete [] data;
94 data = new Vector[m.dim];
95 }
96 dim = m.dim;
97 for (int i=0;i<dim;i++) data[i] = m.data[i];
98 return *this;
99 }
100
101 int operator==(const Matrix& z1, const Matrix& z2)
102 {
103 if (z1.dim<z2.dim) {
104 for (int i=0;i<z1.dim;i++)
105 if (z1.data[i]!=z2.data[i]) return 0;
106 for ( ;i<z2.dim;i++)
107 if (z2.data[i]!=VectorZero) return 0;
108 return -1;
109 }
110 else {
111 for (int i=0;i<z2.dim;i++)
112 if (z1.data[i]!=z2.data[i]) return 0;
113 for ( ;i<z1.dim;i++)
114 if (z1.data[i]!=VectorZero) return 0;
115 return -1;
116 }
117 }
118
119 const Matrix& Matrix::operator+=(const Matrix& z2)
120 {
121 if (z2.dim>dim) Resize(z2.dim);
122 for (int i=0;i<dim;i++) data[i]+=z2.data[i];
123 return *this;
124 }
125 const Matrix& Matrix::operator-=(const Matrix& z2)
126 {
127 if (z2.dim>dim) Resize(z2.dim);
128 for (int i=0;i<dim;i++) data[i]-=z2.data[i];
129 return *this;
130 }
131 const Matrix& Matrix::operator*=(const Real& r)
132 {
133 for (int i=0;i<dim;i++) data[i]*=r;
134 return *this;
135 }
136 const Matrix& Matrix::operator/=(const Real& r)
137 {
138 for (int i=0;i<dim;i++) data[i]/=r;
139 return *this;
140 }
141
142 Vector operator*(const Vector& v, const Matrix& m) {
143 Vector erg;
144
145 erg.Resize(m.Cols());
146 for (int i=0;i<erg.dim;i++) {
147 Real sum = RealZero;
148 for (int j=0;j<v.dim;j++) sum += v(j)*m(j)(i);
149 erg[i] = sum;
150 }
151 return erg;
152 }
153
154 Vector operator*(const Matrix& m, const Vector& v) {
155 Vector erg;
156
157 erg.Resize(m.dim);
158 for (int i=0;i<m.dim;i++) erg[i] = m(i)*v;
159 return erg;
160 }
161
162 const Matrix &Matrix::operator<<=(Matrix &m2) {
163 if (dim) delete [] data;
164 dim = m2.dim;
165 data = m2.data;
166 m2.dim = 0;
167 return *this;
168 }
169
170 const Matrix &Matrix::operator*=(const Matrix &m2) {
171 Matrix help;
172
173 help <<= *this;
174
175 int cols = m2.Cols();
176
177 Resize(help.dim,cols);
178 for (int i=0;i<dim;i++) {
179 for (int j=0;j<cols;j++) {
180 Real sum = RealZero;
181 for (int k=0;k<m2.dim;k++) sum += help(i)(k)*m2(k)(j);
182 (*this)[i][j] = sum;
183 }
184 }
185 return *this;
186 }
187
188 Matrix operator*(const Matrix& m1, const Matrix& m2) {
189 Matrix erg;
190 int cols = m2.Cols();
191
192 erg.Resize(m1.dim,cols);
193 for (int i=0;i<erg.dim;i++) {
194 for (int j=0;j<cols;j++) {
195 Real sum = RealZero;
196 for (int k=0;k<m2.dim;k++) sum += m1(i)(k)*m2(k)(j);
197 erg[i][j] = sum;
198 }
199 }
200 return erg;
201 }
202
203
204 Vector operator/(const Vector &b, const Matrix &a) {
205 Matrix a_tmp = a;
206 Vector b_tmp = b;
207 Vector x;
208
209 if (system_calc(a_tmp,&x,b_tmp)) return VectorZero;
210 else return x;
211 }
212
system_calc(Matrix & a,Vector * x,Vector & b)213 int system_calc(Matrix &a, Vector *x, Vector &b) {
214 int n; // Dimension
215 int *v; // Vertauschungsvektor
216 int i,j,k,help;
217
218 // �berpr�fung der Vektorgr��en
219
220 n = b.dim;
221 if (n!=a.Rows()||n!=a.Cols()) {
222 fprintf(stderr,"Ax=b: Fehler in Dimension A[%d][%d], b[%d]\n", a.Rows(), a.Cols(), n );
223 return -1;
224 }
225
226 // Initialisierung des Vertauschungsvektors
227
228 v = new int[n];
229 for (i=0;i<n;i++) v[i]=i;
230
231
232 for (j=0;j<n-1;j++) {
233 Real pvt = RealZero;
234 int j2 = j;
235
236 // skalierte Pivotsuche
237
238 for (i=j;i<n;i++) {
239 Real a2;
240 Real d = RealZero;
241
242 // Zeilenmaximum f�r Skalierung bestimmen
243
244 for (k=j;k<n;k++) {
245 if (Real(fabs(a(v[i])(k)))>d) d=fabs(a(v[i])(k));
246 }
247 if (d<EPS) {
248 fprintf(stderr,"Ax=b: Matrix singul�r\n");
249 delete v;
250 return -2;
251 }
252
253 // skaliertes Pivotelement berechnen und mit bisher gr��tem vergleichen
254
255 a2 = a(v[i])(j)/d;
256 if (Real(fabs(a2))>pvt) {
257 pvt = fabs(a2);
258 j2 = i; // Pivotzeile merken
259 }
260 }
261
262 // Pivotisierung und fiktiever Zeilentausch
263
264 help = v[j];
265 v[j] = v[j2];
266 v[j2] = help;
267 if (fabs(a(v[j])(j))<EPS) {
268 fprintf(stderr,"Ax=b: Matrix singul�r\n" );
269 delete v;
270 return -3;
271 }
272
273 // v[j] ist Pivotzeile
274
275 for (i=j+1;i<n;i++) {
276 Real q = a(v[i])(j)/a(v[j])(j);
277 for (k=j+1;k<n;k++) {
278 a[v[i]][k] -= q*a(v[j])(k);
279 }
280 b[v[i]] -= q*b(v[j]);
281 }
282 }
283
284 x->Resize(n);
285 for (i=n-1;i>=0;i--) {
286 (*x)[v[i]] = b(v[i]);
287 for (k=i+1;k<n;k++) {
288 (*x)[v[i]] -= a(v[i])(k) * (*x)(v[k]);
289 }
290 (*x)[v[i]] /= a(v[i])(i);
291 }
292
293 delete v;
294 return 0;
295 }
296
297