ROOTPWA
matrix.h
Go to the documentation of this file.
1 // ignore gcc optimization warning:
2 // matrix.h:168:3: warning: assuming signed overflow does not occur when
3 // assuming that (X + c) < X is always false
4 // for more background see http://www.airs.com/blog/archives/120
5 #pragma GCC diagnostic ignored "-Wstrict-overflow"
6 
7 
8 #ifndef __MATRIX_H_
9 #define __MATRIX_H_
10 
11 
12 #include <iostream>
13 #include <sstream>
14 #include <iomanip>
15 #include <limits>
16 #include <cassert>
17 #include <cstring>
18 #include <cmath>
19 #include <complex>
20 
21 #include "Vec.h"
22 
23 
24 inline double conj(const double x) { return x; }
25 
26 
27 template<typename T> class matrix;
28 
29 
30 template<typename T> fourVec operator *= (fourVec& v, const matrix<T>& M);
31 template<typename T> matrix<T> operator * (const T& a, const matrix<T>& M);
32 
33 
34 template<typename T> class matrix {
35 
36 public:
37 
38  matrix();
39  matrix(const int n, const int m);
40  matrix(const int n, const int m, T* p);
41  matrix(const matrix& M);
42  ~matrix();
43  matrix& operator = (const matrix&);
44 
45  T& el (const int i, const int j);
46  T& element(const int i, const int j) { return el(i, j); }
47  const T& el (const int i, const int j) const;
48  const T& element(const int i, const int j) const { return el(i, j); }
49 
50  matrix operator *= (const double a);
51  matrix operator += (const matrix& M);
52  matrix operator -= (const matrix& M);
53  matrix operator * (const matrix& M) const;
54  matrix operator + (const matrix& M) const;
55  matrix operator - (const matrix& M) const;
56  threeVec operator * (const threeVec& v) const;
57  fourVec operator * (const fourVec& v) const;
58  matrix operator * (const T& a) const;
59 
60  matrix conjugate() const;
61  matrix transpose() const;
62  matrix adjoint () const;
63  int status() const { return _status; }
64  void setStatus(const int s) { _status = s; }
65  T trace() const;
66  T det() const;
67  matrix LU() const;
68  matrix inv() ;
69  int nrows() const { return _nrows; }
70  int ncols() const { return _ncols; }
71  matrix zero();
72 
73  const matrix& print(std::ostream& os = std::cout) const;
74  matrix& scan(std::istream& is = std::cin);
75 
76 private:
77 
78  int _nrows;
79  int _ncols;
80  T* _f;
81  int _status;
82 
83  void _create(const int, const int);
84  void _destroy(void);
85  void _copy(const matrix&);
86 
87  matrix _LU(int*, int*) const;
88  matrix _lubksb(int* indx, matrix& b);
89 
90 };
91 
92 
93 template<typename T> class identityMatrix : public matrix<T> {
94 
95 public:
96 
97  identityMatrix(const int n)
98  : matrix<T>(n, n)
99  {
100  for (int i = 0; i < n; ++i)
101  this->el(i, i) = 1;
102  }
103  ~identityMatrix() { }
104 
105 };
106 
107 
108 template<typename T>
109 T
110 matrix<T>::trace() const
111 {
112  T tr = 0;
113  assert(_nrows == _ncols);
114  for (int i = 0; i < _nrows; ++i)
115  tr += el(i, i);
116  return tr;
117 }
118 
119 
120 template<typename T>
121 T
122 matrix<T>::det() const
123 {
124  T dt = 1.0;
125  int* indx = new int[_ncols];
126  int d;
127  assert(_nrows == _ncols);
128  matrix r(_nrows, _ncols);
129  r = _LU(&d, indx);
130  for (int i = 0; i < _nrows; ++i)
131  dt *= r.el(i, i);
132  if (indx)
133  delete[] indx;
134  return d * dt;
135 }
136 
137 
138 template<typename T>
139 double
140 mag(T t)
141 {
142  return abs(std::complex<double>(t));
143 }
144 
145 
146 template<typename T>
147 matrix<T>
148 matrix<T>::_LU(int *d, int *indx) const
149 {
150  int i, j, k, imax = 0;
151  assert(_nrows == _ncols);
152  matrix<T> r(_nrows, _ncols);
153 #define TINY 1.e-20;
154  T big, dum, sum;
155  // T temp;
156  T *vv = new T[_nrows];
157  r = *this;
158  *d = 1;
159  for (i = 0; i < _nrows; ++i) {
160  big = 0.0;
161  for (j = 0; j < _nrows; ++j) {
162  if (mag(r.el(i, j)) > mag(big))
163  big = r.el(i, j);
164  }
165  if (mag(big) == 0.0) {
166  this->print();
167  std::cerr << "singular matrix " << std::endl;
168  r.setStatus(0);
169  return(r);
170  }
171  else
172  r.setStatus(1);
173  vv[i] = 1.0 / big;
174  }
175  for (j = 0; j < _nrows; ++j) {
176  for (i = 0; i < j; ++i) {
177  sum = r.el(i, j);
178  for (k = 0; k < i; ++k)
179  sum -= r.el(i, k) * r.el(k, j);
180  r.el(i, j) = sum;
181  }
182  big = 0.0;
183  for (i = j; i < _nrows; ++i) {
184  sum = r.el(i, j);
185  for (k = 0; k < j; ++k) {
186  sum -= r.el(i, k) * r.el(k, j);
187  }
188  r.el(i, j) = sum;
189  if (mag(vv[i] * sum) >= mag(big)) {
190  big = mag(vv[i] * sum);
191  imax = i;
192  }
193  }
194  if (j != imax) { /* do we need to interchange rows? */
195  for (k = 0; k < _nrows; ++k) { /* Yes, make it so */
196  dum = r.el(imax, k);
197  r.el(imax, k) = r.el(j, k);
198  r.el(j, k) = dum;
199  }
200  *d = -(*d);
201  vv[imax] = vv[j];
202  }
203  indx[j] = imax;
204  if (r.el(j, j) == 0.0)
205  r.el(j, j) = TINY;
206  if (j != _nrows) {
207  dum = 1.0 / r.el(j, j);
208  for (i = j + 1; i < _nrows; ++i)
209  r.el(i, j) *= dum;
210  }
211  }
212  if (vv)
213  delete[] vv;
214  r.setStatus(1);
215  return(r);
216 }
217 
218 
219 template<typename T>
220 matrix<T>
221 matrix<T>::LU() const
222 {
223  int d;
224  int *indx = new int[_nrows];
225  matrix<T> r(_nrows, _ncols);
226  r = _LU(&d, indx);
227  if (indx)
228  delete[] indx;
229  return(r);
230 }
231 
232 
233 template<typename T>
234 matrix<T>
235 matrix<T>::_lubksb(int* indx,
236  matrix& b)
237 {
238  int i, ii = -1, ip, j;
239  T sum;
240  matrix r(b.nrows(), 1);
241  r = b;
242  for (i = 0; i < _nrows; ++i) {
243  ip = indx[i];
244  sum = r.el(ip, 0);
245  r.el(ip, 0) = r.el(i, 0);
246  if (ii > -1) {
247  for (j = ii; j < _nrows; ++j)
248  sum -= element(i, j) * r.element(j, 0);
249  }
250  else if (mag(sum))
251  ii = i;
252  r.el(i, 0) = sum;
253  }
254  for (i = _nrows-1; i > -1; --i) {
255  sum = r.el(i, 0);
256  for (j = i + 1; j < _nrows; ++j)
257  sum -= element(i, j) * r.element(j, 0);
258  r.el(i, 0) = sum / element(i, i);
259  }
260  return(r);
261 }
262 
263 
264 template<typename T>
265 matrix<T>
266 matrix<T>::inv()
267 {
268  int d;
269  int i, j;
270  int *indx = new int[_ncols];
271  matrix r (_nrows, _ncols);
272  matrix lu(_nrows, _ncols);
273  lu = _LU(&d, indx);
274  if (lu.status()) {
275  for (j = 0; j < _nrows; ++j) {
276  matrix col(_nrows, 1);
277  matrix x(_nrows, 1);
278  col.el(j, 0) = 1.0;
279  x = lu._lubksb(indx, col);
280  for (i = 0; i < _nrows; ++i)
281  r.el(i, j) = x.el(i, 0);
282  }
283  r.setStatus(1);
284  if (indx)
285  delete[] indx;
286  } else {
287  r= lu;
288  throw "singular matrix";
289  }
290  return(r);
291 }
292 
293 
294 template<typename T>
295 matrix<T>
296 matrix<T>::conjugate() const
297 {
298  matrix<T> r(_nrows, _ncols);
299  for (int i = 0; i < _nrows; ++i)
300  for (int j = 0; j < _ncols; ++j)
301  r.el(i, j) = conj(el(i, j));
302  return(r);
303 }
304 
305 
306 template<typename T>
307 matrix<T>
308 matrix<T>::transpose() const
309 {
310  matrix<T> r(_ncols, _nrows);
311  for (int i = 0; i < _nrows; ++i)
312  for (int j = 0; j < _ncols; ++j)
313  r.el(j, i) = el(i, j);
314  return(r);
315 }
316 
317 
318 template<typename T>
319 matrix<T>
320 matrix<T>::adjoint() const
321 {
322  matrix<T> r(_ncols, _nrows);
323  r = transpose().conjugate();
324  return(r);
325 }
326 
327 
328 template<typename T>
329 matrix<T>
330 matrix<T>::operator + (const matrix& M) const
331 {
332  assert((_ncols == M._ncols) && (_nrows == M._nrows));
333  matrix<T> r(_nrows, _ncols);
334  for (int i = 0; i < _nrows; ++i)
335  for (int j = 0; j < M._ncols; ++j)
336  r.el(i, j) = el(i, j) + M.el(i, j);
337  return(r);
338 }
339 
340 
341 template<typename T>
342 matrix<T>
343 matrix<T>::operator += (const matrix& M)
344 {
345  assert((_ncols == M._ncols) && (_nrows == M._nrows));
346  for (int i = 0; i < _nrows; ++i)
347  for (int j = 0; j < M._ncols; ++j)
348  el(i, j) += M.el(i, j);
349  return *this;
350 }
351 
352 template<typename T>
353 matrix<T>
354 matrix<T>::operator *= (const double k)
355 {
356  for (int i = 0; i < _nrows; ++i)
357  for (int j = 0; j < _ncols; ++j)
358  el(i, j) *= k;
359  return *this;
360 }
361 
362 
363 template<typename T>
364 matrix<T>
365 matrix<T>::operator - (const matrix& M) const
366 {
367  assert((_ncols == M._ncols) && (_nrows == M._nrows));
368  matrix<T> r(_nrows, _ncols);
369  for (int i = 0; i < _nrows; ++i)
370  for (int j = 0; j < M._ncols; ++j)
371  r.el(i, j) = el(i, j) - M.el(i, j);
372  return(r);
373 }
374 
375 
376 template<typename T>
377 matrix<T>
378 matrix<T>::operator -= (const matrix& M)
379 {
380  assert((_ncols == M._ncols) && (_nrows == M._nrows));
381  for (int i = 0; i < _nrows; ++i)
382  for (int j = 0; j < M._ncols; ++j)
383  el(i, j) -= M.el(i, j);
384  return *this;
385 }
386 
387 
388 template<typename T>
389 matrix<T>
390 matrix<T>::operator * (const matrix& M) const
391 {
392  assert(_ncols == M._nrows);
393  matrix<T> r(_nrows, M._ncols);
394  for (int i = 0; i < _nrows; ++i)
395  for (int j = 0; j < M._ncols; ++j)
396  for (int k = 0; k < _ncols; ++k)
397  r.el(i, j) += el(i, k) * M.el(k, j);
398  return(r);
399 }
400 
401 
402 template<typename T>
403 matrix<T>
404 operator * (const T& a,
405  const matrix<T>& M)
406 {
407  return(M * a);
408 }
409 
410 
411 template<typename T>
412 matrix<T>
413 matrix<T>::operator * (const T& a) const
414 {
415  matrix<T> r(nrows(), ncols());
416  for (int i = 0; i < nrows(); ++i)
417  for (int j = 0; j < ncols(); ++j)
418  r.el(i, j) = a * element(i, j);
419  return(r);
420 }
421 
422 
423 template<typename T>
424 T&
425 matrix<T>::el(const int i,
426  const int j)
427 {
428  assert((i < _nrows) && (j < _ncols));
429  return _f[j + i * _ncols];
430 }
431 
432 
433 template<typename T>
434 const T&
435 matrix<T>::el(const int i,
436  const int j) const
437 {
438  assert((i < _nrows) && (j < _ncols));
439  return _f[j + i * _ncols];
440 }
441 
442 
443 template<typename T>
444 void
445 matrix<T>::_create(const int n,
446  const int m)
447 {
448  _nrows = n;
449  _ncols = m;
450  _f = NULL;
451  if (n * m != 0) {
452  _f = new T[n * m];
453  memset(_f, 0, n * m * sizeof(T));
454  }
455 }
456 
457 
458 template<typename T>
459 void
460 matrix<T>::_copy(const matrix<T>& src)
461 {
462  _nrows = src._nrows;
463  _ncols = src._ncols;
464  _status = src.status();
465  if (_f)
466  delete[] _f;
467  _f = new T[src._nrows * src._ncols];
468  memcpy(_f, src._f, src._nrows * src._ncols * sizeof(T));
469 }
470 
471 
472 template<typename T>
473 void
474 matrix<T>::_destroy(void)
475 {
476  if (_f)
477  delete[] _f;
478  _f = NULL;
479  _nrows = 0;
480  _ncols = 0;
481 }
482 
483 
484 template<typename T>
485 matrix<T>::matrix()
486 {
487  _create(0, 0);
488 }
489 
490 
491 template<typename T>
492 matrix<T>::matrix(const int n,
493  const int m)
494 {
495  _create(n, m);
496 }
497 
498 
499 template<typename T>
500 matrix<T>::matrix(const int n,
501  const int m,
502  T* p)
503 {
504  _create(n, m);
505  _f = p;
506 }
507 
508 
509 template<typename T>
510 matrix<T>::matrix(const matrix<T>& M)
511 {
512  _create(M._nrows, M._ncols);
513  _copy(M);
514 }
515 
516 
517 template<typename T>
518 matrix<T>::~matrix()
519 {
520  _destroy();
521 }
522 
523 
524 template<typename T>
525 matrix<T>&
526 matrix<T>::operator = (const matrix<T>& M)
527 {
528  _copy(M);
529  return *this;
530 }
531 
532 
533 template<typename T>
534 const matrix<T>&
535 matrix<T>::print(std::ostream& os) const
536 {
537  const unsigned int nmbDigits = std::numeric_limits<double>::digits10 + 1;
538  std::ostringstream s;
539  s.precision(nmbDigits);
540  s.setf(std::ios_base::scientific, std::ios_base::floatfield);
541  s << _nrows << " " << _ncols << std::endl;
542  for (int i = 0; i < _nrows; ++i) {
543  for (int j = 0; j < _ncols; ++j)
544  s << _f[j + i * _ncols] << "\t";
545  s << std::endl;
546  }
547  os << s.str();
548  return *this;
549 }
550 
551 
552 template<typename T>
553 matrix<T>&
554 matrix<T>::scan(std::istream& is)
555 {
556  int i = 0, j = 0;
557  is >> i;
558  is >> j;
559  _create(i, j);
560  for (int i = 0; i < _nrows; ++i)
561  for (int j = 0; j < _ncols; ++j) {
562  is >> _f[j + i * _ncols];
563  if (is.eof())
564  throw "trunc matrix input";
565  }
566  return *this;
567 }
568 
569 
570 template<typename T>
571 std::ostream&
572 operator << (std::ostream& os,
573  const matrix<T>& m)
574 {
575  m.print(os);
576  return os;
577 }
578 
579 
580 template<typename T>
581 std::istream&
582 operator >> (std::istream& is,
583  matrix<T>& m)
584 {
585  m.scan(is);
586  return is;
587 }
588 
589 
590 template<typename T>
591 matrix<T>
592 matrix<T>::zero()
593 {
594  memset(_f, 0, _nrows * _ncols * sizeof(T));
595  return *this;
596 }
597 
598 
599 template<typename T>
600 threeVec
601 matrix<T>::operator * (const threeVec& V) const
602 {
603  threeVec R;
604  assert((_nrows == 3) && (_ncols == 3));
605  for (int i = 0; i < 3; ++i)
606  for (int j = 0; j < 3; ++j)
607  R.el(i) += el(i, j) * V.el(j);
608  return(R);
609 }
610 
611 
612 template<typename T>
613 fourVec
614 matrix<T>::operator * (const fourVec& v) const
615 {
616  fourVec r;
617  assert((_nrows == 4) && (_ncols == 4));
618  for (int i = 0; i < 4; ++i)
619  for (int j = 0; j < 4; ++j)
620  r.el(i) += el(i, j) * v.el(j);
621  return(r);
622 }
623 
624 
625 template<typename T>
626 fourVec
627 operator *= (fourVec& v, const matrix<T>& M)
628 {
629  v = M * v;
630  return v;
631 }
632 
633 
634 #endif
Generated by   doxygen 1.8.1.2