N x N symmetric matrix. More...

#include <JMatrixNS.hh>

Inheritance diagram for JMATH::JMatrixNS:

Public Member Functions
	JMatrixNS ()
	Default constructor. More...

	JMatrixNS (const size_t size)
	Constructor (identity matrix). More...

	JMatrixNS (const JMatrixND &A)
	Contructor. More...

JMatrixNS &	operator= (const JMatrixNS &A)
	Assignment operator. More...

void	invert ()
	Invert matrix according LDU decomposition. More...

template<class JVectorND_t >
void	solve (JVectorND_t &u)
	Get solution of equation `A x = b`. More...

void	update (const size_t k, const double value)
	Update inverted matrix at given diagonal element. More...

void	resize (const size_t size)
	Resize matrix. More...

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

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

Detailed Description

N x N symmetric matrix.

Definition at line 28 of file JMatrixNS.hh.

Constructor & Destructor Documentation

◆ JMatrixNS() [1/3]

JMATH::JMatrixNS::JMatrixNS ( )

inline

Default constructor.

Definition at line 34 of file JMatrixNS.hh.

                 :
       JMatrixND()
     {}

◆ JMatrixNS() [2/3]

JMATH::JMatrixNS::JMatrixNS ( const size_t size )

inline

Constructor (identity matrix).

Parameters

size dimension

Definition at line 44 of file JMatrixNS.hh.

                                  :
       JMatrixND(size)
     {}

◆ JMatrixNS() [3/3]

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

inline

Contructor.

Parameters

A matrix

Definition at line 54 of file JMatrixNS.hh.

                                   :
       JMatrixND(A)
     {}

Member Function Documentation

◆ operator=()

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

inline

Assignment operator.

Parameters

A matrix

Definition at line 64 of file JMatrixNS.hh.

     {
       this->set(A);
  
       return *this;
     }

◆ invert()

void JMATH::JMatrixNS::invert ( )

inline

Invert matrix according LDU decomposition.

Definition at line 75 of file JMatrixNS.hh.

     {
       using namespace std;
  
       size_t  i, row, col, root;
  
       double* p0;
       double* p1;
       double* p2;
       double* p3;
  
       const double* c0;
  
       double v0;
       double v1;
       double v2;
       double v3;
  
       double w0;
  
       P.resize(this->size());
  
       // 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[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);
           }
         }
       }
  
  
       JMatrixND_t& A = getInstance();
  
       A.reset();
  
       for (i = 0; i != this->size(); ++i) {
         A(i,i) = 1.0;
       }
  
  
       // forward substitution
  
       col = 0;
       
       for ( ; col + 4 <= this->size(); col += 4) {
           
         for (row = 0; row != this->size(); ++row) {
  
           p0 = A[col + 0];
           p1 = A[col + 1];
           p2 = A[col + 2];
           p3 = A[col + 3];
  
           i  = P[row];
  
           v0 = p0[i];
           v1 = p1[i];
           v2 = p2[i];
           v3 = p3[i];
  
           p0[i] = p0[row];
           p1[i] = p1[row];
           p2[i] = p2[row];
           p3[i] = p3[row];
  
           for (i = 0, c0 = (*this)[row]; i != row; ++c0, ++p0, ++p1, ++p2, ++p3, ++i) {
             v0 -= (*c0) * (*p0);
             v1 -= (*c0) * (*p1);
             v2 -= (*c0) * (*p2);
             v3 -= (*c0) * (*p3);
           }
  
           A[col + 0][row] = v0;
           A[col + 1][row] = v1;
           A[col + 2][row] = v2;
           A[col + 3][row] = v3;
         }
       }
  
       for ( ; col != this->size(); ++col) {
  
         for (row = 0; row != this->size(); ++row) {
  
           p0 = A[col];
           
           i  = P[row];
  
           v0 = p0[i];
  
           p0[i] = p0[row];
  
           for (i = 0, c0 = (*this)[row] + i; i != row; ++c0, ++p0, ++i) {
             v0 -= (*c0) * (*p0);
           }
  
           A[col][row] = v0;
         }
       }
  
       // backward substitution
  
       col = 0;
       
       for ( ; col + 4 <= this->size(); col += 4) {
  
         for (int row = this->size() - 1; row >= 0; --row) {
  
           p0 = A[col + 0] + row;
           p1 = A[col + 1] + row;
           p2 = A[col + 2] + row;
           p3 = A[col + 3] + row;
  
           v0 = *p0;
           v1 = *p1;
           v2 = *p2;
           v3 = *p3;
  
           ++p0;
           ++p1;
           ++p2;
           ++p3;
  
           for (i = row + 1, c0 = (*this)[row] + i; i != this->size(); ++c0, ++p0, ++p1, ++p2, ++p3, ++i) {
             v0 -= (*c0) * (*p0);
             v1 -= (*c0) * (*p1);
             v2 -= (*c0) * (*p2);
             v3 -= (*c0) * (*p3);
           }
  
           w0 = 1.0 / (*this)[row][row];
  
           A[col + 0][row] = v0 * w0;
           A[col + 1][row] = v1 * w0;
           A[col + 2][row] = v2 * w0;
           A[col + 3][row] = v3 * w0;
         }
       }
      
       for ( ; col != this->size(); ++col) {
  
         for (int row = this->size() - 1; row >= 0; --row) {
  
           p0 = A[col] + row;
  
           v0 = *p0;
  
           ++p0;
  
           for (i = row + 1, c0 = (*this)[row] + i; i != this->size(); ++c0, ++p0, ++i) {
             v0 -= (*c0) * (*p0);
           }
  
           w0 = 1.0 / (*this)[row][row];
  
           A[col][row] = v0 * w0;
         }
       }
  
       A.transpose();
  
       this->swap(A);
     }

◆ solve()

template<class JVectorND_t >

void JMATH::JMatrixNS::solve ( JVectorND_t & u )

inline

Get solution of equation A x = b.

Parameters

u	column vector; b on input, x on output(I/O)

Definition at line 308 of file JMatrixNS.hh.

     {
       using namespace std;
  
       size_t  i, row, root;
  
       double* p0;
       double* p1;
       double* p2;
       double* p3;
  
       const double* c0;
  
       double v0;
       double v1;
       double v2;
       double v3;
  
       P.resize(this->size());
  
       // 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[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);
           }
         }
       }
  
       // forward substitution
  
       for (row = 0; row != this->size(); ++row) {
  
         i      = P[row];
         v0     = u[i];
         u[i]   = u[row];
  
         for (i = 0, c0 = (*this)[row] + i; i != row; ++c0, ++i) {
           v0 -= *c0 * u[i];
         }
  
         u[row] = v0;
       }
  
       // backward substitution
  
       for (int row = this->size() - 1; row >= 0; --row) {
  
         v0     = u[row];
  
         for (i = row + 1, c0 = (*this)[row] + i; i < this->size(); ++c0, ++i) {
           v0 -= *c0 * u[i];
         }
  
         u[row] = v0 / (*this)[row][row];
       }
     }

◆ update()

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

inline

Update inverted matrix at given diagonal element.

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

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

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

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

Parameters

k	index of diagonal element
value	value to add

Definition at line 446 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;
     }

◆ getInstance()

JMatrixND_t& JMATH::JMatrixND::getInstance ( )

inlineprotectedinherited

Get work space.

Returns: work space

Definition at line 361 of file JMatrixND.hh.

     {
       return ws;
     }

◆ resize()

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

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

◆ reset()

JMatrixND& JMATH::JMatrixND::reset ( )

inlineinherited

Set matrix to the null matrix.

Returns: this matrix

Definition at line 459 of file JMatrixND.hh.

     {
       JMatrixND_t::reset();
  
       JMatrixND_t& A = getInstance();
  
       A.resize(this->size());
  
       return *this;
     }

◆ setIdentity()

JMatrixND& JMATH::JMatrixND::setIdentity ( )

inlineinherited

Set to identity matrix.

Returns: this matrix

Definition at line 476 of file JMatrixND.hh.

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

◆ negate()

JMatrixND& JMATH::JMatrixND::negate ( )

inlineinherited

Negate matrix.

Returns: this matrix

Definition at line 493 of file JMatrixND.hh.

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

◆ add()

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

inlineinherited

Matrix addition.

Parameters

A matrix

Returns: this matrix

Definition at line 511 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());
       }
  
       return *this;
     }

◆ sub()

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

inlineinherited

Matrix subtraction.

Parameters

A matrix

Returns: this matrix

Definition at line 536 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;
     }

◆ mul() [1/3]

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

inlineinherited

Scale matrix.

Parameters

factor factor

Returns: this matrix

Definition at line 561 of file JMatrixND.hh.

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

◆ mul() [2/3]

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

inlineinherited

Matrix multiplication.

Parameters

A	matrix
B	matrix

Returns: this matrix

Definition at line 598 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;
     }

◆ mul() [3/3]

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

◆ div()

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

inlineinherited

Scale matrix.

Parameters

factor factor

Returns: this matrix

Definition at line 579 of file JMatrixND.hh.

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

◆ equals()

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 695 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());
       }
     }

◆ isIdentity()

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 726 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;
     }

◆ isSymmetric()

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 751 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;
     }

◆ getDot()

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 773 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;
     }

◆ rswap()

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

inlineprotectedinherited

Definition at line 873 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));
     }

◆ cswap()

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

inlineprotectedinherited

Swap columns.

Parameters

c1	column
c2	column

Definition at line 891 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);
       }
     }

◆ clear()

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;
     }

◆ set()

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

inlineinherited

Set matrix.

Parameters

A matrix

Definition at line 155 of file JMatrixND.hh.

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

◆ swap()

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

inlineinherited

Swap matrices.

Parameters

A matrix

Definition at line 168 of file JMatrixND.hh.

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

◆ transpose()

JMatrixND_t& JMATH::JMatrixND_t::transpose ( )

inlineinherited

Transpose.

Returns: this matrix

Definition at line 183 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()

size_t JMATH::JMatrixND_t::size ( ) const

inlineinherited

Get dimension of matrix.

Returns: dimension

Definition at line 202 of file JMatrixND.hh.

     {
       return __n;
     }

◆ capacity()

size_t JMATH::JMatrixND_t::capacity ( ) const

inlineinherited

Get capacity of dimension.

Returns: capacity

Definition at line 213 of file JMatrixND.hh.

     {
       return __m;
     }

◆ empty()

bool JMATH::JMatrixND_t::empty ( ) const

inlineinherited

Check emptyness.

Returns: true if empty; else false

Definition at line 224 of file JMatrixND.hh.

     {
       return __n == 0;
     }

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

     {
       return __p;
     }

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

     {
       return __p;
     }

◆ operator[]() [1/2]

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

     {
       return __p + row*__n;
     }

◆ operator[]() [2/2]

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

inlineinherited

Get row data.

Parameters

row	row number

Returns: pointer to row data

Definition at line 270 of file JMatrixND.hh.

     {
       return __p + row*__n;
     }

◆ operator()() [1/2]

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

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

◆ operator()() [2/2]

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

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

◆ at() [1/2]

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 309 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);
     }

◆ at() [2/2]

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

◆ P

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

private

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 struct was generated from the following file:

JMatrixNS.hh

Public Member Functions

Protected Member Functions

Protected Attributes

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ JMatrixNS() [1/3]

◆ JMatrixNS() [2/3]

◆ JMatrixNS() [3/3]

Member Function Documentation

◆ operator=()

◆ invert()

◆ solve()

◆ update()

◆ getInstance()

◆ resize()

◆ reset()

◆ setIdentity()

◆ negate()

◆ add()

◆ sub()

◆ mul() [1/3]

◆ mul() [2/3]

◆ mul() [3/3]

◆ div()

◆ equals()

◆ isIdentity()

◆ isSymmetric()

◆ getDot()

◆ rswap()

◆ cswap()

◆ clear()

◆ set()

◆ swap()

◆ transpose()

◆ size()

◆ capacity()

◆ empty()

◆ data() [1/2]

◆ data() [2/2]

◆ operator[]() [1/2]

◆ operator[]() [2/2]

◆ operator()() [1/2]

◆ operator()() [2/2]

◆ at() [1/2]

◆ at() [2/2]

Member Data Documentation

◆ P

◆ ws

◆ __p

◆ __n

◆ __m