1 #ifndef __JMATH__JMATRIXNS__
2 #define __JMATH__JMATRIXNS__
17 namespace JPP {
using namespace JMATH; }
79 size_t i, row, col,
root;
95 P.resize(this->
size());
104 for (row = root; ++row != this->
size(); ) {
105 if (fabs((*
this)(row, root)) > fabs(value)) {
106 value = (*this)(row,
root);
116 this->
rswap(root, pivot);
121 for (row = root + 1; row + 4 <= this->
size(); row += 4) {
123 p0 = (*this)[row + 0] +
root;
124 p1 = (*this)[row + 1] +
root;
125 p2 = (*this)[row + 2] +
root;
126 p3 = (*this)[row + 3] +
root;
128 v0 = (*p0++ /= value);
129 v1 = (*
p1++ /= value);
130 v2 = (*
p2++ /= value);
131 v3 = (*
p3++ /= value);
133 c0 = (*this)[
root] + root + 1;
135 for (
i = root + 1;
i + 4 <= this->
size();
i += 4, p0 += 4,
p1 += 4,
p2 += 4,
p3 += 4, c0 += 4) {
136 p0[0] -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
137 p1[0] -= v1 * c0[0];
p1[1] -= v1 * c0[1];
p1[2] -= v1 * c0[2];
p1[3] -= v1 * c0[3];
138 p2[0] -= v2 * c0[0];
p2[1] -= v2 * c0[1];
p2[2] -= v2 * c0[2];
p2[3] -= v2 * c0[3];
139 p3[0] -= v3 * c0[0];
p3[1] -= v3 * c0[1];
p3[2] -= v3 * c0[2];
p3[3] -= v3 * c0[3];
142 for ( ;
i != this->
size(); ++
i, ++p0, ++
p1, ++
p2, ++
p3, ++c0) {
150 for ( ; row != this->
size(); ++row) {
152 p0 = (*this)[row + 0] +
root;
154 v0 = (*p0++ /= value);
156 c0 = (*this)[
root] + root + 1;
158 for (
i = root + 1;
i + 4 <= this->
size();
i += 4, p0 += 4, c0 += 4) {
159 *p0 -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
162 for ( ;
i != this->
size(); ++
i, ++p0, ++c0) {
173 for (
i = 0;
i != this->
size(); ++
i) {
182 for ( ; col + 4 <= this->
size(); col += 4) {
184 for (row = 0; row != this->
size(); ++row) {
203 for (
i = 0, c0 = (*
this)[row];
i != row; ++c0, ++p0, ++
p1, ++
p2, ++
p3, ++
i) {
210 A[col + 0][row] = v0;
211 A[col + 1][row] = v1;
212 A[col + 2][row] = v2;
213 A[col + 3][row] = v3;
217 for ( ; col != this->
size(); ++col) {
219 for (row = 0; row != this->
size(); ++row) {
229 for (
i = 0, c0 = (*
this)[row] +
i;
i != row; ++c0, ++p0, ++
i) {
241 for ( ; col + 4 <= this->
size(); col += 4) {
243 for (
int row = this->
size() - 1; row >= 0; --row) {
245 p0 =
A[col + 0] + row;
246 p1 =
A[col + 1] + row;
247 p2 =
A[col + 2] + row;
248 p3 =
A[col + 3] + row;
260 for (
i = row + 1, c0 = (*
this)[row] +
i;
i != this->
size(); ++c0, ++p0, ++
p1, ++
p2, ++
p3, ++
i) {
267 w0 = 1.0 / (*this)[row][row];
269 A[col + 0][row] = v0 * w0;
270 A[col + 1][row] = v1 * w0;
271 A[col + 2][row] = v2 * w0;
272 A[col + 3][row] = v3 * w0;
276 for ( ; col != this->
size(); ++col) {
278 for (
int row = this->
size() - 1; row >= 0; --row) {
286 for (
i = row + 1, c0 = (*
this)[row] +
i;
i != this->
size(); ++c0, ++p0, ++
i) {
290 w0 = 1.0 / (*this)[row][row];
292 A[col][row] = v0 * w0;
307 template<
class JVectorND_t>
326 P.resize(this->
size());
335 for (row = root; ++row != this->
size(); ) {
336 if (fabs((*
this)(row, root)) > fabs(value)) {
337 value = (*this)(row,
root);
347 this->
rswap(root, pivot);
352 for (row = root + 1; row + 4 <= this->
size(); row += 4) {
354 p0 = (*this)[row + 0] +
root;
355 p1 = (*this)[row + 1] +
root;
356 p2 = (*this)[row + 2] +
root;
357 p3 = (*this)[row + 3] +
root;
359 v0 = (*p0++ /= value);
360 v1 = (*
p1++ /= value);
361 v2 = (*
p2++ /= value);
362 v3 = (*
p3++ /= value);
364 c0 = (*this)[
root] + root + 1;
366 for (
i = root + 1;
i + 4 <= this->
size();
i += 4, p0 += 4,
p1 += 4,
p2 += 4,
p3 += 4, c0 += 4) {
367 p0[0] -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
368 p1[0] -= v1 * c0[0];
p1[1] -= v1 * c0[1];
p1[2] -= v1 * c0[2];
p1[3] -= v1 * c0[3];
369 p2[0] -= v2 * c0[0];
p2[1] -= v2 * c0[1];
p2[2] -= v2 * c0[2];
p2[3] -= v2 * c0[3];
370 p3[0] -= v3 * c0[0];
p3[1] -= v3 * c0[1];
p3[2] -= v3 * c0[2];
p3[3] -= v3 * c0[3];
373 for ( ;
i != this->
size(); ++
i, ++p0, ++
p1, ++
p2, ++
p3, ++c0) {
381 for ( ; row != this->
size(); ++row) {
383 p0 = (*this)[row + 0] +
root;
385 v0 = (*p0++ /= value);
387 c0 = (*this)[
root] + root + 1;
389 for (
i = root + 1;
i + 4 <= this->
size();
i += 4, p0 += 4, c0 += 4) {
390 *p0 -= v0 * c0[0]; p0[1] -= v0 * c0[1]; p0[2] -= v0 * c0[2]; p0[3] -= v0 * c0[3];
393 for ( ;
i != this->
size(); ++
i, ++p0, ++c0) {
401 for (row = 0; row != this->
size(); ++row) {
407 for (
i = 0, c0 = (*
this)[row] +
i;
i != row; ++c0, ++
i) {
416 for (
int row = this->
size() - 1; row >= 0; --row) {
420 for (
i = row + 1, c0 = (*
this)[row] +
i;
i < this->
size(); ++c0, ++
i) {
424 u[row] = v0 / (*this)[row][row];
446 void update(
const size_t k,
const double value)
448 if (k >= this->
size()) {
455 const double del = (*this)(
k,
k);
456 const double din = 1.0 / (1.0 + value * del);
458 double*
root = (*this)[
k];
460 for (row = 0; row != this->
size(); ++row) {
464 const double val = (*this)(row,
k);
466 for (i = 0, col = (*
this)[row]; i != this->
size(); ++
i, ++col) {
468 (*col) -= val * root[
i] * value * din;
472 (*this)(row,
k) *= din;
476 for (i = 0, col = root; i != this->
size(); ++
i, ++col) {
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.
JMatrixNS(const JMatrixND &A)
Contructor.
JMatrixND_t & getInstance()
Get work space.
void update(const size_t k, const double value)
Update inverted matrix at given diagonal element.
JMatrixNS()
Default constructor.
std::vector< int > P
permutation matrix
void set(const JMatrixND_t &A)
Set matrix.
JMatrixNS & operator=(const JMatrixNS &A)
Assignment operator.
void rswap(size_t r1, size_t r2)
void invert()
Invert matrix according LDU decomposition.
size_t size() const
Get dimension of matrix.
Exception for division by zero.
void solve(JVectorND_t &u)
Get solution of equation A x = b.
void swap(JMatrixND_t &A)
Swap matrices.
Exception for accessing an index in a collection that is outside of its range.
source $JPP_DIR setenv csh $JPP_DIR &dev null eval JShellParser o a A
JMatrixNS(const size_t size)
Constructor (identity matrix).
do JPlot2D f $WORKDIR canberra[${EMITTER}] root