Jpp  15.0.0-rc.2
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
JMatrixND.hh
Go to the documentation of this file.
1 #ifndef __JMATH__JMATRIXND__
2 #define __JMATH__JMATRIXND__
3 
4 #include <cstring>
5 #include <limits>
6 #include <cmath>
7 #include <algorithm>
8 
9 #include "JIO/JSerialisable.hh"
10 #include "JLang/JException.hh"
11 #include "JLang/JEquals.hh"
12 #include "JLang/JSingleton.hh"
13 #include "JLang/JManip.hh"
14 #include "JMath/JMath.hh"
15 #include "JMath/JVectorND.hh"
16 
17 /**
18  * \author mdejong
19  */
20 
21 namespace JMATH {}
22 namespace JPP { using namespace JMATH; }
23 
24 namespace JMATH {
25 
26  using JIO::JReader;
27  using JIO::JWriter;
28  using JLANG::JEquals;
29  using JLANG::JSingleton;
30  using JLANG::JException;
33 
34 
35  /**
36  * Basic NxN matrix.
37  */
38  struct JMatrixND_t :
39  public JSingleton<JMatrixND_t>
40  {
41  /**
42  * Default constructor.
43  */
45  __p(NULL),
46  __n(0),
47  __m(0)
48  {}
49 
50 
51  /**
52  * Destructor.
53  */
55  {
56  clear();
57  }
58 
59 
60  /**
61  * Clear memory.
62  */
63  void clear()
64  {
65  if (__p != NULL) {
66  delete [] __p;
67  }
68 
69  __p = NULL;
70  __n = 0;
71  __m = 0;
72  }
73 
74 
75  /**
76  * Resize matrix.
77  *
78  * Note that this method does not maintain data in the matrix.
79  *
80  * \param size dimension
81  */
82  void resize(const size_t size)
83  {
84  if (size != this->__n) {
85 
86  if (size > this->__m) {
87 
88  if (__p != NULL) {
89  delete[] __p;
90  }
91 
92  __m = size;
93  __p = new double[__m*__m];
94 
95  if (__p == NULL) {
96  THROW(JNewException, "JMatrixND::resize(" << size << "): Memory allocation failure.");
97  }
98  }
99 
100  __n = size;
101  }
102  }
103 
104 
105  /**
106  * Set matrix.
107  *
108  * \param A matrix
109  */
110  void set(const JMatrixND_t& A)
111  {
112  this->resize(A.size());
113 
114  memcpy(this->data(), A.data(), A.size() * A.size() * sizeof(double));
115  }
116 
117 
118  /**
119  * Swap matrices
120  *
121  * \param A matrix
122  */
124  {
125  using std::swap;
126 
127  swap(this->__p, A.__p);
128  swap(this->__n, A.__n);
129  swap(this->__m, A.__m);
130  }
131 
132 
133  /**
134  * Transpose.
135  *
136  * \return this matrix
137  */
139  {
140  using std::swap;
141 
142  for (size_t row = 0; row != this->size(); ++row) {
143  for (size_t col = 0; col != row; ++col) {
144  swap((*this)(row,col), (*this)(col,row));
145  }
146  }
147 
148  return *this;
149  }
150 
151 
152  /**
153  * Get dimension of matrix.
154  *
155  * \return dimension
156  */
157  size_t size() const
158  {
159  return __n;
160  }
161 
162 
163  /**
164  * Get capacity of dimension.
165  *
166  * \return capacity
167  */
168  size_t capacity() const
169  {
170  return __m;
171  }
172 
173 
174  /**
175  * Check emptyness.
176  *
177  * \return true if empty; else false
178  */
179  bool empty() const
180  {
181  return __n == 0;
182  }
183 
184 
185  /**
186  * Get pointer to data.
187  *
188  * \return pointer to data.
189  */
190  const double* data() const
191  {
192  return __p;
193  }
194 
195 
196  /**
197  * Get pointer to data.
198  *
199  * \return pointer to data.
200  */
201  double* data()
202  {
203  return __p;
204  }
205 
206 
207  /**
208  * Get row data.
209  *
210  * \param row row number
211  * \return pointer to row data
212  */
213  const double* operator[](size_t row) const
214  {
215  return __p + row*__n;
216  }
217 
218 
219  /**
220  * Get row data.
221  *
222  * \param row row number
223  * \return pointer to row data
224  */
225  double* operator[](size_t row)
226  {
227  return __p + row*__n;
228  }
229 
230 
231  /**
232  * Get matrix element.
233  *
234  * \param row row number
235  * \param col column number
236  * \return matrix element at (row,col)
237  */
238  double operator()(const size_t row, const size_t col) const
239  {
240  return *(__p + row*__n + col);
241  }
242 
243 
244  /**
245  * Get matrix element.
246  *
247  * \param row row number
248  * \param col column number
249  * \return matrix element at (row,col)
250  */
251  double& operator()(const size_t row, const size_t col)
252  {
253  return *(__p + row*__n + col);
254  }
255 
256 
257  /**
258  * Get matrix element.
259  *
260  * \param row row number
261  * \param col column number
262  * \return matrix element at (row,col)
263  */
264  double at(size_t row, size_t col) const
265  {
266  if (row >= this->size() ||
267  col >= this->size()) {
268  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
269  }
270 
271  return (*this)(row, col);
272  }
273 
274 
275  /**
276  * Get matrix element.
277  *
278  * \param row row number
279  * \param col column number
280  * \return matrix element at (row,col)
281  */
282  double& at(size_t row, size_t col)
283  {
284  if (row >= this->size() ||
285  col >= this->size()) {
286  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
287  }
288 
289  return (*this)(row, col);
290  }
291 
292 
293  protected:
294  double* __p; //!< pointer to data
295  size_t __n; //!< dimension of matrix
296  size_t __m; //!< capacity of matrix
297  };
298 
299 
300  /**
301  * NxN matrix.
302  */
303  struct JMatrixND :
304  public JMatrixND_t,
305  public JMath <JMatrixND>,
306  public JEquals<JMatrixND>
307  {
308  using JMath<JMatrixND>::mul;
309 
310  /**
311  * Default constructor.
312  */
314  {}
315 
316 
317  /**
318  * Constructor.
319  *
320  * The matrix is set to zero.
321  *
322  * \param size dimension
323  */
324  JMatrixND(const size_t size)
325  {
326  resize(size);
327  reset();
328  }
329 
330 
331  /**
332  * Copy constructor.
333  *
334  * \param A matrix
335  */
337  {
338  set(A);
339  }
340 
341 
342  /**
343  * Assignment operator.
344  *
345  * \param A matrix
346  */
348  {
349  this->set(A);
350 
351  return *this;
352  }
353 
354 
355  /**
356  * Resize matrix.
357  *
358  * Note that this method does not maintain data in the matrix.
359  *
360  * \param size dimension
361  */
362  void resize(const size_t size)
363  {
364  static_cast<JMatrixND_t&>(*this).resize(size);
365 
366  getInstance().resize(size);
367  }
368 
369 
370  /**
371  * Set matrix to the null matrix.
372  *
373  * \return this matrix
374  */
376  {
378 
379  double* p = A.data();
380 
381  for (size_t i = this->size(); i != 0; --i, ++p) {
382  *p = 0.0;
383  }
384 
385  p = this->data();
386 
387  for (size_t i = this->size(); i != 0; --i, p += this->size()) {
388  memcpy(p, A.data(), this->size() * sizeof(double));
389  }
390 
391  return *this;
392  }
393 
394 
395  /**
396  * Set to identity matrix
397  *
398  * \return this matrix
399  */
401  {
402  reset();
403 
404  for (size_t i = 0; i != this->size(); ++i) {
405  (*this)(i,i) = 1.0;
406  }
407 
408  return *this;
409  }
410 
411 
412  /**
413  * Negate matrix.
414  *
415  * \return this matrix
416  */
418  {
419  double* p = this->data();
420 
421  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
422  *p = -*p;
423  }
424 
425  return *this;
426  }
427 
428 
429  /**
430  * Matrix addition.
431  *
432  * \param A matrix
433  * \return this matrix
434  */
436  {
437  if (this->size() == A.size()) {
438 
439  double* p = this->data();
440  const double* q = A.data();
441 
442  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
443  *p += *q;
444  }
445 
446  } else {
447  THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
448  }
449  }
450 
451 
452  /**
453  * Matrix subtraction.
454  *
455  * \param A matrix
456  * \return this matrix
457  */
459  {
460  if (this->size() == A.size()) {
461 
462  double* p = this->data();
463  const double* q = A.data();
464 
465  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
466  *p -= *q;
467  }
468 
469  } else {
470  THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
471  }
472 
473  return *this;
474  }
475 
476 
477  /**
478  * Scale matrix.
479  *
480  * \param factor factor
481  * \return this matrix
482  */
483  JMatrixND& mul(const double factor)
484  {
485  double* p = this->data();
486 
487  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
488  *p *= factor;
489  }
490 
491  return *this;
492  }
493 
494 
495  /**
496  * Scale matrix.
497  *
498  * \param factor factor
499  * \return this matrix
500  */
501  JMatrixND& div(const double factor)
502  {
503  double* p = this->data();
504 
505  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
506  *p /= factor;
507  }
508 
509  return *this;
510  }
511 
512 
513  /**
514  * Matrix multiplication.
515  *
516  * \param A matrix
517  * \param B matrix
518  * \return this matrix
519  */
521  const JMatrixND& B)
522  {
523  if (A.size() == B.size()) {
524 
525  this->resize(A.size());
526 
527  if (!this->empty()) {
528 
530 
531  C.set(B);
532  C.transpose();
533 
534  size_t i, row;
535 
536  for (row = 0; row + 4 <= this->size(); row += 4) { // process rows by four
537 
538  double* p0 = (*this)[row + 0];
539  double* p1 = (*this)[row + 1];
540  double* p2 = (*this)[row + 2];
541  double* p3 = (*this)[row + 3];
542 
543  for (size_t col = 0; col != this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
544 
545  double w0 = 0.0;
546  double w1 = 0.0;
547  double w2 = 0.0;
548  double w3 = 0.0;
549 
550  const double* a0 = A[row + 0];
551  const double* a1 = A[row + 1];
552  const double* a2 = A[row + 2];
553  const double* a3 = A[row + 3];
554 
555  const double* c0 = C[col];
556 
557  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4, c0 += 4) {
558  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
559  w1 += a1[0] * c0[0] + a1[1] * c0[1] + a1[2] * c0[2] + a1[3] * c0[3];
560  w2 += a2[0] * c0[0] + a2[1] * c0[1] + a2[2] * c0[2] + a2[3] * c0[3];
561  w3 += a3[0] * c0[0] + a3[1] * c0[1] + a3[2] * c0[2] + a3[3] * c0[3];
562  }
563 
564  for ( ; i != this->size(); ++i, ++a0, ++a1, ++a2, ++a3, ++c0) {
565  w0 += (*a0) * (*c0);
566  w1 += (*a1) * (*c0);
567  w2 += (*a2) * (*c0);
568  w3 += (*a3) * (*c0);
569  }
570 
571  *p0 = w0;
572  *p1 = w1;
573  *p2 = w2;
574  *p3 = w3;
575  }
576  }
577 
578  for ( ; row != this->size(); ++row) { // remaining rows
579 
580  double* p0 = (*this)[row + 0];
581 
582  for (size_t col = 0; col != this->size(); ++col, ++p0) {
583 
584  double w0 = 0.0;
585 
586  const double* a0 = A[row + 0];
587  const double* c0 = C[col];
588 
589  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, c0 += 4) {
590  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
591  }
592 
593  for ( ; i != this->size(); ++i, ++a0, ++c0) {
594  w0 += (*a0) * (*c0);
595  }
596 
597  *p0 = w0;
598  }
599  }
600  }
601 
602  } else {
603  THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
604  }
605 
606  return *this;
607  }
608 
609 
610  /**
611  * Equality.
612  *
613  * \param A matrix
614  * \param eps numerical precision
615  * \return true if matrices identical; else false
616  */
617  bool equals(const JMatrixND& A,
618  const double eps = std::numeric_limits<double>::min()) const
619  {
620  if (this->size() == A.size()) {
621 
622  for (size_t row = 0; row != this->size(); ++row) {
623 
624  const double* p = (*this)[row];
625  const double* q = A [row];
626 
627  for (size_t i = this->size(); i != 0; --i, ++p, ++q) {
628  if (fabs(*p - *q) > eps) {
629  return false;
630  }
631  }
632  }
633 
634  return true;
635 
636  } else {
637  THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
638  }
639  }
640 
641 
642  /**
643  * Test identity.
644  *
645  * \param eps numerical precision
646  * \return true if identity matrix; else false
647  */
648  bool isIdentity(const double eps = std::numeric_limits<double>::min()) const
649  {
650  for (size_t i = 0; i != this->size(); ++i) {
651 
652  if (fabs(1.0 - (*this)(i,i)) > eps) {
653  return false;
654  };
655 
656  for (size_t j = 0; j != i; ++j) {
657  if (fabs((*this)(i,j)) > eps || fabs((*this)(j,i)) > eps) {
658  return false;
659  };
660  }
661  }
662 
663  return true;
664  }
665 
666 
667  /**
668  * Test symmetry.
669  *
670  * \param eps numerical precision
671  * \return true if symmetric matrix; else false
672  */
673  bool isSymmetric(const double eps = std::numeric_limits<double>::min()) const
674  {
675  for (size_t i = 0; i != this->size(); ++i) {
676  for (size_t j = 0; j != i; ++j) {
677  if (fabs((*this)(i,j) - (*this)(j,i)) > eps) {
678  return false;
679  };
680  }
681  }
682 
683  return true;
684  }
685 
686 
687  /**
688  * Get dot product.
689  *
690  * The dot product corresponds to \f[ v^{T} \times A \times v \f].
691  *
692  * \param v vector
693  * \return dot product
694  */
695  double getDot(const JVectorND& v) const
696  {
697  double dot = 0.0;
698 
699  for (size_t i = 0; i != v.size(); ++i) {
700 
701  const double* p = (*this)[i];
702 
703  double w = 0.0;
704 
705  for (JVectorND::const_iterator y = v.begin(); y != v.end(); ++y, ++p) {
706  w += (*p) * (*y);
707  }
708 
709  dot += v[i] * w;
710  }
711 
712  return dot;
713  }
714 
715 
716  /**
717  * Read matrix from input.
718  *
719  * \param in reader
720  * \param A matrix
721  * \return reader
722  */
723  friend inline JReader& operator>>(JReader& in, JMatrixND& A)
724  {
725  size_t size;
726 
727  in >> size;
728 
729  A.resize(size);
730 
731  size_t n = A.size() * A.size();
732 
733  for (double* p = A.data(); n != 0; --n, ++p) {
734  in >> *p;
735  }
736 
737  return in;
738  }
739 
740 
741  /**
742  * Write matrix to output.
743  *
744  * \param out writer
745  * \param A matrix
746  * \return writer
747  */
748  friend inline JWriter& operator<<(JWriter& out, const JMatrixND& A)
749  {
750  out << A.size();
751 
752  size_t n = A.size() * A.size();
753 
754  for (const double* p = A.data(); n != 0; --n, ++p) {
755  out << *p;
756  }
757 
758  return out;
759  }
760 
761 
762  /**
763  * Print ASCII formatted output.
764  *
765  * \param out output stream
766  * \param A matrix
767  * \return output stream
768  */
769  friend inline std::ostream& operator<<(std::ostream& out, const JMatrixND& A)
770  {
771  using namespace std;
772 
773  const JFormat format(out, getFormat<JMatrixND>(JFormat_t(10, 3, std::ios::fixed | std::ios::showpos)));
774 
775  for (size_t row = 0; row != A.size(); ++row) {
776 
777  for (size_t col = 0; col != A.size(); ++col) {
778  out << format << A(row,col) << ' ';
779  }
780 
781  out << endl;
782  }
783 
784  return out;
785  }
786 
787 
788  protected:
789  /*
790  * Swap rows.
791  *
792  * \param r1 row
793  * \param r2 row
794  */
795  void rswap(size_t r1, size_t r2)
796  {
797  JMatrixND_t& A = getInstance();
798 
799  memcpy(A.data(), (*this)[r1], this->size() * sizeof(double));
800  memcpy((*this)[r1], (*this)[r2], this->size() * sizeof(double));
801  memcpy((*this)[r2], A.data(), this->size() * sizeof(double));
802  }
803 
804 
805  /**
806  * Swap columns.
807  *
808  * \param c1 column
809  * \param c2 column
810  */
811  void cswap(size_t c1, size_t c2)
812  {
813  using std::swap;
814 
815  double* p1 = this->data() + c1;
816  double* p2 = this->data() + c2;
817 
818  for (size_t i = this->size(); i != 0; --i, p1 += this->size(), p2 += this->size()) {
819  swap(*p1, *p2);
820  }
821  }
822  };
823 }
824 
825 #endif
JMatrixND(const JMatrixND &A)
Copy constructor.
Definition: JMatrixND.hh:336
Exception for failure of memory allocation.
Definition: JException.hh:306
General exception.
Definition: JException.hh:23
data_type w[N+1][M+1]
Definition: JPolint.hh:741
Exceptions.
Simple singleton class.
Definition: JSingleton.hh:18
Interface for binary output.
Auxiliary base class for aritmetic operations of derived class types.
Definition: JMath.hh:110
TPaveText * p1
JMatrixND & sub(const JMatrixND &A)
Matrix subtraction.
Definition: JMatrixND.hh:458
double at(size_t row, size_t col) const
Get matrix element.
Definition: JMatrixND.hh:264
void clear()
Clear memory.
Definition: JMatrixND.hh:63
JMatrixND & mul(const double factor)
Scale matrix.
Definition: JMatrixND.hh:483
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:670
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:375
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:362
const double * operator[](size_t row) const
Get row data.
Definition: JMatrixND.hh:213
JMatrixND & operator=(const JMatrixND &A)
Assignment operator.
Definition: JMatrixND.hh:347
~JMatrixND_t()
Destructor.
Definition: JMatrixND.hh:54
Auxiliary class to temporarily define format specifications.
Definition: JManip.hh:632
double operator()(const size_t row, const size_t col) const
Get matrix element.
Definition: JMatrixND.hh:238
JMatrixND(const size_t size)
Constructor.
Definition: JMatrixND.hh:324
NxN matrix.
Definition: JMatrixND.hh:303
friend JWriter & operator<<(JWriter &out, const JMatrixND &A)
Write matrix to output.
Definition: JMatrixND.hh:748
static const double C
Physics constants.
const int n
Definition: JPolint.hh:660
void set(const JMatrixND_t &A)
Set matrix.
Definition: JMatrixND.hh:110
JMatrixND_t()
Default constructor.
Definition: JMatrixND.hh:44
bool isIdentity(const double eps=std::numeric_limits< double >::min()) const
Test identity.
Definition: JMatrixND.hh:648
Template definition of auxiliary base class for comparison of data structures.
Definition: JEquals.hh:24
double * data()
Get pointer to data.
Definition: JMatrixND.hh:201
JMatrixND & add(const JMatrixND &A)
Matrix addition.
Definition: JMatrixND.hh:435
Interface for binary input.
bool equals(const JMatrixND &A, const double eps=std::numeric_limits< double >::min()) const
Equality.
Definition: JMatrixND.hh:617
void rswap(size_t r1, size_t r2)
Definition: JMatrixND.hh:795
friend JReader & operator>>(JReader &in, JMatrixND &A)
Read matrix from input.
Definition: JMatrixND.hh:723
I/O manipulators.
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:157
TCanvas * c1
Global variables to handle mouse events.
p2
Definition: module-Z:fit.sh:74
JMatrixND & negate()
Negate matrix.
Definition: JMatrixND.hh:417
double & at(size_t row, size_t col)
Get matrix element.
Definition: JMatrixND.hh:282
size_t __n
dimension of matrix
Definition: JMatrixND.hh:295
static data_type & getInstance()
Get unique instance of template class.
Definition: JSingleton.hh:27
double * __p
pointer to data
Definition: JMatrixND.hh:294
bool isSymmetric(const double eps=std::numeric_limits< double >::min()) const
Test symmetry.
Definition: JMatrixND.hh:673
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:82
JMatrixND_t & transpose()
Transpose.
Definition: JMatrixND.hh:138
Base class for data structures with artithmetic capabilities.
void cswap(size_t c1, size_t c2)
Swap columns.
Definition: JMatrixND.hh:811
size_t capacity() const
Get capacity of dimension.
Definition: JMatrixND.hh:168
JMatrixND()
Default constructor.
Definition: JMatrixND.hh:313
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:123
int j
Definition: JPolint.hh:666
Exception for accessing an index in a collection that is outside of its range.
Definition: JException.hh:90
data_type v[N+1][M+1]
Definition: JPolint.hh:740
then fatal Wrong number of arguments fi set_variable DETECTOR $argv[1] set_variable INPUT_FILE $argv[2] eval JPrintDetector a $DETECTOR O IDENTIFIER eval JPrintDetector a $DETECTOR O SUMMARY source JAcoustics sh $DETECTOR_ID CHECK_EXIT_CODE typeset A TRIPODS get_tripods $WORKDIR tripod txt TRIPODS for EMITTER in
Definition: JCanberra.sh:41
double getDot(const JVectorND &v) const
Get dot product.
Definition: JMatrixND.hh:695
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:190
bool empty() const
Check emptyness.
Definition: JMatrixND.hh:179
Data structure for format specifications.
Definition: JManip.hh:522
double & operator()(const size_t row, const size_t col)
Get matrix element.
Definition: JMatrixND.hh:251
size_t __m
capacity of matrix
Definition: JMatrixND.hh:296
p3
Definition: module-Z:fit.sh:74
Basic NxN matrix.
Definition: JMatrixND.hh:38
source $JPP_DIR setenv csh $JPP_DIR eval JShellParser o a A
JMatrixND & setIdentity()
Set to identity matrix.
Definition: JMatrixND.hh:400
double * operator[](size_t row)
Get row data.
Definition: JMatrixND.hh:225
Nx1 matrix.
Definition: JVectorND.hh:21
JMatrixND & div(const double factor)
Scale matrix.
Definition: JMatrixND.hh:501
friend std::ostream & operator<<(std::ostream &out, const JMatrixND &A)
Print ASCII formatted output.
Definition: JMatrixND.hh:769
JMatrixND & mul(const JMatrixND &A, const JMatrixND &B)
Matrix multiplication.
Definition: JMatrixND.hh:520