1 #ifndef __JMATH__JMATRIXNS__
2 #define __JMATH__JMATRIXNS__
15 namespace JPP {
using namespace JMATH; }
72 for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {
79 for (row = root; ++row != this->end(); ) {
80 if (fabs((*row)[k]) > fabs(val)) {
92 std::iter_swap(root, pivot);
94 P.push_back(std::make_pair(std::distance(this->begin(),root),
95 std::distance(this->begin(),pivot)));
100 for (row = root; ++row != this->end(); ) {
102 del = (*row)[k] /= val;
104 for (col = row->begin() + k, cursor = root->begin() + k; ++col != row->end(); ) {
105 (*col) -= del * (*(++cursor));
113 for (k = 0; k != this->size(); ++k) {
115 double& a = this->
at(k,k);
127 for (root = this->begin(), k = 0; root != this->end(); ++root, ++k) {
131 for (col = pivot, j = k + 1; ++col != root->end(); ++j) {
135 for (cursor = pivot, row = root; ++cursor != col; ) {
136 (*col) -= (*cursor) * (*(++row))[j];
139 (*col) *= this->
at(j,j);
146 for (i = size(); i != 0; ) {
148 root = this->begin() + --i;
150 for (j = i; j != 0; ) {
152 double& a = (*root)[--j];
156 for (k = j + 1, cursor = root->begin() + k, row = this->begin() + k; k != i; ++cursor, ++row, ++k) {
157 a -= (*cursor) * (*row)[j];
165 for (root = this->begin(), i = 0; root != this->end(); ++root, ++i) {
169 for (col = root->begin(), j = 0; j != i; ++col, ++j) {
173 for (row = root, cursor = root->begin() + i; ++cursor != root->end(); ) {
174 (*col) += (*cursor) * (*(++row))[j];
178 for (j = i, col = root->begin() + j; j != this->size() - 1; ++col, ++j) {
179 for (k = j + 1, row = this->begin() + k, cursor = root->begin() + k; k != this->size(); ++k, ++cursor, ++row) {
180 (*col) += (*(cursor)) * (*(row))[j];
186 for (JPermutationMatrix::reverse_iterator p = P.rbegin(); p != P.rend(); ++p) {
187 for (row = this->begin(); row != this->end(); ++row) {
188 std::swap((*row)[p->first], (*row)[p->second]);
208 void update(
const unsigned int k,
const double extra)
210 if (k >= this->size()) {
218 const double del = this->
at(k,k);
219 const double din = 1.0 / (1.0 + extra * del);
221 root = this->begin() + k;
223 for (i = 0, row = this->begin(); row != this->end(); ++i, ++row) {
227 const double val = (*row)[k];
229 for (j = 0, col = row->begin(); col != row->end(); ++j, ++col) {
232 (*col) -= val * (*root)[j] * extra * din;
240 for (j = 0, col = root->begin(); col != root->end(); ++j, ++col) {
247 (*root)[k] = del * din;
double at(size_t row, size_t col) const
Get matrix element.
std::vector< double >::iterator col_type
JMatrixNS(const JMatrixND &A)
Contructor.
void update(const unsigned int k, const double extra)
Update inverted matrix at given diagonal element.
JMatrixNS(const unsigned int size)
Constructor (identity matrix).
matrix_type::iterator row_type
Exception for division by zero.
JMatrixNS()
Default constructor.
Exception for accessing an index in a collection that is outside of its range.
void invert()
Invert matrix.