SHS/pip
4
3
Fork 0
Code Issues 67 Pull Requests 1 Actions Releases Wiki Activity

PIVector2D and PIMathMatrix

Browse Source 
completed and optimized

git-svn-id: svn://db.shs.com.ru/pip@657 12ceb7fc-bf1f-11e4-8940-5bc7170c53b5
This commit is contained in:
Бычков Андрей
2018-11-13 13:14:09 +00:00
parent e35313dc1a
commit 95983ca9a4
3 changed files with 119 additions and 94 deletions
Download Patch File Download Diff File Expand all files Collapse all files

99
src_main/math/pimathmatrix.h
View File

@@ -30,16 +30,16 @@
template<typename T>
inline bool _PIMathMatrixNullCompare(const T v) {
return (piAbs(v) < T(1E-100));
return (piAbs(v) < T(1E-200));
}
template<>
inline bool _PIMathMatrixNullCompare<complexf >(const complexf v) {
return (abs(v) < float(1E-100));
return (abs(v) < float(1E-200));
}
template<>
inline bool _PIMathMatrixNullCompare<complexd >(const complexd v) {
return (abs(v) < double(1E-100));
return (abs(v) < double(1E-200));
}
@@ -135,11 +135,6 @@ public:
if (ok != 0) *ok = false;
return *this;
}
for (uint j = 0; j < Cols; ++j) if (smat.m[j][i] != 0) ndet = false;
if (ndet) {
if (ok != 0) *ok = false;
return *this;
}
}
for (uint i = 0; i < Cols; ++i) {
crow = i;
@@ -150,7 +145,7 @@ public:
for (uint k = i; k < Cols; ++k) smat.m[k][j] -= mul * smat.m[k][i];
}
if (i < Cols - 1) {
if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
if (fabs(smat.m[i+1][i+1]) < Type(1E-200)) {
if (ok != 0) *ok = false;
return *this;
}
@@ -177,11 +172,6 @@ public:
if (ok != 0) *ok = false;
return *this;
}
for (uint j = 0; j < Cols; ++j) if (smat.m[j][i] != 0) ndet = false;
if (ndet) {
if (ok != 0) *ok = false;
return *this;
}
}
for (uint i = 0; i < Cols; ++i) {
crow = i;
@@ -197,7 +187,7 @@ public:
}
//cout << i << endl << smat << endl;
if (i < Cols - 1) {
if (fabs(smat.m[i+1][i+1]) < Type(1E-100)) {
if (fabs(smat.m[i+1][i+1]) < Type(1E-200)) {
if (ok != 0) *ok = false;
return *this;
}
@@ -337,26 +327,25 @@ class PIP_EXPORT PIMathMatrix : public PIVector2D<Type> {
typedef PIVector2D<Type> _V2D;
typedef PIMathMatrix<Type> _CMatrix;
typedef PIMathVector<Type> _CMCol;
typedef PIMathVector<Type> _CMRow;
public:
PIMathMatrix(const uint cols = 3, const uint rows = 3) {resize(cols, rows);}
PIMathMatrix(const uint cols, const uint rows, const PIVector<Type> & val) {resize(cols, rows); int i=0; PIMM_FOR_I(c, r) (*this)[r][c] = val[i++];}
PIMathMatrix(const PIVector<PIVector<Type> > & val) {_V2D::cols_ = _V2D::rows_ = 0; if(!val.isEmpty()) {resize(val[0].size(), val.size()); PIMM_FOR_I(c, r) (*this)[r][c] = val[r][c];}}
PIMathMatrix(const PIVector2D<Type> & val) {_V2D::cols_ = _V2D::rows_ = 0; if(!val.isEmpty()) {resize(val.cols(), val.rows()); PIMM_FOR_I(c, r) (*this)[r][c] = val[r][c];}}
PIMathMatrix(const uint cols, const uint rows, const PIVector<Type> & val) {resize(cols, rows); int i=0; PIMM_FOR_I(c, r) _V2D::element(r, c) = val[i++];}
PIMathMatrix(const PIVector<PIVector<Type> > & val) {_V2D::cols_ = _V2D::rows_ = 0; if(!val.isEmpty()) {resize(val[0].size(), val.size()); PIMM_FOR_I(c, r) _V2D::element(r, c) = val[r][c];}}
PIMathMatrix(const PIVector2D<Type> & val) {_V2D::cols_ = _V2D::rows_ = 0; if(!val.isEmpty()) {resize(val.cols(), val.rows()); PIMM_FOR_I(c, r) _V2D::element(r, c) = val.element(r, c);}}
static _CMatrix identity(const uint cols, const uint rows) {_CMatrix tm(cols, rows); for (uint r = 0; r < rows; ++r) for (uint c = 0; c < cols; ++c) tm[r][c] = (c == r ? Type(1) : Type(0)); return tm;}
static _CMatrix identity(const uint cols, const uint rows) {_CMatrix tm(cols, rows); for (uint r = 0; r < rows; ++r) for (uint c = 0; c < cols; ++c) tm.element(r, c) = (c == r ? Type(1) : Type(0)); return tm;}
static _CMatrix matrixRow(const PIMathVector<Type> & val) {return _CMatrix(val.size(), 1, val.toVector());}
static _CMatrix matrixCol(const PIMathVector<Type> & val) {return _CMatrix(1, val.size(), val.toVector());}
_CMatrix & resize(const uint cols, const uint rows, const Type & new_value = Type()) {_V2D::_resizeRaw(rows, cols); PIMM_FOR_A(i) _V2D::mat[i] = new_value; return *this;}
_CMatrix & swapCols(uint r0, uint r1) {Type t; PIMM_FOR_C(i) {t = (*this)[i][r0]; (*this)[i][r0] = (*this)[i][r1]; (*this)[i][r1] = t;} return *this;}
_CMatrix & swapRows(uint c0, uint c1) {Type t; PIMM_FOR_R(i) {t = (*this)[c0][i]; (*this)[c0][i] = (*this)[c1][i]; (*this)[c1][i] = t;} return *this;}
_CMatrix & swapCols(uint r0, uint r1) {PIMM_FOR_C(i) {piSwap(_V2D::element(i, r0), _V2D::element(i, r1));} return *this;}
_CMatrix & swapRows(uint c0, uint c1) {PIMM_FOR_R(i) {piSwap(_V2D::element(c0, i), _V2D::element(c1, i));} return *this;}
_CMatrix & fill(const Type & v) {PIMM_FOR_A(i) _V2D::mat[i] = v; return *this;}
bool isSquare() const {return _V2D::cols_ == _V2D::rows_;}
bool isIdentity() const {PIMM_FOR(c, r) if ((c == r) ? (*this)[c][r] != Type(1) : (*this)[c][r] != Type(0)) return false; return true;}
bool isIdentity() const {PIMM_FOR(c, r) if ((c == r) ? _V2D::element(r, c) != Type(1) : _V2D::element(r, c) != Type(0)) return false; return true;}
bool isNull() const {PIMM_FOR_A(i) if (_V2D::mat[i] != Type(0)) return false; return true;}
_CMatrix & operator =(const PIVector<PIVector<Type> > & v) {*this = PIVector2D<Type>(v); return *this;}
_CMatrix & operator =(const PIVector<PIVector<Type> > & v) {*this = _CMatrix(v); return *this;}
bool operator ==(const _CMatrix & sm) const {PIMM_FOR_A(i) if (_V2D::mat[i] != sm.mat[i]) return false; return true;}
bool operator !=(const _CMatrix & sm) const {return !(*this == sm);}
void operator +=(const _CMatrix & sm) {PIMM_FOR_A(i) _V2D::mat[i] += sm.mat[i];}
@@ -380,7 +369,7 @@ public:
for (uint c = 0; c < _V2D::cols_; ++c)
for (uint r = 0; r < _V2D::rows_; ++r)
if (r == c)
ret *= m[r][c];
ret *= m.element(r, c);
return ret;
}
@@ -391,7 +380,7 @@ public:
return ret;
}
for (uint i = 0; i < _V2D::cols_; ++i) {
ret += (*this)[i][i];
ret += _V2D::element(i, i);
}
if (ok != 0) *ok = true;
return ret;
@@ -408,12 +397,7 @@ public:
Type mul;
for (uint i = 0; i < _V2D::cols_; ++i) {
ndet = true;
for (uint j = 0; j < _V2D::rows_; ++j) if (smat.mat[i][j] != 0) ndet = false;
if (ndet) {
if (ok != 0) *ok = false;
return *this;
}
for (uint j = 0; j < _V2D::cols_; ++j) if (smat.mat[j][i] != 0) ndet = false;
for (uint j = 0; j < _V2D::rows_; ++j) if (smat.element(i, j) != 0) ndet = false;
if (ndet) {
if (ok != 0) *ok = false;
return *this;
@@ -421,14 +405,14 @@ public:
}
for (uint i = 0; i < _V2D::cols_; ++i) {
crow = i;
while (smat.mat[i][i] == Type(0))
while (smat.element(i, i) == Type(0))
smat.swapRows(i, ++crow);
for (uint j = i + 1; j < _V2D::rows_; ++j) {
mul = smat.mat[i][j] / smat.mat[i][i];
for (uint k = i; k < _V2D::cols_; ++k) smat.mat[k][j] -= mul * smat.mat[k][i];
mul = smat.element(i, j) / smat.element(i, i);
for (uint k = i; k < _V2D::cols_; ++k) smat.element(k, j) -= mul * smat.element(k, i);
}
if (i < _V2D::cols_ - 1) {
if (fabs(smat.mat[i+1][i+1]) < Type(1E-100)) {
if (_PIMathMatrixNullCompare(smat.element(i+1, i+1))) {
if (ok != 0) *ok = false;
return *this;
}
@@ -450,12 +434,7 @@ public:
Type mul, iddiv;
for (uint i = 0; i < _V2D::cols_; ++i) {
ndet = true;
for (uint j = 0; j < _V2D::rows_; ++j) if (smat[i][j] != Type(0)) ndet = false;
if (ndet) {
if (ok != 0) *ok = false;
return *this;
}
for (uint j = 0; j < _V2D::cols_; ++j) if (smat[j][i] != Type(0)) ndet = false;
for (uint j = 0; j < _V2D::rows_; ++j) if (smat.element(i, j) != Type(0)) ndet = false;
if (ndet) {
if (ok != 0) *ok = false;
return *this;
@@ -463,35 +442,35 @@ public:
}
for (uint i = 0; i < _V2D::cols_; ++i) {
crow = i;
while (smat[i][i] == Type(0)) {
while (smat.element(i, i) == Type(0)) {
++crow;
smat.swapRows(i, crow);
mtmp.swapRows(i, crow);
if (sv != 0) sv->swap(i, crow);
}
for (uint j = i + 1; j < _V2D::rows_; ++j) {
mul = smat[i][j] / smat[i][i];
for (uint k = i; k < _V2D::cols_; ++k) smat[k][j] -= mul * smat[k][i];
for (uint k = 0; k < _V2D::cols_; ++k) mtmp[k][j] -= mul * mtmp[k][i];
mul = smat.element(i, j) / smat.element(i, i);
for (uint k = i; k < _V2D::cols_; ++k) smat.element(k, j) -= mul * smat.element(k, i);
for (uint k = 0; k < _V2D::cols_; ++k) mtmp.element(k, j) -= mul * mtmp.element(k, i);
if (sv != 0) (*sv)[j] -= mul * (*sv)[i];
}
//cout << i << endl << smat << endl;
if (i < _V2D::cols_ - 1) {
if (_PIMathMatrixNullCompare(smat[i+1][i+1])) {
if (_PIMathMatrixNullCompare(smat.element(i+1, i+1))) {
if (ok != 0) *ok = false;
return *this;
}
}
iddiv = smat[i][i];
for (uint j = i; j < _V2D::cols_; ++j) smat[j][i] /= iddiv;
for (uint j = 0; j < _V2D::cols_; ++j) mtmp[j][i] /= iddiv;
iddiv = smat.element(i, i);
for (uint j = i; j < _V2D::cols_; ++j) smat.element(j, i) /= iddiv;
for (uint j = 0; j < _V2D::cols_; ++j) mtmp.element(j, i) /= iddiv;
if (sv != 0) (*sv)[i] /= iddiv;
}
for (uint i = _V2D::cols_ - 1; i > 0; --i) {
for (uint j = 0; j < i; ++j) {
mul = smat[i][j];
smat[i][j] -= mul;
for (uint k = 0; k < _V2D::cols_; ++k) mtmp[k][j] -= mtmp[k][i] * mul;
mul = smat.element(i, j);
smat.element(i, j) -= mul;
for (uint k = 0; k < _V2D::cols_; ++k) mtmp.element(k, j) -= mul * mtmp.element(k, i);
if (sv != 0) (*sv)[j] -= mul * (*sv)[i];
}
}
@@ -500,7 +479,7 @@ public:
return *this;
}
_CMatrix inverted(bool * ok = 0) {_CMatrix tm(*this); tm.invert(ok); return tm;}
_CMatrix transposed() {_CMatrix tm(_V2D::rows_, _V2D::cols_); PIMM_FOR(c, r) tm[c][r] = (*this)[r][c]; return tm;}
_CMatrix transposed() {_CMatrix tm(_V2D::rows_, _V2D::cols_); PIMM_FOR(c, r) tm.element(c, r) = _V2D::element(r, c); return tm;}
private:
// size_t rows_, cols_;
@@ -510,11 +489,11 @@ private:
#ifdef PIP_STD_IOSTREAM
template<typename Type>
inline std::ostream & operator <<(std::ostream & s, const PIMathMatrix<Type> & m) {s << '{'; for (uint r = 0; r < m.rows(); ++r) { for (uint c = 0; c < m.cols(); ++c) { s << m[c][r]; if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << std::endl << ' ';} s << '}'; return s;}
inline std::ostream & operator <<(std::ostream & s, const PIMathMatrix<Type> & m) {s << '{'; for (uint r = 0; r < m.rows(); ++r) { for (uint c = 0; c < m.cols(); ++c) { s << m.element(r, c); if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << std::endl << ' ';} s << '}'; return s;}
#endif
template<typename Type>
inline PICout operator <<(PICout s, const PIMathMatrix<Type> & m) {s << '{'; for (uint r = 0; r < m.rows(); ++r) { for (uint c = 0; c < m.cols(); ++c) { s << m[r][c]; if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << PICoutManipulators::NewLine << ' ';} s << '}'; return s;}
inline PICout operator <<(PICout s, const PIMathMatrix<Type> & m) {s << '{'; for (uint r = 0; r < m.rows(); ++r) { for (uint c = 0; c < m.cols(); ++c) { s << m.element(r, c); if (c < m.cols() - 1 || r < m.rows() - 1) s << ", ";} if (r < m.rows() - 1) s << PICoutManipulators::NewLine << ' ';} s << '}'; return s;}
/// Multiply matrices {CR x Rows0} on {Cols1 x CR}, result is {Cols1 x Rows0}
template<typename Type>
@@ -528,8 +507,8 @@ inline PIMathMatrix<Type> operator *(const PIMathMatrix<Type> & fm,
for (uint i = 0; i < cols1; ++i) {
t = Type(0);
for (uint k = 0; k < cr; ++k)
t += fm[j][k] * sm[k][i];
tm[j][i] = t;
t += fm.element(j, k) * sm.element(k, i);
tm.element(j, i) = t;
}
}
return tm;
@@ -546,7 +525,7 @@ inline PIMathVector<Type> operator *(const PIMathMatrix<Type> & fm,
for (uint i = 0; i < r; ++i) {
t = Type(0);
for (uint j = 0; j < c; ++j)
t += fm[j][i] * sv[j];
t += fm.element(j, i) * sv[j];
tv[i] = t;
}
return tv;
@@ -558,7 +537,7 @@ typedef PIMathMatrix<double> PIMathMatrixd;
template<typename T>
PIMathMatrix<complex<T> > hermitian(const PIMathMatrix<complex<T> > & m) {
PIMathMatrix<complex<T> > ret(m);
for (uint r = 0; r < ret.rows(); ++r) for (uint c = 0; c < ret.cols(); ++c) ret[r][c].imag(-(ret[r][c].imag()));
for (uint r = 0; r < ret.rows(); ++r) for (uint c = 0; c < ret.cols(); ++c) ret.element(r, c).imag(-(ret.element(r, c).imag()));
return ret.transposed();
}