Jpp  17.2.0
the software that should make you happy
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Classes | Public Types | 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 Types

typedef std::vector< std::pair
< size_t, size_t > > 
JPermutationMatrix
 Type definition of permutation matrix. 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...
 
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...
 

Static Private Member Functions

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

Private Attributes

std::vector< variancebuffer
 

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.

Member Typedef Documentation

typedef std::vector< std::pair<size_t,size_t> > JMATH::JMatrixNS::JPermutationMatrix
inherited

Type definition of permutation matrix.

Definition at line 33 of file JMatrixNS.hh.

Constructor & Destructor Documentation

JFIT::JMatrixNZ::JMatrixNZ ( )
inline

Default contructor.

Definition at line 39 of file JMatrixNZ.hh.

39  :
40  JMatrixNS()
41  {}
JMatrixNS()
Default constructor.
Definition: JMatrixNS.hh:39
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  }
void set(const JVector3D &pos, T __begin, T __end, const double alpha, const double sigma)
Set co-variance matrix.
Definition: JMatrixNZ.hh:85
JMatrixNS()
Default constructor.
Definition: JMatrixNS.hh:39
const double sigma[]
Definition: JQuadrature.cc:74

Member Function Documentation

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  }
data_type w[N+1][M+1]
Definition: JPolint.hh:778
double getKappaC()
Get average R-dependence of arrival time of Cherenkov light (a.k.a.
std::vector< T >::difference_type distance(typename std::vector< T >::const_iterator first, typename PhysicsEvent::const_iterator< T > second)
Specialisation of STL distance.
static double getDot(const variance &first, const variance &second)
Get dot product.
Definition: JMatrixNZ.hh:196
then JShowerPostfit f $INPUT_FILE o $OUTPUT_FILE N
const double sigma[]
Definition: JQuadrature.cc:74
do set_variable OUTPUT_DIRECTORY $WORKDIR T
std::vector< variance > buffer
Definition: JMatrixNZ.hh:186
static const double PI
Mathematical constants.
then usage $script[distance] fi case set_variable R
Definition: JDrawLED.sh:43
const double getInverseSpeedOfLight()
Get inverse speed of light.
int j
Definition: JPolint.hh:703
data_type v[N+1][M+1]
Definition: JPolint.hh:777
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  }
then echo The file $DIR KM3NeT_00000001_00000000 root already please rename or remove it first
void JMATH::JMatrixNS::invert ( )
inlineinherited

Invert matrix according LDU decomposition.

Definition at line 80 of file JMatrixNS.hh.

81  {
82  using namespace std;
83 
84  size_t i, row, col, root;
85 
86  double* p0;
87  double* p1;
88  double* p2;
89  double* p3;
90 
91  const double* c0;
92  const double* c1;
93  const double* c2;
94  const double* c3;
95 
96  const double* a0;
97  const double* a1;
98  const double* a2;
99  const double* a3;
100 
101  double v0;
102  double v1;
103  double v2;
104  double v3;
105 
106  double w0;
107  double w1;
108  double w2;
109  double w3;
110 
112 
113 
114  // LDU decomposition
115 
116  for (root = 0; root != this->size(); ++root) {
117 
118  double value = (*this)(root, root);
119  size_t pivot = root;
120 
121  for (row = root; ++row != this->size(); ) {
122  if (fabs((*this)(row, root)) > fabs(value)) {
123  value = (*this)(row, root);
124  pivot = row;
125  }
126  }
127 
128  if (value == 0.0) {
129  THROW(JDivisionByZero, "LDU decomposition zero pivot at " << root);
130  }
131 
132  if (pivot != root) {
133 
134  this->rswap(root, pivot);
135 
136  P.push_back(make_pair(root, pivot));
137  }
138 
139  for (row = root + 1; row + 4 <= this->size(); row += 4) { // process rows by four
140 
141  p0 = (*this)[row + 0] + root;
142  p1 = (*this)[row + 1] + root;
143  p2 = (*this)[row + 2] + root;
144  p3 = (*this)[row + 3] + root;
145 
146  v0 = (*p0++ /= value);
147  v1 = (*p1++ /= value);
148  v2 = (*p2++ /= value);
149  v3 = (*p3++ /= value);
150 
151  c0 = (*this)[root] + root + 1;
152 
153  for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, p1 += 4, p2 += 4, p3 += 4, c0 += 4) {
154  p0[0] -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
155  p1[0] -= v1 * c0[0]; p1[1] -= v1 * c0[1]; p1[2] -= v1 * c0[2]; p1[3] -= v1 * c0[3];
156  p2[0] -= v2 * c0[0]; p2[1] -= v2 * c0[1]; p2[2] -= v2 * c0[2]; p2[3] -= v2 * c0[3];
157  p3[0] -= v3 * c0[0]; p3[1] -= v3 * c0[1]; p3[2] -= v3 * c0[2]; p3[3] -= v3 * c0[3];
158  }
159 
160  for ( ; i != this->size(); ++i, ++p0, ++p1, ++p2, ++p3, ++c0) {
161  *p0 -= v0 * (*c0);
162  *p1 -= v1 * (*c0);
163  *p2 -= v2 * (*c0);
164  *p3 -= v3 * (*c0);
165  }
166  }
167 
168  for ( ; row != this->size(); ++row) { // remaining rows
169 
170  p0 = (*this)[row + 0] + root;
171 
172  v0 = (*p0++ /= value);
173 
174  c0 = (*this)[root] + root + 1;
175 
176  for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, c0 += 4) {
177  *p0 -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
178  }
179 
180  for ( ; i != this->size(); ++i, ++p0, ++c0) {
181  *p0 -= v0 * (*c0);
182  }
183  }
184  }
185 
186 
187  // invert D
188 
189  for (size_t i = 0; i != this->size(); ++i) {
190 
191  double& pivot = (*this)(i,i);
192 
193  if (pivot == 0.0) {
194  THROW(JDivisionByZero, "Invert D zero pivot at " << i);
195  }
196 
197  pivot = 1.0 / pivot;
198  }
199 
200 
201  // invert U
202 
204 
205  A.set(*this);
206  A.transpose();
207 
208  for (row = 0; row + 4 <= this->size(); row += 4) { // process rows by four
209 
210  p0 = (*this)[row + 0] + row + 0;
211  p1 = (*this)[row + 1] + row + 1;
212  p2 = (*this)[row + 2] + row + 2;
213  p3 = (*this)[row + 3] + row + 3;
214 
215  v0 = *(p0++);
216  v1 = *(p1++);
217  v2 = *(p2++);
218  v3 = *(p3++);
219 
220  for (col = row + 1; col + 4 <= this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
221 
222  w0 = (*p0) * -v0;
223  w1 = (*p1) * -v1;
224  w2 = (*p2) * -v2;
225  w3 = (*p3) * -v3;
226 
227  c0 = (*this)[row + 0] + row + 0 + 1;
228  c1 = (*this)[row + 1] + row + 1 + 1;
229  c2 = (*this)[row + 2] + row + 2 + 1;
230  c3 = (*this)[row + 3] + row + 3 + 1;
231 
232  a0 = A[col + 0] + row + 0 + 1;
233  a1 = A[col + 1] + row + 1 + 1;
234  a2 = A[col + 2] + row + 2 + 1;
235  a3 = A[col + 3] + row + 3 + 1;
236 
237  for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
238  w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
239  w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3];
240  w2 -= c2[0]*a2[0] + c2[1]*a2[1] + c2[2]*a2[2] + c2[3]*a2[3];
241  w3 -= c3[0]*a3[0] + c3[1]*a3[1] + c3[2]*a3[2] + c3[3]*a3[3];
242  }
243 
244  for ( ; c0 != p0; ++c0, ++c1, ++c2, ++c3, ++a0, ++a1, ++a2, ++a3) {
245  w0 -= (*c0) * (*a0);
246  w1 -= (*c1) * (*a1);
247  w2 -= (*c2) * (*a2);
248  w3 -= (*c3) * (*a3);
249  }
250 
251  *p0 = w0 * (*this)(col + 0, col + 0);
252  *p1 = w1 * (*this)(col + 1, col + 1);
253  *p2 = w2 * (*this)(col + 2, col + 2);
254  *p3 = w3 * (*this)(col + 3, col + 3);
255  }
256 
257  // last 3 columns
258 
259  w0 = (*p0) * -v0;
260  w1 = (*p1) * -v1;
261  w2 = (*p2) * -v2;
262 
263  c0 = (*this)[row + 0] + row + 0 + 1;
264  c1 = (*this)[row + 1] + row + 1 + 1;
265  c2 = (*this)[row + 2] + row + 2 + 1;
266 
267  a0 = A[col + 0] + row + 0 + 1;
268  a1 = A[col + 1] + row + 1 + 1;
269  a2 = A[col + 2] + row + 2 + 1;
270 
271  for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
272  w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
273  w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3];
274  w2 -= c2[0]*a2[0] + c2[1]*a2[1] + c2[2]*a2[2] + c2[3]*a2[3];
275  }
276 
277  for ( ; c0 != p0; ++c0, ++c1, ++c2, ++a0, ++a1, ++a2) {
278  w0 -= (*c0) * (*a0);
279  w1 -= (*c1) * (*a1);
280  w2 -= (*c2) * (*a2);
281  }
282 
283  *p0 = w0 * (*this)(col + 0, col + 0);
284  *p1 = w1 * (*this)(col + 1, col + 1);
285  *p2 = w2 * (*this)(col + 2, col + 2);
286 
287  ++col; ++p0; ++p1;
288 
289  // last 2 columns
290 
291  w0 = (*p0) * -v0;
292  w1 = (*p1) * -v1;
293 
294  c0 = (*this)[row + 0] + row + 0 + 1;
295  c1 = (*this)[row + 1] + row + 1 + 1;
296 
297  a0 = A[col + 0] + row + 0 + 1;
298  a1 = A[col + 1] + row + 1 + 1;
299 
300  for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
301  w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
302  w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3];
303  }
304 
305  for ( ; c0 != p0; ++c0, ++c1, ++a0, ++a1) {
306  w0 -= (*c0) * (*a0);
307  w1 -= (*c1) * (*a1);
308  }
309 
310  *p0 = w0 * (*this)(col + 0, col + 0);
311  *p1 = w1 * (*this)(col + 1, col + 1);
312 
313  ++col; ++p0;
314 
315  // last column
316 
317  w0 = (*p0) * -v0;
318 
319  c0 = (*this)[row + 0] + row + 0 + 1;
320 
321  a0 = A[col + 0] + row + 0 + 1;
322 
323  for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
324  w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
325  }
326 
327  for ( ; c0 != p0; ++c0, ++a0) {
328  w0 -= (*c0) * (*a0);
329  }
330 
331  *p0 = w0 * (*this)(col + 0, col + 0);
332  }
333 
334  for ( ; row != this->size(); ++row) { // remaining rows
335 
336  p0 = (*this)[row + 0] + row + 0;
337 
338  v0 = *(p0++);
339 
340  for (col = row + 1; col != this->size(); ++col, ++p0) {
341 
342  w0 = *p0;
343 
344  w0 *= -v0;
345 
346  c0 = (*this)[row + 0] + row + 1;
347 
348  a0 = A[col + 0] + row + 0 + 1;
349 
350  for ( ; c0 != p0; ++c0, ++a0) {
351  w0 -= (*c0) * (*a0);
352  }
353 
354  *p0 = w0 * (*this)(col + 0, col + 0);
355  }
356  }
357 
358 
359  // solve A^-1 x L = U^-1
360 
361  double* buffer = A.data();
362 
363  for (col = this->size() - 1; col != 0; --col) {
364 
365  for (i = col; i < this->size(); ++i) {
366 
367  buffer[i] = (*this)(i, col - 1);
368 
369  (*this)(i, col - 1) = 0.0;
370  }
371 
372  for (row = 0; row + 4 < this->size(); row += 4) {
373 
374  p0 = (*this)[row + 0] + col - 1;
375  p1 = (*this)[row + 1] + col - 1;
376  p2 = (*this)[row + 2] + col - 1;
377  p3 = (*this)[row + 3] + col - 1;
378 
379  c0 = p0 + 1;
380  c1 = p1 + 1;
381  c2 = p2 + 1;
382  c3 = p3 + 1;
383 
384  a0 = buffer + col;
385 
386  w0 = 0.0;
387  w1 = 0.0;
388  w2 = 0.0;
389  w3 = 0.0;
390 
391  for (i = col; i + 4 <= this->size(); i += 4, c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4) {
392  w0 += c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
393  w1 += c1[0]*a0[0] + c1[1]*a0[1] + c1[2]*a0[2] + c1[3]*a0[3];
394  w2 += c2[0]*a0[0] + c2[1]*a0[1] + c2[2]*a0[2] + c2[3]*a0[3];
395  w3 += c3[0]*a0[0] + c3[1]*a0[1] + c3[2]*a0[2] + c3[3]*a0[3];
396  }
397 
398  for ( ; i != this->size(); ++i, ++c0, ++c1, ++c2, ++c3, ++a0) {
399  w0 += (*c0) * (*a0);
400  w1 += (*c1) * (*a0);
401  w2 += (*c2) * (*a0);
402  w3 += (*c3) * (*a0);
403  }
404 
405  *p0 -= w0;
406  *p1 -= w1;
407  *p2 -= w2;
408  *p3 -= w3;
409  }
410 
411  for ( ; row < this->size(); ++row) {
412 
413  p0 = (*this)[row] + col - 1;
414 
415  c0 = p0 + 1;
416 
417  a0 = buffer + col;
418 
419  w0 = 0.0;
420 
421  for (i = col; i + 4 <= this->size(); i += 4, c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4) {
422  w0 += c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
423  }
424 
425  for ( ; i != this->size(); ++i, ++c0, ++a0) {
426  w0 += (*c0) * (*a0);
427  }
428 
429  *p0 -= w0;
430  }
431  }
432 
433 
434  for (JPermutationMatrix::reverse_iterator p = P.rbegin(); p != P.rend(); ++p) {
435  cswap(p->first, p->second);
436  }
437  }
TPaveText * p1
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:336
void rswap(size_t r1, size_t r2)
Definition: JMatrixND.hh:856
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
TCanvas * c1
Global variables to handle mouse events.
p2
Definition: module-Z:fit.sh:74
do JPlot2D f $WORKDIR detector root
std::vector< std::pair< size_t, size_t > > JPermutationMatrix
Type definition of permutation matrix.
Definition: JMatrixNS.hh:33
void cswap(size_t c1, size_t c2)
Swap columns.
Definition: JMatrixND.hh:874
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
then $DIR JPlotNPE PDG P
Definition: JPlotNPE-PDG.sh:62
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 457 of file JMatrixNS.hh.

458  {
459  if (k >= this->size()) {
460  THROW(JIndexOutOfRange, "JMatrixNS::update(): index out of range " << k);
461  }
462 
463  size_t i, row;
464  double* col;
465 
466  const double del = (*this)(k,k);
467  const double din = 1.0 / (1.0 + value * del);
468 
469  double* root = (*this)[k];
470 
471  for (row = 0; row != this->size(); ++row) {
472 
473  if (row != k) {
474 
475  const double val = (*this)(row, k);
476 
477  for (i = 0, col = (*this)[row]; i != this->size(); ++i, ++col) {
478  if (i != k) {
479  (*col) -= val * root[i] * value * din;
480  }
481  }
482 
483  (*this)(row, k) *= din;
484  }
485  }
486 
487  for (i = 0, col = root; i != this->size(); ++i, ++col) {
488  if (i != k) {
489  (*col) *= din;
490  }
491  }
492 
493  root[k] = del * din;
494  }
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:696
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
do JPlot2D f $WORKDIR detector root
JMatrixND_t& JMATH::JMatrixND::getInstance ( )
inlineprotectedinherited

Get work space.

Returns
work space

Definition at line 336 of file JMatrixND.hh.

337  {
338  return ws;
339  }
JMatrixND_t ws
Definition: JMatrixND.hh:329
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 421 of file JMatrixND.hh.

422  {
423  static_cast<JMatrixND_t&>(*this).resize(size);
424 
426  }
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:336
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
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 434 of file JMatrixND.hh.

435  {
437 
438  A.resize(this->size());
439 
440  double* p = A.data();
441 
442  for (size_t i = this->size(); i != 0; --i, ++p) {
443  *p = 0.0;
444  }
445 
446  p = this->data();
447 
448  for (size_t i = this->size(); i != 0; --i, p += this->size()) {
449  memcpy(p, A.data(), this->size() * sizeof(double));
450  }
451 
452  return *this;
453  }
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:336
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
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:210
Basic NxN matrix.
Definition: JMatrixND.hh:36
JMatrixND& JMATH::JMatrixND::setIdentity ( )
inlineinherited

Set to identity matrix.

Returns
this matrix

Definition at line 461 of file JMatrixND.hh.

462  {
463  reset();
464 
465  for (size_t i = 0; i != this->size(); ++i) {
466  (*this)(i,i) = 1.0;
467  }
468 
469  return *this;
470  }
JMatrixND & reset()
Set matrix to the null matrix.
Definition: JMatrixND.hh:434
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
JMatrixND& JMATH::JMatrixND::negate ( )
inlineinherited

Negate matrix.

Returns
this matrix

Definition at line 478 of file JMatrixND.hh.

479  {
480  double* p = this->data();
481 
482  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
483  *p = -*p;
484  }
485 
486  return *this;
487  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:210
JMatrixND& JMATH::JMatrixND::add ( const JMatrixND A)
inlineinherited

Matrix addition.

Parameters
Amatrix
Returns
this matrix

Definition at line 496 of file JMatrixND.hh.

497  {
498  if (this->size() == A.size()) {
499 
500  double* p = this->data();
501  const double* q = A.data();
502 
503  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
504  *p += *q;
505  }
506 
507  } else {
508  THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
509  }
510  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:210
JMatrixND& JMATH::JMatrixND::sub ( const JMatrixND A)
inlineinherited

Matrix subtraction.

Parameters
Amatrix
Returns
this matrix

Definition at line 519 of file JMatrixND.hh.

520  {
521  if (this->size() == A.size()) {
522 
523  double* p = this->data();
524  const double* q = A.data();
525 
526  for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
527  *p -= *q;
528  }
529 
530  } else {
531  THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
532  }
533 
534  return *this;
535  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:210
JMatrixND& JMATH::JMatrixND::mul ( const double  factor)
inlineinherited

Scale matrix.

Parameters
factorfactor
Returns
this matrix

Definition at line 544 of file JMatrixND.hh.

545  {
546  double* p = this->data();
547 
548  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
549  *p *= factor;
550  }
551 
552  return *this;
553  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:210
JMatrixND& JMATH::JMatrixND::mul ( const JMatrixND A,
const JMatrixND B 
)
inlineinherited

Matrix multiplication.

Parameters
Amatrix
Bmatrix
Returns
this matrix

Definition at line 581 of file JMatrixND.hh.

583  {
584  if (A.size() == B.size()) {
585 
586  this->resize(A.size());
587 
588  if (!this->empty()) {
589 
591 
592  C.set(B);
593  C.transpose();
594 
595  size_t i, row;
596 
597  for (row = 0; row + 4 <= this->size(); row += 4) { // process rows by four
598 
599  double* p0 = (*this)[row + 0];
600  double* p1 = (*this)[row + 1];
601  double* p2 = (*this)[row + 2];
602  double* p3 = (*this)[row + 3];
603 
604  for (size_t col = 0; col != this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
605 
606  double w0 = 0.0;
607  double w1 = 0.0;
608  double w2 = 0.0;
609  double w3 = 0.0;
610 
611  const double* a0 = A[row + 0];
612  const double* a1 = A[row + 1];
613  const double* a2 = A[row + 2];
614  const double* a3 = A[row + 3];
615 
616  const double* c0 = C[col];
617 
618  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4, c0 += 4) {
619  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
620  w1 += a1[0] * c0[0] + a1[1] * c0[1] + a1[2] * c0[2] + a1[3] * c0[3];
621  w2 += a2[0] * c0[0] + a2[1] * c0[1] + a2[2] * c0[2] + a2[3] * c0[3];
622  w3 += a3[0] * c0[0] + a3[1] * c0[1] + a3[2] * c0[2] + a3[3] * c0[3];
623  }
624 
625  for ( ; i != this->size(); ++i, ++a0, ++a1, ++a2, ++a3, ++c0) {
626  w0 += (*a0) * (*c0);
627  w1 += (*a1) * (*c0);
628  w2 += (*a2) * (*c0);
629  w3 += (*a3) * (*c0);
630  }
631 
632  *p0 = w0;
633  *p1 = w1;
634  *p2 = w2;
635  *p3 = w3;
636  }
637  }
638 
639  for ( ; row != this->size(); ++row) { // remaining rows
640 
641  double* p0 = (*this)[row + 0];
642 
643  for (size_t col = 0; col != this->size(); ++col, ++p0) {
644 
645  double w0 = 0.0;
646 
647  const double* a0 = A[row + 0];
648  const double* c0 = C[col];
649 
650  for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, c0 += 4) {
651  w0 += a0[0] * c0[0] + a0[1] * c0[1] + a0[2] * c0[2] + a0[3] * c0[3];
652  }
653 
654  for ( ; i != this->size(); ++i, ++a0, ++c0) {
655  w0 += (*a0) * (*c0);
656  }
657 
658  *p0 = w0;
659  }
660  }
661  }
662 
663  } else {
664  THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
665  }
666 
667  return *this;
668  }
TPaveText * p1
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:421
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:336
static const double C
Physics constants.
void set(const JMatrixND_t &A)
Set matrix.
Definition: JMatrixND.hh:130
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
p2
Definition: module-Z:fit.sh:74
JMatrixND_t & transpose()
Transpose.
Definition: JMatrixND.hh:158
bool empty() const
Check emptyness.
Definition: JMatrixND.hh:199
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 562 of file JMatrixND.hh.

563  {
564  double* p = this->data();
565 
566  for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
567  *p /= factor;
568  }
569 
570  return *this;
571  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:210
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 678 of file JMatrixND.hh.

680  {
681  if (this->size() == A.size()) {
682 
683  for (size_t row = 0; row != this->size(); ++row) {
684 
685  const double* p = (*this)[row];
686  const double* q = A [row];
687 
688  for (size_t i = this->size(); i != 0; --i, ++p, ++q) {
689  if (fabs(*p - *q) > eps) {
690  return false;
691  }
692  }
693  }
694 
695  return true;
696 
697  } else {
698  THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
699  }
700  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
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 709 of file JMatrixND.hh.

710  {
711  for (size_t i = 0; i != this->size(); ++i) {
712 
713  if (fabs(1.0 - (*this)(i,i)) > eps) {
714  return false;
715  };
716 
717  for (size_t j = 0; j != i; ++j) {
718  if (fabs((*this)(i,j)) > eps || fabs((*this)(j,i)) > eps) {
719  return false;
720  };
721  }
722  }
723 
724  return true;
725  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
int j
Definition: JPolint.hh:703
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 734 of file JMatrixND.hh.

735  {
736  for (size_t i = 0; i != this->size(); ++i) {
737  for (size_t j = 0; j != i; ++j) {
738  if (fabs((*this)(i,j) - (*this)(j,i)) > eps) {
739  return false;
740  };
741  }
742  }
743 
744  return true;
745  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
int j
Definition: JPolint.hh:703
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 756 of file JMatrixND.hh.

757  {
758  double dot = 0.0;
759 
760  for (size_t i = 0; i != v.size(); ++i) {
761 
762  const double* p = (*this)[i];
763 
764  double w = 0.0;
765 
766  for (JVectorND::const_iterator y = v.begin(); y != v.end(); ++y, ++p) {
767  w += (*p) * (*y);
768  }
769 
770  dot += v[i] * w;
771  }
772 
773  return dot;
774  }
data_type w[N+1][M+1]
Definition: JPolint.hh:778
void JMATH::JMatrixND::rswap ( size_t  r1,
size_t  r2 
)
inlineprotectedinherited

Definition at line 856 of file JMatrixND.hh.

857  {
859 
860  A.resize(this->size());
861 
862  memcpy(A.data(), (*this)[r1], this->size() * sizeof(double));
863  memcpy((*this)[r1], (*this)[r2], this->size() * sizeof(double));
864  memcpy((*this)[r2], A.data(), this->size() * sizeof(double));
865  }
JMatrixND_t & getInstance()
Get work space.
Definition: JMatrixND.hh:336
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
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:210
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 874 of file JMatrixND.hh.

875  {
876  using std::swap;
877 
878  double* p1 = this->data() + c1;
879  double* p2 = this->data() + c2;
880 
881  for (size_t i = this->size(); i != 0; --i, p1 += this->size(), p2 += this->size()) {
882  swap(*p1, *p2);
883  }
884  }
TPaveText * p1
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
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:143
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:210
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:315
double * __p
pointer to data
Definition: JMatrixND.hh:314
size_t __m
capacity of matrix
Definition: JMatrixND.hh:316
void JMATH::JMatrixND_t::set ( const JMatrixND_t A)
inlineinherited

Set matrix.

Parameters
Amatrix

Definition at line 130 of file JMatrixND.hh.

131  {
132  this->resize(A.size());
133 
134  memcpy(this->data(), A.data(), A.size() * A.size() * sizeof(double));
135  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
void resize(const size_t size)
Resize matrix.
Definition: JMatrixND.hh:102
const double * data() const
Get pointer to data.
Definition: JMatrixND.hh:210
void JMATH::JMatrixND_t::swap ( JMatrixND_t A)
inlineinherited

Swap matrices.

Parameters
Amatrix

Definition at line 143 of file JMatrixND.hh.

144  {
145  using std::swap;
146 
147  swap(this->__p, A.__p);
148  swap(this->__n, A.__n);
149  swap(this->__m, A.__m);
150  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:315
double * __p
pointer to data
Definition: JMatrixND.hh:314
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:143
size_t __m
capacity of matrix
Definition: JMatrixND.hh:316
JMatrixND_t& JMATH::JMatrixND_t::transpose ( )
inlineinherited

Transpose.

Returns
this matrix

Definition at line 158 of file JMatrixND.hh.

159  {
160  using std::swap;
161 
162  for (size_t row = 0; row != this->size(); ++row) {
163  for (size_t col = 0; col != row; ++col) {
164  swap((*this)(row,col), (*this)(col,row));
165  }
166  }
167 
168  return *this;
169  }
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
void swap(JMatrixND_t &A)
Swap matrices.
Definition: JMatrixND.hh:143
size_t JMATH::JMatrixND_t::size ( ) const
inlineinherited

Get dimension of matrix.

Returns
dimension

Definition at line 177 of file JMatrixND.hh.

178  {
179  return __n;
180  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:315
size_t JMATH::JMatrixND_t::capacity ( ) const
inlineinherited

Get capacity of dimension.

Returns
capacity

Definition at line 188 of file JMatrixND.hh.

189  {
190  return __m;
191  }
size_t __m
capacity of matrix
Definition: JMatrixND.hh:316
bool JMATH::JMatrixND_t::empty ( ) const
inlineinherited

Check emptyness.

Returns
true if empty; else false

Definition at line 199 of file JMatrixND.hh.

200  {
201  return __n == 0;
202  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:315
const double* JMATH::JMatrixND_t::data ( ) const
inlineinherited

Get pointer to data.

Returns
pointer to data.

Definition at line 210 of file JMatrixND.hh.

211  {
212  return __p;
213  }
double * __p
pointer to data
Definition: JMatrixND.hh:314
double* JMATH::JMatrixND_t::data ( )
inlineinherited

Get pointer to data.

Returns
pointer to data.

Definition at line 221 of file JMatrixND.hh.

222  {
223  return __p;
224  }
double * __p
pointer to data
Definition: JMatrixND.hh:314
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 233 of file JMatrixND.hh.

234  {
235  return __p + row*__n;
236  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:315
double * __p
pointer to data
Definition: JMatrixND.hh:314
double* JMATH::JMatrixND_t::operator[] ( size_t  row)
inlineinherited

Get row data.

Parameters
rowrow number
Returns
pointer to row data

Definition at line 245 of file JMatrixND.hh.

246  {
247  return __p + row*__n;
248  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:315
double * __p
pointer to data
Definition: JMatrixND.hh:314
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 258 of file JMatrixND.hh.

259  {
260  return *(__p + row*__n + col);
261  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:315
double * __p
pointer to data
Definition: JMatrixND.hh:314
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 271 of file JMatrixND.hh.

272  {
273  return *(__p + row*__n + col);
274  }
size_t __n
dimension of matrix
Definition: JMatrixND.hh:315
double * __p
pointer to data
Definition: JMatrixND.hh:314
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 284 of file JMatrixND.hh.

285  {
286  if (row >= this->size() ||
287  col >= this->size()) {
288  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
289  }
290 
291  return (*this)(row, col);
292  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177
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 302 of file JMatrixND.hh.

303  {
304  if (row >= this->size() ||
305  col >= this->size()) {
306  THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
307  }
308 
309  return (*this)(row, col);
310  }
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Definition: JException.hh:696
size_t size() const
Get dimension of matrix.
Definition: JMatrixND.hh:177

Member Data Documentation

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

Definition at line 186 of file JMatrixNZ.hh.

JMatrixND_t JMATH::JMatrixND::ws
protectedinherited

Definition at line 329 of file JMatrixND.hh.

double* JMATH::JMatrixND_t::__p
protectedinherited

pointer to data

Definition at line 314 of file JMatrixND.hh.

size_t JMATH::JMatrixND_t::__n
protectedinherited

dimension of matrix

Definition at line 315 of file JMatrixND.hh.

size_t JMATH::JMatrixND_t::__m
protectedinherited

capacity of matrix

Definition at line 316 of file JMatrixND.hh.


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