1 #ifndef __JMATH__JMATRIXND__
2 #define __JMATH__JMATRIXND__
25 namespace JPP {
using namespace JMATH; }
65 public JMath <JMatrixND>,
100 matrix_type::resize(size);
102 for (
row_type row = this->begin(); row != this->end(); ++row) {
106 for (
col_type col = row->begin(); col != row->end(); ++col) {
120 double at(
size_t row,
size_t col)
const
122 return (*
this)[row][col];
133 double&
at(
size_t row,
size_t col)
135 return (*
this)[row][col];
172 for (
row_type row = this->begin(); row != this->end(); ++row) {
173 for (
col_type col = row->begin(); col != row->end(); ++col) {
189 for (
size_t i = 0; i != this->size(); ++i) {
193 for (
size_t j = 0; j != i; ++j) {
194 at(i,j) =
at(j,i) = 0.0;
209 for (
size_t i = 0; i != this->size(); ++i) {
210 for (
size_t j = 0; j != i; ++j) {
211 std::swap(
at(i,j),
at(j,i));
226 for (
row_type row = this->begin(); row != this->end(); ++row) {
227 for (
col_type col = row->begin(); col != row->end(); ++col) {
244 if (this->size() == A.size()) {
246 for (
size_t i = 0; i != this->size(); ++i) {
247 for (
size_t j = 0; j != this->size(); ++j) {
248 at(i,j) += A.
at(i,j);
253 THROW(
JException,
"JMatrixND::add() inconsistent matrix dimensions " << this->size() <<
' ' << A.size());
268 if (this->size() == A.size()) {
270 for (
size_t i = 0; i != this->size(); ++i) {
271 for (
size_t j = 0; j != this->size(); ++j) {
272 at(i,j) -= A.
at(i,j);
277 THROW(
JException,
"JMatrixND::sub() inconsistent matrix dimensions " << this->size() <<
' ' << A.size());
292 for (
row_type row = this->begin(); row != this->end(); ++row) {
293 for (
col_type col = row->begin(); col != row->end(); ++col) {
310 for (
row_type row = this->begin(); row != this->end(); ++row) {
311 for (
col_type col = row->begin(); col != row->end(); ++col) {
330 if (A.size() == B.size()) {
334 for (
size_t i = 0; i != A.size(); ++i) {
335 for (
size_t j = 0; j != A.size(); ++j) {
339 for (
size_t k = 0; k != A.size(); ++k) {
340 this->
at(i,j) += A[i][k] * B[k][j];
347 THROW(
JException,
"JMatrixND::mul() inconsistent matrix dimensions " << A.size() <<
' ' << B.size());
362 const double eps = std::numeric_limits<double>::min())
const
364 if (this->size() == A.size()) {
366 for (
const_row_type row1 = this->begin(), row2 = A.begin(); row1 != this->end(); ++row1, ++row2) {
367 for (
const_col_type col1 = row1->begin(), col2 = row2->begin(); col1 != row1->end(); ++col1, ++col2) {
368 if (fabs(*col1 - *col2) > eps) {
377 THROW(
JException,
"JMatrixND::equals() inconsistent matrix dimensions " << this->size() <<
' ' << A.size());
388 bool isIdentity(
const double eps = std::numeric_limits<double>::min())
const
390 for (
size_t i = 0; i != this->size(); ++i) {
392 if (fabs(1.0 -
at(i,i)) > eps) {
396 for (
size_t j = 0; j != i; ++j) {
397 if (fabs(
at(i,j)) > eps || fabs(
at(j,i)) > eps) {
420 for (
size_t i = 0; i != A.size(); ++i) {
439 for (
size_t i = 0; i != A.size(); ++i) {
440 if (A[i][i] <= 0.0) {
460 for (
size_t i = 0; i != A.size(); ++i) {
485 for (
row_type row = matrix.begin(); row != matrix.end(); ++row) {
486 for (
col_type col = row->begin(); col != row->end(); ++col) {
504 int size = matrix.size();
508 for (
const_row_type row = matrix.begin(); row != matrix.end(); ++row) {
509 for (
const_col_type col = row->begin(); col != row->end(); ++col) {
533 for (
const_col_type col = row->begin(); col != row->end(); ++col) {
534 out <<
FIXED(10,2) << *col <<
' ';
583 for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {
590 for (row = root; ++row != this->end(); ) {
591 if (fabs((*row)[k]) > fabs(val)) {
603 iter_swap(root, pivot);
611 for (row = root; ++row != this->end(); ) {
613 del = (*row)[k] /= val;
615 for (col = row->begin() + k, cursor = root->begin() + k; ++col != row->end(); )
616 (*col) -= del * (*(++cursor));
660 dot += (*yl) * (*col) * (*yr);
double & at(size_t row, size_t col)
Get matrix element.
double at(size_t row, size_t col) const
Get matrix element.
Interface for binary output.
friend JReader & operator>>(JReader &in, JMatrixND &matrix)
Read matrix from input.
Auxiliary base class for aritmetic operations of derived class types.
std::vector< double >::iterator col_type
int getSign() const
Get sign corresponding to permutation matrix.
double getDot(const JNeutrinoDirection &first, const JNeutrinoDirection &second)
Dot product.
JMatrixND & transpose()
Transpose.
#define THROW(JException_t, A)
Marco for throwing exception with std::ostream compatible message.
Auxiliary data structure for floating point format specification.
const bool isPositiveSemiDefinite() const
Test positive semi-definiteness.
bool isPositiveDefinite() const
Test positive definiteness.
const JMatrixHelper & set(const JMatrixND &A)
Set matrix helper.
void resize(const size_t size)
Resize matrix.
JMatrixND & add(const JMatrixND &A)
Matrix addition.
std::vector< std::pair< int, int > > JPermutationMatrix
Type definition of permutation matrix.
I/O formatting auxiliaries.
matrix_type::const_iterator const_row_type
JMatrixND & setIdentity()
Set to identity matrix.
JMatrixND & reset()
Set matrix to the null matrix.
matrix_type::iterator row_type
std::vector< double >::const_iterator const_col_type
Template definition of auxiliary base class for comparison of data structures.
friend std::ostream & operator<<(std::ostream &out, const JMatrixND &A)
Print ASCII formatted output.
friend JWriter & operator<<(JWriter &out, const JMatrixND &matrix)
Write matrix to output.
Interface for binary input.
JMatrixND(const size_t size)
Constructor.
JMatrixND & negate()
Negate matrix.
Auxiliary class for matrix decomposition.
void decompose(const JMatrixND &A)
LU decomposition.
bool equals(const JMatrixND &A, const double eps=std::numeric_limits< double >::min()) const
Equality.
Exception for division by zero.
static data_type & getInstance()
Get unique instance of template class.
Base class for data structures with artithmetic capabilities.
JMatrixND & sub(const JMatrixND &A)
Matrix subtraction.
double getDeterminant() const
Get determinant of matrix.
JMatrixND & mul(const JMatrixND &A, const JMatrixND &B)
Matrix multiplication.
JMatrixND & div(const double factor)
Scale matrix.
JMatrixND & mul(const double factor)
Scale matrix.
bool isIdentity(const double eps=std::numeric_limits< double >::min()) const
Test identity.
Type definition of a NxN matrix.
Auxiliary class to temporarily modify format specifications.
double operator()(size_t row, size_t col) const
Get matrix element.
JMatrixND()
Default constructor.
std::vector< std::vector< double > > matrix_type
double & operator()(size_t row, size_t col)
Get matrix element.