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:

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...

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 JNullType &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

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.

                 :
       JMatrixNS()
     {}

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

pos	reference position [m]
__begin	begin of data
__end	end of data
alpha	angular resolution [deg]
sigma	time resolution [ns]

Definition at line 59 of file JMatrixNZ.hh.

                                       :
       JMatrixNS()
     {
       set(pos, __begin, __end, alpha, sigma);
     }

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

pos	reference position [m]
__begin	begin of data
__end	end of data
alpha	angular resolution [deg]
sigma	time resolution [ns]

Definition at line 85 of file JMatrixNZ.hh.

     {
       using namespace std;
       using namespace JPP;
 
 
       const int N = distance(__begin, __end);
 
       this ->resize(N);
       buffer.resize(N);
 
       const double ta = alpha * PI / 180.0;
       const double ct = cos(ta);
       const double st = sin(ta);
 
 
       // angular resolution
 
       for (T i = __begin; i != __end; ++i) {
         
         const double dx = i->getX() - pos.getX(); 
         const double dy = i->getY() - pos.getY(); 
         const double dz = i->getZ() - pos.getZ(); 
 
         const double R  = sqrt(dx*dx + dy*dy);
         
         double x  = ta * getKappaC() * getInverseSpeedOfLight();
         double y  = ta * getKappaC() * getInverseSpeedOfLight();
         double v  = ta * getInverseSpeedOfLight();
         double w  = ta * getInverseSpeedOfLight();
 
         if (R != 0.0) {
           x *= dx / R;
           y *= dy / R;
         }
 
         x *=  (dz*ct - dx*st);
         y *=  (dz*ct - dy*st);
         v *= -(dx*ct + dz*st);
         w *= -(dy*ct + dz*st);
         
         buffer[distance(__begin,i)] = variance(x, y, v, w);
       }
 
       for (int i = 0; i != N; ++i) {
 
         for (int j = 0; j != i; ++j) {
           (*this)(i, j) = getDot(buffer[i], buffer[j]);
           (*this)(j, i) = (*this)(i, j);
         }
 
         (*this)(i, i) = getDot(buffer[i], buffer[i]) +  sigma * sigma;
       }
     }

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

inlinestaticprivate

Get dot product.

Parameters

first	first variance
second	second variance

Returns: dot product

Definition at line 196 of file JMatrixNZ.hh.

     {
       return (first.x * second.x  +
               first.y * second.y  +
               first.v * second.v  +
               first.w * second.w);
               
     }

void JMATH::JMatrixNS::invert ( )

inlineinherited

Invert matrix according LDU decomposition.

Definition at line 80 of file JMatrixNS.hh.

     {
       using namespace std;
 
       size_t  i, row, col, root;
 
       double* p0;
       double* p1;
       double* p2;
       double* p3;
 
       const double* c0;
       const double* c1;
       const double* c2;
       const double* c3;
 
       const double* a0;
       const double* a1;
       const double* a2;
       const double* a3;
 
       double v0;
       double v1;
       double v2;
       double v3;
 
       double w0;
       double w1;
       double w2;
       double w3;          
 
       JPermutationMatrix P;
 
 
       // LDU decomposition
 
       for (root = 0; root != this->size(); ++root) {
 
         double value = (*this)(root, root);
         size_t pivot = root;
 
         for (row = root; ++row != this->size(); ) {
           if (fabs((*this)(row, root)) > fabs(value)) {
             value = (*this)(row, root);
             pivot = row;
           }
         }
 
         if (value == 0.0) {
           THROW(JDivisionByZero, "LDU decomposition zero pivot at " << root);
         }
 
         if (pivot != root) {
 
           this->rswap(root, pivot);
 
           P.push_back(make_pair(root, pivot));
         }
 
         for (row = root + 1; row + 4 <= this->size(); row += 4) {            // process rows by four
 
           p0 = (*this)[row + 0]  +  root;
           p1 = (*this)[row + 1]  +  root;
           p2 = (*this)[row + 2]  +  root;
           p3 = (*this)[row + 3]  +  root;
 
           v0 = (*p0++ /= value);
           v1 = (*p1++ /= value);
           v2 = (*p2++ /= value);
           v3 = (*p3++ /= value);
 
           c0 = (*this)[root] + root + 1;
 
           for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, p1 += 4, p2 += 4, p3 += 4, c0 += 4) {
             p0[0] -= v0 * c0[0];  p0[1] -= v0 * c0[1];  p0[2] -= v0 * c0[2];  p0[3] -= v0 * c0[3];
             p1[0] -= v1 * c0[0];  p1[1] -= v1 * c0[1];  p1[2] -= v1 * c0[2];  p1[3] -= v1 * c0[3];
             p2[0] -= v2 * c0[0];  p2[1] -= v2 * c0[1];  p2[2] -= v2 * c0[2];  p2[3] -= v2 * c0[3];
             p3[0] -= v3 * c0[0];  p3[1] -= v3 * c0[1];  p3[2] -= v3 * c0[2];  p3[3] -= v3 * c0[3];
           }
 
           for ( ; i != this->size(); ++i, ++p0, ++p1, ++p2, ++p3, ++c0) {
             *p0 -= v0 * (*c0);
             *p1 -= v1 * (*c0);
             *p2 -= v2 * (*c0);
             *p3 -= v3 * (*c0);
           }
         }
 
         for ( ; row != this->size(); ++row) {                                // remaining rows
 
           p0 = (*this)[row + 0]  +  root;
 
           v0 = (*p0++ /= value);
 
           c0 = (*this)[root] + root + 1;
 
           for (i = root + 1; i + 4 <= this->size(); i += 4, p0 += 4, c0 += 4) {
             *p0 -= v0 * c0[0];  p0[1] -= v0 * c0[1];  p0[2] -= v0 * c0[2];  p0[3] -= v0 * c0[3];
           }
 
           for ( ; i != this->size(); ++i, ++p0, ++c0) {
             *p0 -= v0 * (*c0);
           }
         }
       }
 
 
       // invert D
 
       for (size_t i = 0; i != this->size(); ++i) {
 
         double& pivot = (*this)(i,i);
 
         if (pivot == 0.0) {
           THROW(JDivisionByZero, "Invert D zero pivot at " << i);
         }
 
         pivot = 1.0 / pivot;
       }
 
 
       // invert U
 
       JMatrixND_t& A = getInstance();
 
       A.set(*this);
       A.transpose();
 
       for (row = 0; row + 4 <= this->size(); row += 4) {                     // process rows by four
 
         p0 = (*this)[row + 0]  +  row + 0;
         p1 = (*this)[row + 1]  +  row + 1;
         p2 = (*this)[row + 2]  +  row + 2;
         p3 = (*this)[row + 3]  +  row + 3;
         
         v0 = *(p0++);
         v1 = *(p1++);
         v2 = *(p2++);
         v3 = *(p3++);
 
         for (col = row + 1; col + 4 <= this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
 
           w0 = (*p0) * -v0;
           w1 = (*p1) * -v1;
           w2 = (*p2) * -v2;
           w3 = (*p3) * -v3;
 
           c0 = (*this)[row + 0]  +  row + 0 + 1;
           c1 = (*this)[row + 1]  +  row + 1 + 1;
           c2 = (*this)[row + 2]  +  row + 2 + 1;
           c3 = (*this)[row + 3]  +  row + 3 + 1;
 
           a0 = A[col + 0]  +  row + 0 + 1;
           a1 = A[col + 1]  +  row + 1 + 1;
           a2 = A[col + 2]  +  row + 2 + 1;
           a3 = A[col + 3]  +  row + 3 + 1;
 
           for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
             w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3]; 
             w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3]; 
             w2 -= c2[0]*a2[0] + c2[1]*a2[1] + c2[2]*a2[2] + c2[3]*a2[3]; 
             w3 -= c3[0]*a3[0] + c3[1]*a3[1] + c3[2]*a3[2] + c3[3]*a3[3]; 
           }
 
           for ( ; c0 != p0; ++c0, ++c1, ++c2, ++c3, ++a0, ++a1, ++a2, ++a3) {
             w0 -= (*c0) * (*a0);
             w1 -= (*c1) * (*a1);
             w2 -= (*c2) * (*a2);
             w3 -= (*c3) * (*a3);
           }
 
           *p0 = w0 * (*this)(col + 0, col + 0);
           *p1 = w1 * (*this)(col + 1, col + 1);
           *p2 = w2 * (*this)(col + 2, col + 2);
           *p3 = w3 * (*this)(col + 3, col + 3);
         }
 
         // last 3 columns
 
         w0 = (*p0) * -v0;
         w1 = (*p1) * -v1;
         w2 = (*p2) * -v2;
 
         c0 = (*this)[row + 0]  +  row + 0 + 1;
         c1 = (*this)[row + 1]  +  row + 1 + 1;
         c2 = (*this)[row + 2]  +  row + 2 + 1;
 
         a0 = A[col + 0]  +  row + 0 + 1;
         a1 = A[col + 1]  +  row + 1 + 1;
         a2 = A[col + 2]  +  row + 2 + 1;
 
         for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
           w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3]; 
           w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3]; 
           w2 -= c2[0]*a2[0] + c2[1]*a2[1] + c2[2]*a2[2] + c2[3]*a2[3]; 
         }
 
         for ( ; c0 != p0; ++c0, ++c1, ++c2, ++a0, ++a1, ++a2) {
           w0 -= (*c0) * (*a0);
           w1 -= (*c1) * (*a1);
           w2 -= (*c2) * (*a2);
         }
           
         *p0 = w0 * (*this)(col + 0, col + 0);
         *p1 = w1 * (*this)(col + 1, col + 1);
         *p2 = w2 * (*this)(col + 2, col + 2);
 
         ++col; ++p0; ++p1;
 
         // last 2 columns
 
         w0 = (*p0) * -v0;
         w1 = (*p1) * -v1;
 
         c0 = (*this)[row + 0]  +  row + 0 + 1;
         c1 = (*this)[row + 1]  +  row + 1 + 1;
 
         a0 = A[col + 0]  +  row + 0 + 1;
         a1 = A[col + 1]  +  row + 1 + 1;
 
         for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
           w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3]; 
           w1 -= c1[0]*a1[0] + c1[1]*a1[1] + c1[2]*a1[2] + c1[3]*a1[3]; 
         }
 
         for ( ; c0 != p0; ++c0, ++c1, ++a0, ++a1) {
           w0 -= (*c0) * (*a0);
           w1 -= (*c1) * (*a1);
         }
           
         *p0 = w0 * (*this)(col + 0, col + 0);
         *p1 = w1 * (*this)(col + 1, col + 1);
 
         ++col; ++p0;
 
         // last column
 
         w0 = (*p0) * -v0;
 
         c0 = (*this)[row + 0]  +  row + 0 + 1;
 
         a0 = A[col + 0]  +  row + 0 + 1;
 
         for ( ; c0 + 4 < p0; c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4) {
           w0 -= c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3]; 
         }
 
         for ( ; c0 != p0; ++c0, ++a0) {
           w0 -= (*c0) * (*a0);
         }
           
         *p0 = w0 * (*this)(col + 0, col + 0);
       }
 
       for ( ; row != this->size(); ++row) {                                  // remaining rows
 
         p0 = (*this)[row + 0]  +  row + 0;
 
         v0 = *(p0++);
 
         for (col = row + 1; col != this->size(); ++col, ++p0) {
 
           w0 = *p0;
 
           w0 *= -v0;
 
           c0 = (*this)[row + 0]  +  row + 1;
 
           a0 = A[col + 0]  +  row + 0 + 1;
 
           for ( ; c0 != p0; ++c0, ++a0) {
             w0 -= (*c0) * (*a0);
           }
           
           *p0 = w0 * (*this)(col + 0, col + 0);
         }
       }
 
 
       // solve A^-1 x L = U^-1
 
       double* buffer = A.data();
 
       for (col = this->size() - 1; col != 0; --col) {
 
         for (i = col; i < this->size(); ++i) {
 
           buffer[i] = (*this)(i, col - 1);
 
           (*this)(i, col - 1) = 0.0;
         }
 
         for (row = 0; row + 4 < this->size(); row += 4) {
 
           p0 = (*this)[row + 0] + col - 1;
           p1 = (*this)[row + 1] + col - 1;
           p2 = (*this)[row + 2] + col - 1;
           p3 = (*this)[row + 3] + col - 1;
 
           c0 = p0 + 1;
           c1 = p1 + 1;
           c2 = p2 + 1;
           c3 = p3 + 1;
 
           a0 = buffer + col;
 
           w0 = 0.0;
           w1 = 0.0;
           w2 = 0.0;
           w3 = 0.0;
 
           for (i = col; i + 4 <= this->size(); i += 4, c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4) {
             w0 += c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
             w1 += c1[0]*a0[0] + c1[1]*a0[1] + c1[2]*a0[2] + c1[3]*a0[3];
             w2 += c2[0]*a0[0] + c2[1]*a0[1] + c2[2]*a0[2] + c2[3]*a0[3];
             w3 += c3[0]*a0[0] + c3[1]*a0[1] + c3[2]*a0[2] + c3[3]*a0[3];
           }
 
           for ( ; i != this->size(); ++i, ++c0, ++c1, ++c2, ++c3, ++a0) {
             w0 += (*c0) * (*a0);
             w1 += (*c1) * (*a0);
             w2 += (*c2) * (*a0);
             w3 += (*c3) * (*a0);
           }
 
           *p0 -= w0;
           *p1 -= w1;
           *p2 -= w2;
           *p3 -= w3;
         }
 
         for ( ; row < this->size(); ++row) {
 
           p0 = (*this)[row] + col - 1;
 
           c0 = p0 + 1;
 
           a0 = buffer + col;
 
           w0 = 0.0;
 
           for (i = col; i + 4 <= this->size(); i += 4, c0 += 4, c1 += 4, c2 += 4, c3 += 4, a0 += 4) {
             w0 += c0[0]*a0[0] + c0[1]*a0[1] + c0[2]*a0[2] + c0[3]*a0[3];
           }
 
           for ( ; i != this->size(); ++i, ++c0, ++a0) {
             w0 += (*c0) * (*a0);
           }
 
           *p0 -= w0;
         }
       }
 
 
       for (JPermutationMatrix::reverse_iterator p = P.rbegin(); p != P.rend(); ++p) {
         cswap(p->first, p->second);
       }
     }

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

k	index of diagonal element
value	value to add

Definition at line 457 of file JMatrixNS.hh.

     {
       if (k >= this->size()) {
         THROW(JIndexOutOfRange, "JMatrixNS::update(): index out of range " << k);
       }
 
       size_t  i, row;
       double* col;
 
       const double del  =  (*this)(k,k);
       const double din  =  1.0 / (1.0 + value * del);
 
       double* root = (*this)[k];
 
       for (row = 0; row != this->size(); ++row) {
 
         if (row != k) {
 
           const double val = (*this)(row, k);
 
           for (i = 0, col = (*this)[row]; i != this->size(); ++i, ++col) {
             if (i != k) {
               (*col) -= val * root[i] * value * din;
             }
           }
 
           (*this)(row, k) *= din;
         }
       }
 
       for (i = 0, col = root; i != this->size(); ++i, ++col) {
         if (i != k) {
           (*col) *= din;
         }
       }
 
       root[k] = del * din;
     }

JMatrixND_t& JMATH::JMatrixND::getInstance ( )

inlineprotectedinherited

Get work space.

Returns: work space

Definition at line 336 of file JMatrixND.hh.

     {
       return ws;
     }

void JMATH::JMatrixND::resize ( const size_t size )

inlineinherited

Resize matrix.

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

Parameters

size dimension

Definition at line 421 of file JMatrixND.hh.

     {
       static_cast<JMatrixND_t&>(*this).resize(size);
 
       getInstance().resize(size);
     }

JMatrixND& JMATH::JMatrixND::reset ( )

inlineinherited

Set matrix to the null matrix.

Returns: this matrix

Definition at line 434 of file JMatrixND.hh.

     {
       JMatrixND_t& A = getInstance();
 
       A.resize(this->size());
       
       double* p = A.data();
 
       for (size_t i = this->size(); i != 0; --i, ++p) {
         *p = 0.0;
       }
 
       p = this->data();
       
       for (size_t i = this->size(); i != 0; --i, p += this->size()) {
         memcpy(p, A.data(), this->size() * sizeof(double));
       }
 
       return *this;
     }

JMatrixND& JMATH::JMatrixND::setIdentity ( )

inlineinherited

Set to identity matrix.

Returns: this matrix

Definition at line 461 of file JMatrixND.hh.

     {
       reset();
 
       for (size_t i = 0; i != this->size(); ++i) {
         (*this)(i,i) = 1.0;
       }
 
       return *this;
     }

JMatrixND& JMATH::JMatrixND::negate ( )

inlineinherited

Negate matrix.

Returns: this matrix

Definition at line 478 of file JMatrixND.hh.

     {
       double* p = this->data();
 
       for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
         *p = -*p;
       }
 
       return *this;
     }

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

inlineinherited

Matrix addition.

Parameters

A matrix

Returns: this matrix

Definition at line 496 of file JMatrixND.hh.

     {
       if (this->size() == A.size()) {
 
         double*       p = this->data();
         const double* q = A.data();
 
         for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
           *p += *q;
         }
 
       } else {
         THROW(JException, "JMatrixND::add() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
       }
     }

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

inlineinherited

Matrix subtraction.

Parameters

A matrix

Returns: this matrix

Definition at line 519 of file JMatrixND.hh.

     {
       if (this->size() == A.size()) {
 
         double*       p = this->data();
         const double* q = A.data();
 
         for (size_t i = this->size()*this->size(); i != 0; --i, ++p, ++q) {
           *p -= *q;
         }
 
       } else {
         THROW(JException, "JMatrixND::sub() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
       }
 
       return *this;
     }

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

inlineinherited

Scale matrix.

Parameters

factor factor

Returns: this matrix

Definition at line 544 of file JMatrixND.hh.

     {
       double* p = this->data();
 
       for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
         *p *= factor;
       }
 
       return *this;
     }

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

inlineinherited

Matrix multiplication.

Parameters

A	matrix
B	matrix

Returns: this matrix

Definition at line 581 of file JMatrixND.hh.

     {
       if (A.size() == B.size()) {
 
         this->resize(A.size());
 
         if (!this->empty()) {
 
           JMatrixND_t& C = getInstance();
           
           C.set(B);
           C.transpose();
 
           size_t i, row;
 
           for (row = 0; row + 4 <= this->size(); row += 4) {              // process rows by four
 
             double* p0 = (*this)[row + 0];
             double* p1 = (*this)[row + 1];
             double* p2 = (*this)[row + 2];
             double* p3 = (*this)[row + 3];
 
             for (size_t col = 0; col != this->size(); ++col, ++p0, ++p1, ++p2, ++p3) {
 
               double w0 = 0.0;
               double w1 = 0.0;
               double w2 = 0.0;
               double w3 = 0.0;
 
               const double* a0 =  A[row + 0];
               const double* a1 =  A[row + 1];
               const double* a2 =  A[row + 2];
               const double* a3 =  A[row + 3];
 
               const double* c0 =  C[col];
 
               for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, a1 += 4, a2 += 4, a3 += 4, c0 += 4) {
                 w0  +=  a0[0] * c0[0]  +  a0[1] * c0[1]  +  a0[2] * c0[2]  +  a0[3] * c0[3];
                 w1  +=  a1[0] * c0[0]  +  a1[1] * c0[1]  +  a1[2] * c0[2]  +  a1[3] * c0[3];
                 w2  +=  a2[0] * c0[0]  +  a2[1] * c0[1]  +  a2[2] * c0[2]  +  a2[3] * c0[3];
                 w3  +=  a3[0] * c0[0]  +  a3[1] * c0[1]  +  a3[2] * c0[2]  +  a3[3] * c0[3];
               }
 
               for ( ; i != this->size(); ++i, ++a0, ++a1, ++a2, ++a3, ++c0) {
                 w0  +=  (*a0) * (*c0);
                 w1  +=  (*a1) * (*c0);
                 w2  +=  (*a2) * (*c0);
                 w3  +=  (*a3) * (*c0);
               }
 
               *p0 = w0;
               *p1 = w1;
               *p2 = w2;
               *p3 = w3;
             }
           }
 
           for ( ; row != this->size(); ++row) {                  // remaining rows
 
             double* p0 = (*this)[row + 0];
 
             for (size_t col = 0; col != this->size(); ++col, ++p0) {
 
               double w0 = 0.0;
 
               const double* a0 =  A[row + 0];
               const double* c0 =  C[col];
 
               for (i = 0; i + 4 <= this->size(); i += 4, a0 += 4, c0 += 4) {
                 w0  +=  a0[0] * c0[0]  +  a0[1] * c0[1]  +  a0[2] * c0[2]  +  a0[3] * c0[3];
               }
 
               for ( ; i != this->size(); ++i, ++a0, ++c0) {
                 w0  +=  (*a0) * (*c0);
               }
 
               *p0 = w0;
             }
           }
         }
 
       } else {
         THROW(JException, "JMatrixND::mul() inconsistent matrix dimensions " << A.size() << ' ' << B.size());
       }
 
       return *this;
     }

JMatrixND & JMATH::JMath< JMatrixND , JNullType >::mul ( const JNullType & object )

inlineinherited

Multiply with object.

Parameters

object object

Returns: result object

Definition at line 354 of file JMath.hh.

     {
       return static_cast<JFirst_t&>(*this) = JFirst_t().mul(static_cast<const JFirst_t&>(*this), object);
     }

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

inlineinherited

Scale matrix.

Parameters

factor factor

Returns: this matrix

Definition at line 562 of file JMatrixND.hh.

     {
       double* p = this->data();
 
       for (size_t i = this->size()*this->size(); i != 0; --i, ++p) {
         *p /= factor;
       }
 
       return *this;
     }

bool JMATH::JMatrixND::equals	(	const JMatrixND &	A,
		const double	eps = `std::numeric_limits<double>::min()`
	)		const

inlineinherited

Equality.

Parameters

A	matrix
eps	numerical precision

Returns: true if matrices identical; else false

Definition at line 678 of file JMatrixND.hh.

     {
       if (this->size() == A.size()) {
 
         for (size_t row = 0; row != this->size(); ++row) {
 
           const double* p = (*this)[row];
           const double* q = A      [row];
 
           for (size_t i = this->size(); i != 0; --i, ++p, ++q) {
             if (fabs(*p - *q) > eps) {
               return false;
             }
           }
         }
 
         return true;
 
       } else {
         THROW(JException, "JMatrixND::equals() inconsistent matrix dimensions " << this->size() << ' ' << A.size());
       }
     }

bool JMATH::JMatrixND::isIdentity ( const double eps = std::numeric_limits<double>::min() ) const

inlineinherited

Test identity.

Parameters

eps	numerical precision

Returns: true if identity matrix; else false

Definition at line 709 of file JMatrixND.hh.

     {
       for (size_t i = 0; i != this->size(); ++i) {
 
         if (fabs(1.0 - (*this)(i,i)) > eps) {
           return false;
         };
 
         for (size_t j = 0; j != i; ++j) {
           if (fabs((*this)(i,j)) > eps || fabs((*this)(j,i)) > eps) {
             return false;
           };
         }
       }
 
       return true;
     }

bool JMATH::JMatrixND::isSymmetric ( const double eps = std::numeric_limits<double>::min() ) const

inlineinherited

Test symmetry.

Parameters

eps	numerical precision

Returns: true if symmetric matrix; else false

Definition at line 734 of file JMatrixND.hh.

     {
       for (size_t i = 0; i != this->size(); ++i) {
         for (size_t j = 0; j != i; ++j) {
           if (fabs((*this)(i,j) - (*this)(j,i)) > eps) {
             return false;
           };
         }
       }
 
       return true;
     }

double JMATH::JMatrixND::getDot ( const JVectorND & v ) const

inlineinherited

Get dot product.

The dot product corresponds to

$v^{T} \times A \times v$

.

Parameters

v vector

Returns: dot product

Definition at line 756 of file JMatrixND.hh.

     {
       double dot = 0.0;
 
       for (size_t i = 0; i != v.size(); ++i) {
 
         const double* p = (*this)[i];
 
         double w = 0.0;
 
         for (JVectorND::const_iterator y = v.begin(); y != v.end(); ++y, ++p) {
           w += (*p) * (*y);
         }
 
         dot += v[i] * w;
       }
 
       return dot;
     }

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

inlineprotectedinherited

Definition at line 856 of file JMatrixND.hh.

     {
       JMatrixND_t& A = getInstance();
 
       A.resize(this->size());
       
       memcpy(A.data(),    (*this)[r1], this->size() * sizeof(double));
       memcpy((*this)[r1], (*this)[r2], this->size() * sizeof(double));
       memcpy((*this)[r2], A.data(),    this->size() * sizeof(double));
     }

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

inlineprotectedinherited

Swap columns.

Parameters

c1	column
c2	column

Definition at line 874 of file JMatrixND.hh.

     {
       using std::swap;
 
       double* p1 = this->data() + c1;
       double* p2 = this->data() + c2;
 
       for (size_t i = this->size(); i != 0; --i, p1 += this->size(), p2 += this->size()) {
         swap(*p1, *p2);
       }
     }

void JMATH::JMatrixND_t::clear ( )

inlineinherited

Clear memory.

Definition at line 83 of file JMatrixND.hh.

     {
       if (__p != NULL) {
         delete [] __p;
       }
 
       __p = NULL;
       __n = 0;
       __m = 0;
     }

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

inlineinherited

Set matrix.

Parameters

A matrix

Definition at line 130 of file JMatrixND.hh.

     {
       this->resize(A.size());
 
       memcpy(this->data(), A.data(), A.size() * A.size() * sizeof(double));
     }

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

inlineinherited

Swap matrices.

Parameters

A matrix

Definition at line 143 of file JMatrixND.hh.

     {
       using std::swap;
 
       swap(this->__p, A.__p);
       swap(this->__n, A.__n);
       swap(this->__m, A.__m);
     }

JMatrixND_t& JMATH::JMatrixND_t::transpose ( )

inlineinherited

Transpose.

Returns: this matrix

Definition at line 158 of file JMatrixND.hh.

     {
       using std::swap;
 
       for (size_t row = 0; row != this->size(); ++row) {
         for (size_t col = 0; col != row; ++col) {
           swap((*this)(row,col), (*this)(col,row));
         }
       }
 
       return *this;
     }

size_t JMATH::JMatrixND_t::size ( ) const

inlineinherited

Get dimension of matrix.

Returns: dimension

Definition at line 177 of file JMatrixND.hh.

     {
       return __n;
     }

size_t JMATH::JMatrixND_t::capacity ( ) const

inlineinherited

Get capacity of dimension.

Returns: capacity

Definition at line 188 of file JMatrixND.hh.

     {
       return __m;
     }

bool JMATH::JMatrixND_t::empty ( ) const

inlineinherited

Check emptyness.

Returns: true if empty; else false

Definition at line 199 of file JMatrixND.hh.

     {
       return __n == 0;
     }

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

inlineinherited

Get pointer to data.

Returns: pointer to data.

Definition at line 210 of file JMatrixND.hh.

     {
       return __p;
     }

double* JMATH::JMatrixND_t::data ( )

inlineinherited

Get pointer to data.

Returns: pointer to data.

Definition at line 221 of file JMatrixND.hh.

     {
       return __p;
     }

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

inlineinherited

Get row data.

Parameters

row	row number

Returns: pointer to row data

Definition at line 233 of file JMatrixND.hh.

     {
       return __p + row*__n;
     }

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

inlineinherited

Get row data.

Parameters

row	row number

Returns: pointer to row data

Definition at line 245 of file JMatrixND.hh.

     {
       return __p + row*__n;
     }

double JMATH::JMatrixND_t::operator()	(	const size_t	row,
		const size_t	col
	)		const

inlineinherited

Get matrix element.

Parameters

row	row number
col	column number

Returns: matrix element at (row,col)

Definition at line 258 of file JMatrixND.hh.

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

double& JMATH::JMatrixND_t::operator()	(	const size_t	row,
		const size_t	col
	)

inlineinherited

Get matrix element.

Parameters

row	row number
col	column number

Returns: matrix element at (row,col)

Definition at line 271 of file JMatrixND.hh.

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

double JMATH::JMatrixND_t::at	(	size_t	row,
		size_t	col
	)		const

inlineinherited

Get matrix element.

Parameters

row	row number
col	column number

Returns: matrix element at (row,col)

Definition at line 284 of file JMatrixND.hh.

     {
       if (row >= this->size() ||
           col >= this->size()) {
         THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
       }
 
       return (*this)(row, col);
     }

double& JMATH::JMatrixND_t::at	(	size_t	row,
		size_t	col
	)

inlineinherited

Get matrix element.

Parameters

row	row number
col	column number

Returns: matrix element at (row,col)

Definition at line 302 of file JMatrixND.hh.

     {
       if (row >= this->size() ||
           col >= this->size()) {
         THROW(JIndexOutOfRange, "JMatrixND::at(" << row << "," << col << "): index out of range.");
       }
 
       return (*this)(row, col);
     }

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:

JMatrixNZ.hh

Classes

Public Types

Public Member Functions

Protected Member Functions

Protected Attributes

Static Private Member Functions

Private Attributes

Detailed Description

Member Typedef Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation