From 95983ca9a4d0d808049a8835aca155b6674104b0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=D0=91=D1=8B=D1=87=D0=BA=D0=BE=D0=B2=20=D0=90=D0=BD=D0=B4?= =?UTF-8?q?=D1=80=D0=B5=D0=B9?= Date: Tue, 13 Nov 2018 13:14:09 +0000 Subject: [PATCH] PIVector2D and PIMathMatrix completed and optimized git-svn-id: svn://db.shs.com.ru/pip@657 12ceb7fc-bf1f-11e4-8940-5bc7170c53b5 --- src_main/containers/pivector2d.h | 110 ++++++++++++++++++++++--------- src_main/math/pimathmatrix.h | 99 +++++++++++----------------- src_main/math/pimathsolver.cpp | 4 +- 3 files changed, 119 insertions(+), 94 deletions(-) diff --git a/src_main/containers/pivector2d.h b/src_main/containers/pivector2d.h index 3f929ac4..eb40dcfa 100644 --- a/src_main/containers/pivector2d.h +++ b/src_main/containers/pivector2d.h @@ -73,63 +73,110 @@ public: class Row { friend class PIVector2D; private: - inline Row(PIVector2D * p, size_t row) : p_(p) {st_ = p_->cols_ * row;} - PIVector2D * p_; - size_t st_; + inline Row(PIVector2D * p, size_t row) : p_(&(p->mat)) {st_ = p->cols_ * row; sz_ = p->cols_;} + PIVector * p_; + size_t st_, sz_; public: - inline size_t size() const {return p_->cols_;} - inline T & operator [](size_t index) {return p_->mat[st_ + index];} - inline const T & operator [](size_t index) const {return p_->mat[st_ + index];} - inline T * data(size_t index = 0) {return p_->mat.data(st_ + index);} - inline const T * data(size_t index = 0) const {return p_->mat.data(st_ + index);} + inline size_t size() const {return sz_;} + inline T & operator [](size_t index) {return (*p_)[st_ + index];} + inline const T & operator [](size_t index) const {return (*p_)[st_ + index];} + inline T * data(size_t index = 0) {return p_->data(st_ + index);} + inline const T * data(size_t index = 0) const {return p_->data(st_ + index);} inline Row & operator =(const Row & other) { if (p_ == other.p_ && st_ == other.st_) return *this; - size_t sz = piMin(p_->cols_, other.p_->cols_); - p_->copyRow(st_, other.data(), sz); + size_t sz = piMin(sz_, other.sz_); + copyRow(st_, other.data(), sz, p_); return *this; } inline Row & operator =(const PIVector & other) { - size_t sz = piMin(p_->cols_, other.size()); - p_->copyRow(st_, other.data(), sz); + size_t sz = piMin(sz, other.size()); + copyRow(st_, other.data(), sz, p_); return *this; } - inline PIVector toVector() const {return PIVector(p_->mat.data(st_), p_->cols_);} + inline PIVector toVector() const {return PIVector(p_->data(st_), sz_);} + }; + + class Col { + friend class PIVector2D; + private: + inline Col(PIVector2D * p, size_t row) : p_(&(p->mat)) {step_ = p->cols_; row_ = row; sz_ = p->rows_;} + PIVector * p_; + size_t step_, row_, sz_; + public: + inline size_t size() const {return sz_;} + inline T & operator [](size_t index) {return (*p_)[index * step_ + row_];} + inline const T & operator [](size_t index) const {return (*p_)[index * step_ + row_];} + inline T * data(size_t index = 0) {return p_->data(index * step_ + row_);} + inline const T * data(size_t index = 0) const {return p_->data(index * step_ + row_);} + inline Col & operator =(const Col & other) { + if (p_ == other.p_ && row_ == other.row_) return *this; + size_t sz = piMin(sz_, other.sz_); + for (int i=0; i & other) { + size_t sz = piMin(sz_, other.size()); + for (int i=0; i toVector() const { + PIVector ret; + ret.reserve(sz_); + for (size_t i=0; i; private: - inline RowConst(const PIVector2D * p, size_t row) : p_(p) {st_ = p_->cols_ * row;} - const PIVector2D * p_; - size_t st_; + inline RowConst(const PIVector2D * p, size_t row) : p_(&(p->mat)) {st_ = p->cols_ * row; sz_ = p->cols_;} + const PIVector * p_; + size_t st_, sz_; public: - inline size_t size() const {return p_->cols_;} - inline const T & operator [](size_t index) const {return p_->mat[st_ + index];} - inline const T * data(size_t index = 0) const {return p_->mat.data(st_ + index);} - inline PIVector toVector() const {return PIVector(p_->mat.data(st_), p_->cols_);} + inline size_t size() const {return sz_;} + inline const T & operator [](size_t index) const {return (*p_)[st_ + index];} + inline const T * data(size_t index = 0) const {return p_->data(st_ + index);} + inline PIVector toVector() const {return PIVector(p_->data(st_), sz_);} }; + class ColConst { + friend class PIVector2D; + private: + inline ColConst(const PIVector2D * p, size_t row) : p_(&(p->mat)) {step_ = p->cols_; row_ = row; sz_ = p->rows_;} + const PIVector * p_; + size_t step_, row_, sz_; + public: + inline size_t size() const {return p_->rows_;} + inline const T & operator [](size_t index) const {return (*p_)[index * step_ + row_];} + inline const T * data(size_t index = 0) const {return p_->data(index * step_ + row_);} + inline PIVector toVector() const { + PIVector ret; + ret.reserve(sz_); + for (int i=0; i col(size_t index) { - PIVector ret; - ret.reserve(rows_); - for (int i=0; i & setRow(size_t row, const Row & other) { size_t sz = piMin(cols_, other.p_->cols_); - copyRow(cols_ * row, other.data(), sz); + copyRow(cols_ * row, other.data(), sz, &mat); return *this; } inline PIVector2D & setRow(size_t row, const PIVector & other) { size_t sz = piMin(cols_, other.size()); - copyRow(cols_ * row, other.data(), sz); + copyRow(cols_ * row, other.data(), sz, &mat); return *this; } @@ -156,12 +203,11 @@ public: } protected: - inline void copyRow(size_t start, const T * data, size_t size) { + inline void copyRow(size_t start, const T * data, size_t size, PIVector * dst) { for (size_t i = 0; i < size; i++) - mat[start + i] = data[i]; + (*dst)[start + i] = data[i]; } -//private: size_t rows_, cols_; PIVector mat; }; @@ -187,7 +233,7 @@ inline PICout operator <<(PICout s, const PIVector2D & v) { } #define __PIVECTOR2D_SIMPLE_TYPE__(T) \ - template<> inline void PIVector2D::copyRow(size_t start, const T * data, size_t size) {memcpy(mat.data(start), data, size * sizeof(T));} \ + template<> inline void PIVector2D::copyRow(size_t start, const T * data, size_t size, PIVector * dst) {memcpy(dst->data(start), data, size * sizeof(T));} \ template<> inline PIVector2D & PIVector2D::_resizeRaw(size_t r, size_t c) {rows_ = r; cols_ = c; mat._resizeRaw(r*c); return *this;} __PIVECTOR2D_SIMPLE_TYPE__(bool) diff --git a/src_main/math/pimathmatrix.h b/src_main/math/pimathmatrix.h index 9382bfd6..eeae390e 100644 --- a/src_main/math/pimathmatrix.h +++ b/src_main/math/pimathmatrix.h @@ -30,16 +30,16 @@ template inline bool _PIMathMatrixNullCompare(const T v) { - return (piAbs(v) < T(1E-100)); + return (piAbs(v) < T(1E-200)); } template<> inline bool _PIMathMatrixNullCompare(const complexf v) { - return (abs(v) < float(1E-100)); + return (abs(v) < float(1E-200)); } template<> inline bool _PIMathMatrixNullCompare(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 { typedef PIVector2D _V2D; typedef PIMathMatrix _CMatrix; typedef PIMathVector _CMCol; - typedef PIMathVector _CMRow; public: PIMathMatrix(const uint cols = 3, const uint rows = 3) {resize(cols, rows);} - PIMathMatrix(const uint cols, const uint rows, const PIVector & val) {resize(cols, rows); int i=0; PIMM_FOR_I(c, r) (*this)[r][c] = val[i++];} - PIMathMatrix(const PIVector > & 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 & 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 & val) {resize(cols, rows); int i=0; PIMM_FOR_I(c, r) _V2D::element(r, c) = val[i++];} + PIMathMatrix(const PIVector > & 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 & 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 & val) {return _CMatrix(val.size(), 1, val.toVector());} static _CMatrix matrixCol(const PIMathVector & 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 > & v) {*this = PIVector2D(v); return *this;} + _CMatrix & operator =(const PIVector > & 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 -inline std::ostream & operator <<(std::ostream & s, const PIMathMatrix & 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 & 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 -inline PICout operator <<(PICout s, const PIMathMatrix & 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 & 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 @@ -528,8 +507,8 @@ inline PIMathMatrix operator *(const PIMathMatrix & 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 operator *(const PIMathMatrix & 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 PIMathMatrixd; template PIMathMatrix > hermitian(const PIMathMatrix > & m) { PIMathMatrix > 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(); } diff --git a/src_main/math/pimathsolver.cpp b/src_main/math/pimathsolver.cpp index c1e0143c..721ea7be 100644 --- a/src_main/math/pimathsolver.cpp +++ b/src_main/math/pimathsolver.cpp @@ -98,7 +98,7 @@ void PIMathSolver::fromTF(const TransferFunction & TF) { a1 /= a0; b1 /= a0; - d[0] = b1[0]; // ������������ ������ d + d[0] = b1[0]; for (uint i = 1; i < size + 1; ++i) { sum = 0.; for (uint m = 0; m < i; ++m) @@ -106,7 +106,7 @@ void PIMathSolver::fromTF(const TransferFunction & TF) { d[i] = b1[i] - sum; } - for (uint i = 0; i < size - 1; ++i) // ��������� ������� � + for (uint i = 0; i < size - 1; ++i) for (uint j = 0; j < size; ++j) A[j][i] = (j == i + 1); for (uint i = 0; i < size; ++i)