Jpp  18.2.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Protected Member Functions | Protected Attributes | Private Attributes | List of all members
JMATH::JMatrixNS Struct Reference

N x N symmetric matrix. More...

#include <JMatrixNS.hh>

Inheritance diagram for JMATH::JMatrixNS:
JMATH::JMatrixND JMATH::JMatrixND_t JMATH::JMath< JMatrixND > JLANG::JEquals< JFirst_t, JSecond_t > JFIT::JMatrixNZ

Public Member Functions

 JMatrixNS ()
 Default constructor. More...
 
 JMatrixNS (const size_t size)
 Constructor (identity matrix). More...
 
 JMatrixNS (const JMatrixND &A)
 Contructor. More...
 
JMatrixNSoperator= (const JMatrixNS &A)
 Assignment operator. More...
 
void invert ()
 Invert matrix according LDU decomposition. More...
 
template<class JVectorND_t >
void solve (JVectorND_t &u)
 Get solution of equation A x = b. More...
 
void update (const size_t k, const double value)
 Update inverted matrix at given diagonal element. More...
 
void resize (const size_t size)
 Resize matrix. More...
 
JMatrixNDreset ()
 Set matrix to the null matrix. More...
 
JMatrixNDsetIdentity ()
 Set to identity matrix. More...
 
JMatrixNDnegate ()
 Negate matrix. More...
 
JMatrixNDadd (const JMatrixND &A)
 Matrix addition. More...
 
JMatrixNDsub (const JMatrixND &A)
 Matrix subtraction. More...
 
JMatrixNDmul (const double factor)
 Scale matrix. More...
 
JMatrixNDmul (const JMatrixND &A, const JMatrixND &B)
 Matrix multiplication. More...
 
JMatrixNDmul (const JNullType &object)
 Multiply with object. More...
 
JMatrixNDdiv (const double factor)
 Scale matrix. More...
 
bool equals (const JMatrixND &A, const double eps=std::numeric_limits< double >::min()) const
 Equality. More...
 
bool isIdentity (const double eps=std::numeric_limits< double >::min()) const
 Test identity. More...
 
bool isSymmetric (const double eps=std::numeric_limits< double >::min()) const
 Test symmetry. More...
 
double getDot (const JVectorND &v) const
 Get dot product. More...
 
void clear ()
 Clear memory. More...
 
void set (const JMatrixND_t &A)
 Set matrix. More...
 
void swap (JMatrixND_t &A)
 Swap matrices. More...
 
JMatrixND_ttranspose ()
 Transpose. More...
 
size_t size () const
 Get dimension of matrix. More...
 
size_t capacity () const
 Get capacity of dimension. More...
 
bool empty () const
 Check emptyness. More...
 
const double * data () const
 Get pointer to data. More...
 
double * data ()
 Get pointer to data. More...
 
const double * operator[] (size_t row) const
 Get row data. More...
 
double * operator[] (size_t row)
 Get row data. More...
 
double operator() (const size_t row, const size_t col) const
 Get matrix element. More...
 
double & operator() (const size_t row, const size_t col)
 Get matrix element. More...
 
double at (size_t row, size_t col) const
 Get matrix element. More...
 
double & at (size_t row, size_t col)
 Get matrix element. More...
 

Protected Member Functions

JMatrixND_tgetInstance ()
 Get work space. More...
 
void rswap (size_t r1, size_t r2)
 
void cswap (size_t c1, size_t c2)
 Swap columns. More...
 

Protected Attributes

JMatrixND_t ws
 
double * __p
 pointer to data More...
 
size_t __n
 dimension of matrix More...
 
size_t __m
 capacity of matrix More...
 

Private Attributes

std::vector< int > P
 permutation matrix More...
 

Detailed Description

N x N symmetric matrix.

Definition at line 28 of file JMatrixNS.hh.

Constructor & Destructor Documentation

JMATH::JMatrixNS::JMatrixNS ( )
inline

Default constructor.

Definition at line 34 of file JMatrixNS.hh.

34  :
35  JMatrixND()
36  {}
JMatrixND()
Default constructor.
Definition: JMatrixND.hh:370
JMATH::JMatrixNS::JMatrixNS ( const size_t  size)
inline

Constructor (identity matrix).

Parameters
sizedimension

Definition at line 44 of file JMatrixNS.hh.

44  :
46  {}
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
JMatrixND()
Default constructor.
Definition: JMatrixND.hh:370
JMATH::JMatrixNS::JMatrixNS ( const JMatrixND A)
inline

Contructor.

Parameters
Amatrix

Definition at line 54 of file JMatrixNS.hh.

54  :
55  JMatrixND(A)
56  {}
JMatrixND()
Default constructor.
Definition: JMatrixND.hh:370

Member Function Documentation

JMatrixNS& JMATH::JMatrixNS::operator= ( const JMatrixNS A)
inline

Assignment operator.

Parameters
Amatrix

Definition at line 64 of file JMatrixNS.hh.

65  {
66  this->set(A);
67 
68  return *this;
69  }
void set(const JMatrixND_t &A)
Set matrix.
Definition: JMatrixND.hh:152
void JMATH::JMatrixNS::invert ( )
inline

Invert matrix according LDU decomposition.

Definition at line 75 of file JMatrixNS.hh.

76  {
77  using namespace std;
78 
79  size_t i, row, col, root;
80 
81  double* p0;
82  double* p1;
83  double* p2;
84  double* p3;
85 
86  const double* c0;
87 
88  double v0;
89  double v1;
90  double v2;
91  double v3;
92 
93  double w0;
94 
95  P.resize(this->size());
96 
97  // LDU decomposition
98 
99  for (root = 0; root != this->size(); ++root) {
100 
101  double value = (*this)(root, root);
102  size_t pivot = root;
103 
104  for (row = root; ++row != this->size(); ) {
105  if (fabs((*this)(row, root)) > fabs(value)) {
106  value = (*this)(row, root);
107  pivot = row;
108  }
109  }
110 
111  if (value == 0.0) {
112  THROW(JDivisionByZero, "LDU decomposition zero pivot at " << root);
113  }
114 
115  if (pivot != root) {
116  this->rswap(root, pivot);
117  }
118 
119  P[root] = pivot;
120 
121  for (row = root + 1; row + 4 <= this->size(); row += 4) { // process rows by four
122 
123  p0 = (*this)[row + 0] + root;
124  p1 = (*this)[row + 1] + root;
125  p2 = (*this)[row + 2] + root;
126  p3 = (*this)[row + 3] + root;
127 
128  v0 = (*p0++ /= value);
129  v1 = (*p1++ /= value);
130  v2 = (*p2++ /= value);
131  v3 = (*p3++ /= value);
132 
133  c0 = (*this)[root] + root + 1;
134 
135  for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, p1 += 4, p2 += 4, p3 += 4, c0 += 4) {
136  p0[0] -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
137  p1[0] -= v1 * c0[0]; p1[1] -= v1 * c0[1]; p1[2] -= v1 * c0[2]; p1[3] -= v1 * c0[3];
138  p2[0] -= v2 * c0[0]; p2[1] -= v2 * c0[1]; p2[2] -= v2 * c0[2]; p2[3] -= v2 * c0[3];
139  p3[0] -= v3 * c0[0]; p3[1] -= v3 * c0[1]; p3[2] -= v3 * c0[2]; p3[3] -= v3 * c0[3];
140  }
141 
142  for ( ; i != this->size(); ++i, ++p0, ++p1, ++p2, ++p3, ++c0) {
143  *p0 -= v0 * (*c0);
144  *p1 -= v1 * (*c0);
145  *p2 -= v2 * (*c0);
146  *p3 -= v3 * (*c0);
147  }
148  }
149 
150  for ( ; row != this->size(); ++row) { // remaining rows
151 
152  p0 = (*this)[row + 0] + root;
153 
154  v0 = (*p0++ /= value);
155 
156  c0 = (*this)[root] + root + 1;
157 
158  for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, c0 += 4) {
159  *p0 -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
160  }
161 
162  for ( ; i != this->size(); ++i, ++p0, ++c0) {
163  *p0 -= v0 * (*c0);
164  }
165  }
166  }
167 
168 
170 
171  A.reset();
172 
173  for (i = 0; i != this->size(); ++i) {
174  A(i,i) = 1.0;
175  }
176 
177 
178  // forward substitution
179 
180  col = 0;
181 
182  for ( ; col + 4 <= this->size(); col += 4) {
183 
184  for (row = 0; row != this->size(); ++row) {
185 
186  p0 = A[col + 0];
187  p1 = A[col + 1];
188  p2 = A[col + 2];
189  p3 = A[col + 3];
190 
191  i = P[row];
192 
193  v0 = p0[i];
194  v1 = p1[i];
195  v2 = p2[i];
196  v3 = p3[i];
197 
198  p0[i] = p0[row];
199  p1[i] = p1[row];
200  p2[i] = p2[row];
201  p3[i] = p3[row];
202 
203  for (i = 0, c0 = (*this)[row]; i != row; ++c0, ++p0, ++p1, ++p2, ++p3, ++i) {
204  v0 -= (*c0) * (*p0);
205  v1 -= (*c0) * (*p1);
206  v2 -= (*c0) * (*p2);
207  v3 -= (*c0) * (*p3);
208  }
209 
210  A[col + 0][row] = v0;
211  A[col + 1][row] = v1;
212  A[col + 2][row] = v2;
213  A[col + 3][row] = v3;
214  }
215  }
216 
217  for ( ; col != this->size(); ++col) {
218 
219  for (row = 0; row != this->size(); ++row) {
220 
221  p0 = A[col];
222 
223  i = P[row];
224 
225  v0 = p0[i];
226 
227  p0[i] = p0[row];
228 
229  for (i = 0, c0 = (*this)[row] + i; i != row; ++c0, ++p0, ++i) {
230  v0 -= (*c0) * (*p0);
231  }
232 
233  A[col][row] = v0;
234  }
235  }
236 
237  // backward substitution
238 
239  col = 0;
240 
241  for ( ; col + 4 <= this->size(); col += 4) {
242 
243  for (int row = this->size() - 1; row >= 0; --row) {
244 
245  p0 = A[col + 0] + row;
246  p1 = A[col + 1] + row;
247  p2 = A[col + 2] + row;
248  p3 = A[col + 3] + row;
249 
250  v0 = *p0;
251  v1 = *p1;
252  v2 = *p2;
253  v3 = *p3;
254 
255  ++p0;
256  ++p1;
257  ++p2;
258  ++p3;
259 
260  for (i = row + 1, c0 = (*this)[row] + i; i != this->size(); ++c0, ++p0, ++p1, ++p2, ++p3, ++i) {
261  v0 -= (*c0) * (*p0);
262  v1 -= (*c0) * (*p1);
263  v2 -= (*c0) * (*p2);
264  v3 -= (*c0) * (*p3);
265  }
266 
267  w0 = 1.0 / (*this)[row][row];
268 
269  A[col + 0][row] = v0 * w0;
270  A[col + 1][row] = v1 * w0;
271  A[col + 2][row] = v2 * w0;
272  A[col + 3][row] = v3 * w0;
273  }
274  }
275 
276  for ( ; col != this->size(); ++col) {
277 
278  for (int row = this->size() - 1; row >= 0; --row) {
279 
280  p0 = A[col] + row;
281 
282  v0 = *p0;
283 
284  ++p0;
285 
286  for (i = row + 1, c0 = (*this)[row] + i; i != this->size(); ++c0, ++p0, ++i) {
287  v0 -= (*c0) * (*p0);
288  }
289 
290  w0 = 1.0 / (*this)[row][row];
291 
292  A[col][row] = v0 * w0;
293  }
294  }
295 
296  A.transpose();
297 
298  this->swap(A);
299  }
TPaveText * p1
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:358
std::vector< int > P
permutation matrix
Definition: JMatrixNS.hh:486
void rswap(size_t r1, size_t r2)
Definition: JMatrixND.hh:870
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
p2
Definition: module-Z:fit.sh:74
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:165
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
p3
Definition: module-Z:fit.sh:74
Basic NxN matrix.
Definition: JMatrixND.hh:36
do JPlot2D f $WORKDIR canberra[${EMITTER}] root
Definition: JCanberra.sh:132
template<class JVectorND_t >
void JMATH::JMatrixNS::solve ( JVectorND_t &  u)
inline

Get solution of equation A x = b.

Parameters
ucolumn vector; b on input, x on output(I/O)

Definition at line 308 of file JMatrixNS.hh.

309  {
310  using namespace std;
311 
312  size_t i, row, root;
313 
314  double* p0;
315  double* p1;
316  double* p2;
317  double* p3;
318 
319  const double* c0;
320 
321  double v0;
322  double v1;
323  double v2;
324  double v3;
325 
326  P.resize(this->size());
327 
328  // LDU decomposition
329 
330  for (root = 0; root != this->size(); ++root) {
331 
332  double value = (*this)(root, root);
333  size_t pivot = root;
334 
335  for (row = root; ++row != this->size(); ) {
336  if (fabs((*this)(row, root)) > fabs(value)) {
337  value = (*this)(row, root);
338  pivot = row;
339  }
340  }
341 
342  if (value == 0.0) {
343  THROW(JDivisionByZero, "LDU decomposition zero pivot at " << root);
344  }
345 
346  if (pivot != root) {
347  this->rswap(root, pivot);
348  }
349 
350  P[root] = pivot;
351 
352  for (row = root + 1; row + 4 <= this->size(); row += 4) { // process rows by four
353 
354  p0 = (*this)[row + 0] + root;
355  p1 = (*this)[row + 1] + root;
356  p2 = (*this)[row + 2] + root;
357  p3 = (*this)[row + 3] + root;
358 
359  v0 = (*p0++ /= value);
360  v1 = (*p1++ /= value);
361  v2 = (*p2++ /= value);
362  v3 = (*p3++ /= value);
363 
364  c0 = (*this)[root] + root + 1;
365 
366  for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, p1 += 4, p2 += 4, p3 += 4, c0 += 4) {
367  p0[0] -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
368  p1[0] -= v1 * c0[0]; p1[1] -= v1 * c0[1]; p1[2] -= v1 * c0[2]; p1[3] -= v1 * c0[3];
369  p2[0] -= v2 * c0[0]; p2[1] -= v2 * c0[1]; p2[2] -= v2 * c0[2]; p2[3] -= v2 * c0[3];
370  p3[0] -= v3 * c0[0]; p3[1] -= v3 * c0[1]; p3[2] -= v3 * c0[2]; p3[3] -= v3 * c0[3];
371  }
372 
373  for ( ; i != this->size(); ++i, ++p0, ++p1, ++p2, ++p3, ++c0) {
374  *p0 -= v0 * (*c0);
375  *p1 -= v1 * (*c0);
376  *p2 -= v2 * (*c0);
377  *p3 -= v3 * (*c0);
378  }
379  }
380 
381  for ( ; row != this->size(); ++row) { // remaining rows
382 
383  p0 = (*this)[row + 0] + root;
384 
385  v0 = (*p0++ /= value);
386 
387  c0 = (*this)[root] + root + 1;
388 
389  for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, c0 += 4) {
390  *p0 -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
391  }
392 
393  for ( ; i != this->size(); ++i, ++p0, ++c0) {
394  *p0 -= v0 * (*c0);
395  }
396  }
397  }
398 
399  // forward substitution
400 
401  for (row = 0; row != this->size(); ++row) {
402 
403  i = P[row];
404  v0 = u[i];
405  u[i] = u[row];
406 
407  for (i = 0, c0 = (*this)[row] + i; i != row; ++c0, ++i) {
408  v0 -= *c0 * u[i];
409  }
410 
411  u[row] = v0;
412  }
413 
414  // backward substitution
415 
416  for (int row = this->size() - 1; row >= 0; --row) {
417 
418  v0 = u[row];
419 
420  for (i = row + 1, c0 = (*this)[row] + i; i < this->size(); ++c0, ++i) {
421  v0 -= *c0 * u[i];
422  }
423 
424  u[row] = v0 / (*this)[row][row];
425  }
426  }
TPaveText * p1
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
std::vector< int > P
permutation matrix
Definition: JMatrixNS.hh:486
void rswap(size_t r1, size_t r2)
Definition: JMatrixND.hh:870
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
p2
Definition: module-Z:fit.sh:74
double u[N+1]
Definition: JPolint.hh:821
p3
Definition: module-Z:fit.sh:74
do JPlot2D f $WORKDIR canberra[${EMITTER}] root
Definition: JCanberra.sh:132
void JMATH::JMatrixNS::update ( const size_t  k,
const double  value 
)
inline

Update inverted matrix at given diagonal element.

If A is the original matrix and A' is such that for each (i,j) != (k,k):

      A'[i][j] = A[i][j]
      A'[k][k] = A[k][k] + value

then JMatrixNS::update(k, value) will modify the already inverted matrix of A such that it will be the inverse of A'.

See also arXiv.
This algorithm can be considered a special case of the Sherman–Morrison formula.

Parameters
kindex of diagonal element
valuevalue to add

Definition at line 446 of file JMatrixNS.hh.

447  {
448  if (k >= this->size()) {
449  THROW(JIndexOutOfRange, "JMatrixNS::update(): index out of range " << k);
450  }
451 
452  size_t i, row;
453  double* col;
454 
455  const double del = (*this)(k,k);
456  const double din = 1.0 / (1.0 + value * del);
457 
458  double* root = (*this)[k];
459 
460  for (row = 0; row != this->size(); ++row) {
461 
462  if (row != k) {
463 
464  const double val = (*this)(row, k);
465 
466  for (i = 0, col = (*this)[row]; i != this->size(); ++i, ++col) {
467  if (i != k) {
468  (*col) -= val * root[i] * value * din;
469  }
470  }
471 
472  (*this)(row, k) *= din;
473  }
474  }
475 
476  for (i = 0, col = root; i != this->size(); ++i, ++col) {
477  if (i != k) {
478  (*col) *= din;
479  }
480  }
481 
482  root[k] = del * din;
483  }
then fatal No hydrophone data file $HYDROPHONE_TXT fi sort gr k
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
do JPlot2D f $WORKDIR canberra[${EMITTER}] root
Definition: JCanberra.sh:132
JMatrixND_t& JMATH::JMatrixND::getInstance ( )
inlineprotectedinherited

Get work space.

Returns
work space

Definition at line 358 of file JMatrixND.hh.

359  {
360  return ws;
361  }
JMatrixND_t ws
Definition: JMatrixND.hh:351
void JMATH::JMatrixND::resize ( const size_t  size)
inlineinherited

Resize matrix.

Note that this method does not maintain data in the matrix.

Parameters
sizedimension

Definition at line 443 of file JMatrixND.hh.

444  {
445  static_cast<JMatrixND_t&>(*this).resize(size);
446 
448  }
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:358
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:102
Basic NxN matrix.
Definition: JMatrixND.hh:36
JMatrixND& JMATH::JMatrixND::reset ( )
inlineinherited

Set matrix to the null matrix.

Returns
this matrix

Definition at line 456 of file JMatrixND.hh.

457  {
459 
461 
462  A.resize(this->size());
463 
464  return *this;
465  }
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:358
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:102
JMatrixND_t & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:130
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
Basic NxN matrix.
Definition: JMatrixND.hh:36
JMatrixND& JMATH::JMatrixND::setIdentity ( )
inlineinherited

Set to identity matrix.

Returns
this matrix

Definition at line 473 of file JMatrixND.hh.

474  {
475  reset();
476 
477  for (size_t i = 0; i != this->size(); ++i) {
478  (*this)(i,i) = 1.0;
479  }
480 
481  return *this;
482  }
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:456
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
JMatrixND& JMATH::JMatrixND::negate ( )
inlineinherited

Negate matrix.

Returns
this matrix

Definition at line 490 of file JMatrixND.hh.

491  {
492  double* p = this->data();
493 
494  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
495  *p = -*p;
496  }
497 
498  return *this;
499  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:232
JMatrixND& JMATH::JMatrixND::add ( const JMatrixND A)
inlineinherited

Matrix addition.

Parameters
Amatrix
Returns
this matrix

Definition at line 508 of file JMatrixND.hh.

509  {
510  if (this->size() == A.size()) {
511 
512  double* p = this->data();
513  const double* q = A.data();
514 
515  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
516  *p += *q;
517  }
518 
519  } else {
520  THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
521  }
522 
523  return *this;
524  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:232
JMatrixND& JMATH::JMatrixND::sub ( const JMatrixND A)
inlineinherited

Matrix subtraction.

Parameters
Amatrix
Returns
this matrix

Definition at line 533 of file JMatrixND.hh.

534  {
535  if (this->size() == A.size()) {
536 
537  double* p = this->data();
538  const double* q = A.data();
539 
540  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
541  *p -= *q;
542  }
543 
544  } else {
545  THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
546  }
547 
548  return *this;
549  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:232
JMatrixND& JMATH::JMatrixND::mul ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 558 of file JMatrixND.hh.

559  {
560  double* p = this->data();
561 
562  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
563  *p *= factor;
564  }
565 
566  return *this;
567  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:232
JMatrixND& JMATH::JMatrixND::mul ( const JMatrixND A,
const JMatrixND B 
)
inlineinherited

Matrix multiplication.

Parameters
Amatrix
Bmatrix
Returns
this matrix

Definition at line 595 of file JMatrixND.hh.

597  {
598  if (A.size() == B.size()) {
599 
600  this->resize(A.size());
601 
602  if (!this->empty()) {
603 
605 
606  C.set(B);
607  C.transpose();
608 
609  size_t i, row;
610 
611  for (row = 0; row + 4 <= this->size(); row += 4) { // process rows by four
612 
613  double* p0 = (*this)[row + 0];
614  double* p1 = (*this)[row + 1];
615  double* p2 = (*this)[row + 2];
616  double* p3 = (*this)[row + 3];
617 
618  for (size_t col = 0; col != this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
619 
620  double w0 = 0.0;
621  double w1 = 0.0;
622  double w2 = 0.0;
623  double w3 = 0.0;
624 
625  const double* a0 = A[row + 0];
626  const double* a1 = A[row + 1];
627  const double* a2 = A[row + 2];
628  const double* a3 = A[row + 3];
629 
630  const double* c0 = C[col];
631 
632  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4, c0 += 4) {
633  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
634  w1 += a1[0] * c0[0] + a1[1] * c0[1] + a1[2] * c0[2] + a1[3] * c0[3];
635  w2 += a2[0] * c0[0] + a2[1] * c0[1] + a2[2] * c0[2] + a2[3] * c0[3];
636  w3 += a3[0] * c0[0] + a3[1] * c0[1] + a3[2] * c0[2] + a3[3] * c0[3];
637  }
638 
639  for ( ; i != this->size(); ++i, ++a0, ++a1, ++a2, ++a3, ++c0) {
640  w0 += (*a0) * (*c0);
641  w1 += (*a1) * (*c0);
642  w2 += (*a2) * (*c0);
643  w3 += (*a3) * (*c0);
644  }
645 
646  *p0 = w0;
647  *p1 = w1;
648  *p2 = w2;
649  *p3 = w3;
650  }
651  }
652 
653  for ( ; row != this->size(); ++row) { // remaining rows
654 
655  double* p0 = (*this)[row + 0];
656 
657  for (size_t col = 0; col != this->size(); ++col, ++p0) {
658 
659  double w0 = 0.0;
660 
661  const double* a0 = A[row + 0];
662  const double* c0 = C[col];
663 
664  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, c0 += 4) {
665  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
666  }
667 
668  for ( ; i != this->size(); ++i, ++a0, ++c0) {
669  w0 += (*a0) * (*c0);
670  }
671 
672  *p0 = w0;
673  }
674  }
675  }
676 
677  } else {
678  THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
679  }
680 
681  return *this;
682  }
TPaveText * p1
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:443
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:358
static const double C
Physics constants.
void set(const JMatrixND_t &A)
Set matrix.
Definition: JMatrixND.hh:152
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
p2
Definition: module-Z:fit.sh:74
JMatrixND_t & transpose()
Transpose.
Definition: JMatrixND.hh:180
bool empty() const
Check emptyness.
Definition: JMatrixND.hh:221
p3
Definition: module-Z:fit.sh:74
Basic NxN matrix.
Definition: JMatrixND.hh:36
JMatrixND & JMATH::JMath< JMatrixND , JNullType >::mul ( const JNullType object)
inlineinherited

Multiply with object.

Parameters
objectobject
Returns
result object

Definition at line 354 of file JMath.hh.

355  {
356  return static_cast<JFirst_t&>(*this) = JFirst_t().mul(static_cast<const JFirst_t&>(*this), object);
357  }
JMatrixND& JMATH::JMatrixND::div ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 576 of file JMatrixND.hh.

577  {
578  double* p = this->data();
579 
580  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
581  *p /= factor;
582  }
583 
584  return *this;
585  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:232
bool JMATH::JMatrixND::equals ( const JMatrixND A,
const double  eps = std::numeric_limits<double>::min() 
) const
inlineinherited

Equality.

Parameters
Amatrix
epsnumerical precision
Returns
true if matrices identical; else false

Definition at line 692 of file JMatrixND.hh.

694  {
695  if (this->size() == A.size()) {
696 
697  for (size_t row = 0; row != this->size(); ++row) {
698 
699  const double* p = (*this)[row];
700  const double* q = A [row];
701 
702  for (size_t i = this->size(); i != 0; --i, ++p, ++q) {
703  if (fabs(*p - *q) > eps) {
704  return false;
705  }
706  }
707  }
708 
709  return true;
710 
711  } else {
712  THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
713  }
714  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
bool JMATH::JMatrixND::isIdentity ( const double  eps = std::numeric_limits<double>::min()) const
inlineinherited

Test identity.

Parameters
epsnumerical precision
Returns
true if identity matrix; else false

Definition at line 723 of file JMatrixND.hh.

724  {
725  for (size_t i = 0; i != this->size(); ++i) {
726 
727  if (fabs(1.0 - (*this)(i,i)) > eps) {
728  return false;
729  };
730 
731  for (size_t j = 0; j != i; ++j) {
732  if (fabs((*this)(i,j)) > eps || fabs((*this)(j,i)) > eps) {
733  return false;
734  };
735  }
736  }
737 
738  return true;
739  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
int j
Definition: JPolint.hh:748
bool JMATH::JMatrixND::isSymmetric ( const double  eps = std::numeric_limits<double>::min()) const
inlineinherited

Test symmetry.

Parameters
epsnumerical precision
Returns
true if symmetric matrix; else false

Definition at line 748 of file JMatrixND.hh.

749  {
750  for (size_t i = 0; i != this->size(); ++i) {
751  for (size_t j = 0; j != i; ++j) {
752  if (fabs((*this)(i,j) - (*this)(j,i)) > eps) {
753  return false;
754  };
755  }
756  }
757 
758  return true;
759  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
int j
Definition: JPolint.hh:748
double JMATH::JMatrixND::getDot ( const JVectorND v) const
inlineinherited

Get dot product.

The dot product corresponds to

\[ v^{T} \times A \times v \]

.

Parameters
vvector
Returns
dot product

Definition at line 770 of file JMatrixND.hh.

771  {
772  double dot = 0.0;
773 
774  for (size_t i = 0; i != v.size(); ++i) {
775 
776  const double* p = (*this)[i];
777 
778  double w = 0.0;
779 
780  for (JVectorND::const_iterator y = v.begin(); y != v.end(); ++y, ++p) {
781  w += (*p) * (*y);
782  }
783 
784  dot += v[i] * w;
785  }
786 
787  return dot;
788  }
data_type w[N+1][M+1]
Definition: JPolint.hh:823
void JMATH::JMatrixND::rswap ( size_t  r1,
size_t  r2 
)
inlineprotectedinherited

Definition at line 870 of file JMatrixND.hh.

871  {
873 
874  A.resize(this->size());
875 
876  memcpy(A.data(), (*this)[r1], this->size() * sizeof(double));
877  memcpy((*this)[r1], (*this)[r2], this->size() * sizeof(double));
878  memcpy((*this)[r2], A.data(), this->size() * sizeof(double));
879  }
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:358
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:102
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:232
Basic NxN matrix.
Definition: JMatrixND.hh:36
void JMATH::JMatrixND::cswap ( size_t  c1,
size_t  c2 
)
inlineprotectedinherited

Swap columns.

Parameters
c1column
c2column

Definition at line 888 of file JMatrixND.hh.

889  {
890  using std::swap;
891 
892  double* p1 = this->data() + c1;
893  double* p2 = this->data() + c2;
894 
895  for (size_t i = this->size(); i != 0; --i, p1 += this->size(), p2 += this->size()) {
896  swap(*p1, *p2);
897  }
898  }
TPaveText * p1
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
TCanvas * c1
Global variables to handle mouse events.
p2
Definition: module-Z:fit.sh:74
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:165
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:232
void JMATH::JMatrixND_t::clear ( )
inlineinherited

Clear memory.

Definition at line 83 of file JMatrixND.hh.

84  {
85  if (__p != NULL) {
86  delete [] __p;
87  }
88 
89  __p = NULL;
90  __n = 0;
91  __m = 0;
92  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:337
double * __p
pointer to data
Definition: JMatrixND.hh:336
size_t __m
capacity of matrix
Definition: JMatrixND.hh:338
void JMATH::JMatrixND_t::set ( const JMatrixND_t A)
inlineinherited

Set matrix.

Parameters
Amatrix

Definition at line 152 of file JMatrixND.hh.

153  {
154  this->resize(A.size());
155 
156  memcpy(this->data(), A.data(), A.size() * A.size() * sizeof(double));
157  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:102
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:232
void JMATH::JMatrixND_t::swap ( JMatrixND_t A)
inlineinherited

Swap matrices.

Parameters
Amatrix

Definition at line 165 of file JMatrixND.hh.

166  {
167  using std::swap;
168 
169  swap(this->__p, A.__p);
170  swap(this->__n, A.__n);
171  swap(this->__m, A.__m);
172  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:337
double * __p
pointer to data
Definition: JMatrixND.hh:336
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:165
size_t __m
capacity of matrix
Definition: JMatrixND.hh:338
JMatrixND_t& JMATH::JMatrixND_t::transpose ( )
inlineinherited

Transpose.

Returns
this matrix

Definition at line 180 of file JMatrixND.hh.

181  {
182  using std::swap;
183 
184  for (size_t row = 0; row != this->size(); ++row) {
185  for (size_t col = 0; col != row; ++col) {
186  swap((*this)(row,col), (*this)(col,row));
187  }
188  }
189 
190  return *this;
191  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:165
size_t JMATH::JMatrixND_t::size ( ) const
inlineinherited

Get dimension of matrix.

Returns
dimension

Definition at line 199 of file JMatrixND.hh.

200  {
201  return __n;
202  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:337
size_t JMATH::JMatrixND_t::capacity ( ) const
inlineinherited

Get capacity of dimension.

Returns
capacity

Definition at line 210 of file JMatrixND.hh.

211  {
212  return __m;
213  }
size_t __m
capacity of matrix
Definition: JMatrixND.hh:338
bool JMATH::JMatrixND_t::empty ( ) const
inlineinherited

Check emptyness.

Returns
true if empty; else false

Definition at line 221 of file JMatrixND.hh.

222  {
223  return __n == 0;
224  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:337
const double* JMATH::JMatrixND_t::data ( ) const
inlineinherited

Get pointer to data.

Returns
pointer to data.

Definition at line 232 of file JMatrixND.hh.

233  {
234  return __p;
235  }
double * __p
pointer to data
Definition: JMatrixND.hh:336
double* JMATH::JMatrixND_t::data ( )
inlineinherited

Get pointer to data.

Returns
pointer to data.

Definition at line 243 of file JMatrixND.hh.

244  {
245  return __p;
246  }
double * __p
pointer to data
Definition: JMatrixND.hh:336
const double* JMATH::JMatrixND_t::operator[] ( size_t  row) const
inlineinherited

Get row data.

Parameters
rowrow number
Returns
pointer to row data

Definition at line 255 of file JMatrixND.hh.

256  {
257  return __p + row*__n;
258  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:337
double * __p
pointer to data
Definition: JMatrixND.hh:336
double* JMATH::JMatrixND_t::operator[] ( size_t  row)
inlineinherited

Get row data.

Parameters
rowrow number
Returns
pointer to row data

Definition at line 267 of file JMatrixND.hh.

268  {
269  return __p + row*__n;
270  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:337
double * __p
pointer to data
Definition: JMatrixND.hh:336
double JMATH::JMatrixND_t::operator() ( const size_t  row,
const size_t  col 
) const
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 280 of file JMatrixND.hh.

281  {
282  return *(__p + row*__n + col);
283  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:337
double * __p
pointer to data
Definition: JMatrixND.hh:336
double& JMATH::JMatrixND_t::operator() ( const size_t  row,
const size_t  col 
)
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 293 of file JMatrixND.hh.

294  {
295  return *(__p + row*__n + col);
296  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:337
double * __p
pointer to data
Definition: JMatrixND.hh:336
double JMATH::JMatrixND_t::at ( size_t  row,
size_t  col 
) const
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 306 of file JMatrixND.hh.

307  {
308  if (row >= this->size() ||
309  col >= this->size()) {
310  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
311  }
312 
313  return (*this)(row, col);
314  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199
double& JMATH::JMatrixND_t::at ( size_t  row,
size_t  col 
)
inlineinherited

Get matrix element.

Parameters
rowrow number
colcolumn number
Returns
matrix element at (row,col)

Definition at line 324 of file JMatrixND.hh.

325  {
326  if (row >= this->size() ||
327  col >= this->size()) {
328  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
329  }
330 
331  return (*this)(row, col);
332  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:199

Member Data Documentation

std::vector<int> JMATH::JMatrixNS::P
private

permutation matrix

Definition at line 486 of file JMatrixNS.hh.

JMatrixND_t JMATH::JMatrixND::ws
protectedinherited

Definition at line 351 of file JMatrixND.hh.

double* JMATH::JMatrixND_t::__p
protectedinherited

pointer to data

Definition at line 336 of file JMatrixND.hh.

size_t JMATH::JMatrixND_t::__n
protectedinherited

dimension of matrix

Definition at line 337 of file JMatrixND.hh.

size_t JMATH::JMatrixND_t::__m
protectedinherited

capacity of matrix

Definition at line 338 of file JMatrixND.hh.


The documentation for this struct was generated from the following file: