Jpp  debug
the software that should make you happy
Classes | Public Member Functions | Protected Member Functions | Protected Attributes | Static Private Member Functions | Private Attributes | List of all members
JFIT::JMatrixNZ Class Reference

Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z). More...

#include <JMatrixNZ.hh>

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

Classes

struct  variance
 Auxiliary data structure for co-variance calculation. More...
 

Public Member Functions

 JMatrixNZ ()
 Default contructor. More...
 
template<class T >
 JMatrixNZ (const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
 Constructor. More...
 
template<class T >
void set (const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
 Set co-variance matrix. 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...
 
JMatrixND & reset ()
 Set matrix to the null matrix. More...
 
JMatrixND & setIdentity ()
 Set to identity matrix. More...
 
JMatrixND & negate ()
 Negate matrix. More...
 
JMatrixND & add (const JMatrixND &A)
 Matrix addition. More...
 
JMatrixND & sub (const JMatrixND &A)
 Matrix subtraction. More...
 
JMatrixND & mul (const double factor)
 Scale matrix. More...
 
JMatrixND & mul (const JMatrixND &A, const JMatrixND &B)
 Matrix multiplication. More...
 
JMatrixND & mul (const JSecond_t &object)
 Multiply with object. More...
 
JMatrixND & div (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_t & transpose ()
 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_t & getInstance ()
 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...
 

Static Private Member Functions

static double getDot (const variance &first, const variance &second)
 Get dot product. More...
 

Private Attributes

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

Detailed Description

Determination of the co-variance matrix of hits for a track along z-axis (JFIT::JLine1Z).

In this, the given angular and time resolution are taken into account.

Definition at line 28 of file JMatrixNZ.hh.

Constructor & Destructor Documentation

◆ JMatrixNZ() [1/2]

JFIT::JMatrixNZ::JMatrixNZ ( )
inline

Default contructor.

Definition at line 39 of file JMatrixNZ.hh.

39  :
40  JMatrixNS()
41  {}
JMATH::JMatrixNS::JMatrixNS
JMatrixNS()
Default constructor.
Definition: JMatrixNS.hh:34

◆ JMatrixNZ() [2/2]

template<class T >
JFIT::JMatrixNZ::JMatrixNZ ( const JVector3D &  pos,
T  __begin,
T  __end,
const double  alpha,
const double  sigma 
)
inline

Constructor.

The template argument T refers to an iterator of a data structure which should have the following member methods:

  • double getX(); // [m]
  • double getY(); // [m]
  • double getZ(); // [m]
Parameters
posreference position [m]
__beginbegin of data
__endend of data
alphaangular resolution [deg]
sigmatime resolution [ns]

Definition at line 59 of file JMatrixNZ.hh.

63  :
64  JMatrixNS()
65  {
66  set(pos, __begin, __end, alpha, sigma);
67  }
JFIT::JMatrixNZ::set
void set(const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
Set co-variance matrix.
Definition: JMatrixNZ.hh:85
GAUSS_HERMITE::sigma
const double sigma[]
Definition: JQuadrature.cc:74

Member Function Documentation

◆ set() [1/2]

template<class T >
void JFIT::JMatrixNZ::set ( const JVector3D &  pos,
T  __begin,
T  __end,
const double  alpha,
const double  sigma 
)
inline

Set co-variance matrix.

The template argument T refers to an iterator of a data structure which should have the following member methods:

  • double getX(); // [m]
  • double getY(); // [m]
  • double getZ(); // [m]
Parameters
posreference position [m]
__beginbegin of data
__endend of data
alphaangular resolution [deg]
sigmatime resolution [ns]

Definition at line 85 of file JMatrixNZ.hh.

90  {
91  using namespace std;
92  using namespace JPP;
93 
94 
95  const int N = distance(__begin, __end);
96 
97  this ->resize(N);
98  buffer.resize(N);
99 
100  const double ta = alpha * PI / 180.0;
101  const double ct = cos(ta);
102  const double st = sin(ta);
103 
104 
105  // angular resolution
106 
107  for (T i = __begin; i != __end; ++i) {
108 
109  const double dx = i->getX() - pos.getX();
110  const double dy = i->getY() - pos.getY();
111  const double dz = i->getZ() - pos.getZ();
112 
113  const double R = sqrt(dx*dx + dy*dy);
114 
115  double x = ta * getKappaC() * getInverseSpeedOfLight();
116  double y = ta * getKappaC() * getInverseSpeedOfLight();
117  double v = ta * getInverseSpeedOfLight();
118  double w = ta * getInverseSpeedOfLight();
119 
120  if (R != 0.0) {
121  x *= dx / R;
122  y *= dy / R;
123  }
124 
125  x *= (dz*ct - dx*st);
126  y *= (dz*ct - dy*st);
127  v *= -(dx*ct + dz*st);
128  w *= -(dy*ct + dz*st);
129 
130  buffer[distance(__begin,i)] = variance(x, y, v, w);
131  }
132 
133  for (int i = 0; i != N; ++i) {
134 
135  for (int j = 0; j != i; ++j) {
136  (*this)(i, j) = getDot(buffer[i], buffer[j]);
137  (*this)(j, i) = (*this)(i, j);
138  }
139 
140  (*this)(i, i) = getDot(buffer[i], buffer[i]) + sigma * sigma;
141  }
142  }
distance
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
Definition: PhysicsEvent.hh:434
JFIT::JMatrixNZ::buffer
std::vector< variance > buffer
Definition: JMatrixNZ.hh:186
JFIT::JMatrixNZ::getDot
static double getDot(const variance &first, const variance &second)
Get dot product.
Definition: JMatrixNZ.hh:196
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition: JVector3D.hh:104
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition: JVector3D.hh:115
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition: JVector3D.hh:94
JMATH::PI
static const double PI
Mathematical constants.
Definition: JMath/JConstants.hh:20
JPHYSICS::getKappaC
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
Definition: JPhysics/JConstants.hh:189
JPHYSICS::getInverseSpeedOfLight
const double getInverseSpeedOfLight()
Get inverse speed of light.
Definition: JPhysics/JConstants.hh:123
JPP
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JAAnetToolkit.hh:43
JTOOLS::w
data_type w[N+1][M+1]
Definition: JPolint.hh:867
JTOOLS::v
data_type v[N+1][M+1]
Definition: JPolint.hh:866
JTOOLS::j
int j
Definition: JPolint.hh:792
std
Definition: JSTDTypes.hh:14

◆ getDot() [1/2]

static double JFIT::JMatrixNZ::getDot ( const variance &  first,
const variance &  second 
)
inlinestaticprivate

Get dot product.

Parameters
firstfirst variance
secondsecond variance
Returns
dot product

Definition at line 196 of file JMatrixNZ.hh.

197  {
198  return (first.x * second.x +
199  first.y * second.y +
200  first.v * second.v +
201  first.w * second.w);
202 
203  }

◆ invert()

void JMATH::JMatrixNS::invert ( )
inlineinherited

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 
169  JMatrixND_t& A = getInstance();
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  }
p1
TPaveText * p1
Definition: JDrawModule3D.cc:35
THROW
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:712
root
Definition: root.py:1
JMATH::JMatrixND_t
Basic NxN matrix.
Definition: JMatrixND.hh:36
JMATH::JMatrixND_t::reset
JMatrixND_t & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:130
JMATH::JMatrixND_t::size
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:202
JMATH::JMatrixND_t::swap
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:168
JMATH::JMatrixND_t::transpose
JMatrixND_t & transpose()
Transpose.
Definition: JMatrixND.hh:183
JMATH::JMatrixND::getInstance
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:361
JMATH::JMatrixND::rswap
void rswap(size_t r1, size_t r2)
Definition: JMatrixND.hh:873
JMATH::JMatrixNS::P
std::vector< int > P
permutation matrix
Definition: JMatrixNS.hh:486

◆ solve()

template<class JVectorND_t >
void JMATH::JMatrixNS::solve ( JVectorND_t &  u)
inlineinherited

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  }
JTOOLS::u
double u[N+1]
Definition: JPolint.hh:865

◆ update()

void JMATH::JMatrixNS::update ( const size_t  k,
const double  value 
)
inlineinherited

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  }

◆ getInstance()

JMatrixND_t& JMATH::JMatrixND::getInstance ( )
inlineprotectedinherited

Get work space.

Returns
work space

Definition at line 361 of file JMatrixND.hh.

362  {
363  return ws;
364  }
JMATH::JMatrixND::ws
JMatrixND_t ws
Definition: JMatrixND.hh:354

◆ resize()

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 446 of file JMatrixND.hh.

447  {
448  static_cast<JMatrixND_t&>(*this).resize(size);
449 
450  getInstance().resize(size);
451  }
JMATH::JMatrixND_t::resize
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:102

◆ reset()

JMatrixND& JMATH::JMatrixND::reset ( )
inlineinherited

Set matrix to the null matrix.

Returns
this matrix

Definition at line 459 of file JMatrixND.hh.

460  {
461  JMatrixND_t::reset();
462 
463  JMatrixND_t& A = getInstance();
464 
465  A.resize(this->size());
466 
467  return *this;
468  }

◆ setIdentity()

JMatrixND& JMATH::JMatrixND::setIdentity ( )
inlineinherited

Set to identity matrix.

Returns
this matrix

Definition at line 476 of file JMatrixND.hh.

477  {
478  reset();
479 
480  for (size_t i = 0; i != this->size(); ++i) {
481  (*this)(i,i) = 1.0;
482  }
483 
484  return *this;
485  }
JMATH::JMatrixND::reset
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:459

◆ negate()

JMatrixND& JMATH::JMatrixND::negate ( )
inlineinherited

Negate matrix.

Returns
this matrix

Definition at line 493 of file JMatrixND.hh.

494  {
495  double* p = this->data();
496 
497  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
498  *p = -*p;
499  }
500 
501  return *this;
502  }
JMATH::JMatrixND_t::data
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:235

◆ add()

JMatrixND& JMATH::JMatrixND::add ( const JMatrixND &  A)
inlineinherited

Matrix addition.

Parameters
Amatrix
Returns
this matrix

Definition at line 511 of file JMatrixND.hh.

512  {
513  if (this->size() == A.size()) {
514 
515  double* p = this->data();
516  const double* q = A.data();
517 
518  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
519  *p += *q;
520  }
521 
522  } else {
523  THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
524  }
525 
526  return *this;
527  }

◆ sub()

JMatrixND& JMATH::JMatrixND::sub ( const JMatrixND &  A)
inlineinherited

Matrix subtraction.

Parameters
Amatrix
Returns
this matrix

Definition at line 536 of file JMatrixND.hh.

537  {
538  if (this->size() == A.size()) {
539 
540  double* p = this->data();
541  const double* q = A.data();
542 
543  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
544  *p -= *q;
545  }
546 
547  } else {
548  THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
549  }
550 
551  return *this;
552  }

◆ mul() [1/3]

JMatrixND& JMATH::JMatrixND::mul ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 561 of file JMatrixND.hh.

562  {
563  double* p = this->data();
564 
565  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
566  *p *= factor;
567  }
568 
569  return *this;
570  }

◆ mul() [2/3]

JMatrixND& JMATH::JMatrixND::mul ( const JMatrixND &  A,
const JMatrixND &  B 
)
inlineinherited

Matrix multiplication.

Parameters
Amatrix
Bmatrix
Returns
this matrix

Definition at line 598 of file JMatrixND.hh.

600  {
601  if (A.size() == B.size()) {
602 
603  this->resize(A.size());
604 
605  if (!this->empty()) {
606 
607  JMatrixND_t& C = getInstance();
608 
609  C.set(B);
610  C.transpose();
611 
612  size_t i, row;
613 
614  for (row = 0; row + 4 <= this->size(); row += 4) { // process rows by four
615 
616  double* p0 = (*this)[row + 0];
617  double* p1 = (*this)[row + 1];
618  double* p2 = (*this)[row + 2];
619  double* p3 = (*this)[row + 3];
620 
621  for (size_t col = 0; col != this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
622 
623  double w0 = 0.0;
624  double w1 = 0.0;
625  double w2 = 0.0;
626  double w3 = 0.0;
627 
628  const double* a0 = A[row + 0];
629  const double* a1 = A[row + 1];
630  const double* a2 = A[row + 2];
631  const double* a3 = A[row + 3];
632 
633  const double* c0 = C[col];
634 
635  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4, c0 += 4) {
636  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
637  w1 += a1[0] * c0[0] + a1[1] * c0[1] + a1[2] * c0[2] + a1[3] * c0[3];
638  w2 += a2[0] * c0[0] + a2[1] * c0[1] + a2[2] * c0[2] + a2[3] * c0[3];
639  w3 += a3[0] * c0[0] + a3[1] * c0[1] + a3[2] * c0[2] + a3[3] * c0[3];
640  }
641 
642  for ( ; i != this->size(); ++i, ++a0, ++a1, ++a2, ++a3, ++c0) {
643  w0 += (*a0) * (*c0);
644  w1 += (*a1) * (*c0);
645  w2 += (*a2) * (*c0);
646  w3 += (*a3) * (*c0);
647  }
648 
649  *p0 = w0;
650  *p1 = w1;
651  *p2 = w2;
652  *p3 = w3;
653  }
654  }
655 
656  for ( ; row != this->size(); ++row) { // remaining rows
657 
658  double* p0 = (*this)[row + 0];
659 
660  for (size_t col = 0; col != this->size(); ++col, ++p0) {
661 
662  double w0 = 0.0;
663 
664  const double* a0 = A[row + 0];
665  const double* c0 = C[col];
666 
667  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, c0 += 4) {
668  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
669  }
670 
671  for ( ; i != this->size(); ++i, ++a0, ++c0) {
672  w0 += (*a0) * (*c0);
673  }
674 
675  *p0 = w0;
676  }
677  }
678  }
679 
680  } else {
681  THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
682  }
683 
684  return *this;
685  }
JPHYSICS::C
static const double C
Physics constants.
Definition: JPhysics/JConstants.hh:25
JMATH::JMatrixND_t::empty
bool empty() const
Check emptyness.
Definition: JMatrixND.hh:224
JMATH::JMatrixND::resize
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:446

◆ mul() [3/3]

JMatrixND & JMATH::JMath< JMatrixND , JSecond_t >::mul ( const JSecond_t &  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  }

◆ div()

JMatrixND& JMATH::JMatrixND::div ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 579 of file JMatrixND.hh.

580  {
581  double* p = this->data();
582 
583  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
584  *p /= factor;
585  }
586 
587  return *this;
588  }

◆ equals()

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 695 of file JMatrixND.hh.

697  {
698  if (this->size() == A.size()) {
699 
700  for (size_t row = 0; row != this->size(); ++row) {
701 
702  const double* p = (*this)[row];
703  const double* q = A [row];
704 
705  for (size_t i = this->size(); i != 0; --i, ++p, ++q) {
706  if (fabs(*p - *q) > eps) {
707  return false;
708  }
709  }
710  }
711 
712  return true;
713 
714  } else {
715  THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
716  }
717  }

◆ isIdentity()

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 726 of file JMatrixND.hh.

727  {
728  for (size_t i = 0; i != this->size(); ++i) {
729 
730  if (fabs(1.0 - (*this)(i,i)) > eps) {
731  return false;
732  };
733 
734  for (size_t j = 0; j != i; ++j) {
735  if (fabs((*this)(i,j)) > eps || fabs((*this)(j,i)) > eps) {
736  return false;
737  };
738  }
739  }
740 
741  return true;
742  }

◆ isSymmetric()

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 751 of file JMatrixND.hh.

752  {
753  for (size_t i = 0; i != this->size(); ++i) {
754  for (size_t j = 0; j != i; ++j) {
755  if (fabs((*this)(i,j) - (*this)(j,i)) > eps) {
756  return false;
757  };
758  }
759  }
760 
761  return true;
762  }

◆ getDot() [2/2]

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 773 of file JMatrixND.hh.

774  {
775  double dot = 0.0;
776 
777  for (size_t i = 0; i != v.size(); ++i) {
778 
779  const double* p = (*this)[i];
780 
781  double w = 0.0;
782 
783  for (JVectorND::const_iterator y = v.begin(); y != v.end(); ++y, ++p) {
784  w += (*p) * (*y);
785  }
786 
787  dot += v[i] * w;
788  }
789 
790  return dot;
791  }

◆ rswap()

void JMATH::JMatrixND::rswap ( size_t  r1,
size_t  r2 
)
inlineprotectedinherited

Definition at line 873 of file JMatrixND.hh.

874  {
875  JMatrixND_t& A = getInstance();
876 
877  A.resize(this->size());
878 
879  memcpy(A.data(), (*this)[r1], this->size() * sizeof(double));
880  memcpy((*this)[r1], (*this)[r2], this->size() * sizeof(double));
881  memcpy((*this)[r2], A.data(), this->size() * sizeof(double));
882  }

◆ cswap()

void JMATH::JMatrixND::cswap ( size_t  c1,
size_t  c2 
)
inlineprotectedinherited

Swap columns.

Parameters
c1column
c2column

Definition at line 891 of file JMatrixND.hh.

892  {
893  using std::swap;
894 
895  double* p1 = this->data() + c1;
896  double* p2 = this->data() + c2;
897 
898  for (size_t i = this->size(); i != 0; --i, p1 += this->size(), p2 += this->size()) {
899  swap(*p1, *p2);
900  }
901  }
c1
TCanvas * c1
Global variables to handle mouse events.
Definition: JDrawModule3D.cc:34

◆ clear()

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  }
JMATH::JMatrixND_t::__p
double * __p
pointer to data
Definition: JMatrixND.hh:339
JMATH::JMatrixND_t::__m
size_t __m
capacity of matrix
Definition: JMatrixND.hh:341
JMATH::JMatrixND_t::__n
size_t __n
dimension of matrix
Definition: JMatrixND.hh:340

◆ set() [2/2]

void JMATH::JMatrixND_t::set ( const JMatrixND_t &  A)
inlineinherited

Set matrix.

Parameters
Amatrix

Definition at line 155 of file JMatrixND.hh.

156  {
157  this->resize(A.size());
158 
159  memcpy(this->data(), A.data(), A.size() * A.size() * sizeof(double));
160  }

◆ swap()

void JMATH::JMatrixND_t::swap ( JMatrixND_t &  A)
inlineinherited

Swap matrices.

Parameters
Amatrix

Definition at line 168 of file JMatrixND.hh.

169  {
170  using std::swap;
171 
172  swap(this->__p, A.__p);
173  swap(this->__n, A.__n);
174  swap(this->__m, A.__m);
175  }

◆ transpose()

JMatrixND_t& JMATH::JMatrixND_t::transpose ( )
inlineinherited

Transpose.

Returns
this matrix

Definition at line 183 of file JMatrixND.hh.

184  {
185  using std::swap;
186 
187  for (size_t row = 0; row != this->size(); ++row) {
188  for (size_t col = 0; col != row; ++col) {
189  swap((*this)(row,col), (*this)(col,row));
190  }
191  }
192 
193  return *this;
194  }

◆ size()

size_t JMATH::JMatrixND_t::size ( ) const
inlineinherited

Get dimension of matrix.

Returns
dimension

Definition at line 202 of file JMatrixND.hh.

203  {
204  return __n;
205  }

◆ capacity()

size_t JMATH::JMatrixND_t::capacity ( ) const
inlineinherited

Get capacity of dimension.

Returns
capacity

Definition at line 213 of file JMatrixND.hh.

214  {
215  return __m;
216  }

◆ empty()

bool JMATH::JMatrixND_t::empty ( ) const
inlineinherited

Check emptyness.

Returns
true if empty; else false

Definition at line 224 of file JMatrixND.hh.

225  {
226  return __n == 0;
227  }

◆ data() [1/2]

const double* JMATH::JMatrixND_t::data ( ) const
inlineinherited

Get pointer to data.

Returns
pointer to data.

Definition at line 235 of file JMatrixND.hh.

236  {
237  return __p;
238  }

◆ data() [2/2]

double* JMATH::JMatrixND_t::data ( )
inlineinherited

Get pointer to data.

Returns
pointer to data.

Definition at line 246 of file JMatrixND.hh.

247  {
248  return __p;
249  }

◆ operator[]() [1/2]

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 258 of file JMatrixND.hh.

259  {
260  return __p + row*__n;
261  }

◆ operator[]() [2/2]

double* JMATH::JMatrixND_t::operator[] ( size_t  row)
inlineinherited

Get row data.

Parameters
rowrow number
Returns
pointer to row data

Definition at line 270 of file JMatrixND.hh.

271  {
272  return __p + row*__n;
273  }

◆ operator()() [1/2]

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 283 of file JMatrixND.hh.

284  {
285  return *(__p + row*__n + col);
286  }

◆ operator()() [2/2]

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 296 of file JMatrixND.hh.

297  {
298  return *(__p + row*__n + col);
299  }

◆ at() [1/2]

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 309 of file JMatrixND.hh.

310  {
311  if (row >= this->size() ||
312  col >= this->size()) {
313  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
314  }
315 
316  return (*this)(row, col);
317  }

◆ at() [2/2]

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 327 of file JMatrixND.hh.

328  {
329  if (row >= this->size() ||
330  col >= this->size()) {
331  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
332  }
333 
334  return (*this)(row, col);
335  }

Member Data Documentation

◆ buffer

std::vector<variance> JFIT::JMatrixNZ::buffer
private

Definition at line 186 of file JMatrixNZ.hh.

◆ P

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

permutation matrix

Definition at line 486 of file JMatrixNS.hh.

◆ ws

JMatrixND_t JMATH::JMatrixND::ws
protectedinherited

Definition at line 354 of file JMatrixND.hh.

◆ __p

double* JMATH::JMatrixND_t::__p
protectedinherited

pointer to data

Definition at line 339 of file JMatrixND.hh.

◆ __n

size_t JMATH::JMatrixND_t::__n
protectedinherited

dimension of matrix

Definition at line 340 of file JMatrixND.hh.

◆ __m

size_t JMATH::JMatrixND_t::__m
protectedinherited

capacity of matrix

Definition at line 341 of file JMatrixND.hh.

The documentation for this class was generated from the following file:
Generated by doxygen 1.9.1