From 18ba239ff3cb97dae4cf807aa386e53b75404012 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Wed, 17 Dec 2025 15:39:34 +0000 Subject: [PATCH 01/16] Refactor Matrix to have a RectArray but not be a RectArray. This removes the inheritance of Matrix from RectArray, alongside changing the internal representation of RectArray to be a continuous std::vector of elements rather than a T ** with dodgy pointer arithmetic for index range adjustments. --- src/core/matrix.h | 140 +++++++++++++++++----- src/core/matrix.imp | 281 ++++++++++++++++++-------------------------- src/core/recarray.h | 173 ++++++++++----------------- 3 files changed, 289 insertions(+), 305 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 04920da69f..5d96ac7426 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -36,49 +36,77 @@ class SingularMatrixException final : public std::runtime_error { template Vector operator*(const Vector &, const Matrix &); -template class Matrix : public RectArray { - friend Vector operator* <>(const Vector &, const Matrix &); +template class Matrix { + friend Vector operator* <>(const Vector &, const Matrix &); + + RectArray m_data; public: /// @name Lifecycle //@{ - Matrix(); - Matrix(unsigned int rows, unsigned int cols); - Matrix(unsigned int rows, unsigned int cols, int minrows); - Matrix(int rl, int rh, int cl, int ch); - Matrix(const Matrix &); - ~Matrix() override; - - Matrix &operator=(const Matrix &); - Matrix &operator=(const T &); + Matrix() = default; + Matrix(unsigned int rows, unsigned int cols) : m_data(rows, cols) {} + Matrix(unsigned int rows, unsigned int cols, int minrows) + : m_data(minrows, minrows + rows - 1, 1, cols) + { + } + Matrix(int rl, int rh, int cl, int ch) : m_data(rl, rh, cl, ch) {} + Matrix(const Matrix &) = default; + Matrix(Matrix &&) noexcept = default; + ~Matrix() = default; + + Matrix &operator=(const Matrix &) = default; + Matrix &operator=(Matrix &&) noexcept = default; + Matrix &operator=(const T &); //@} + // element access + T &operator()(int r, int c) { return m_data(r, c); } + const T &operator()(int r, int c) const { return m_data(r, c); } + + // bounds / dimensions + int MinRow() const { return m_data.MinRow(); } + int MaxRow() const { return m_data.MaxRow(); } + int MinCol() const { return m_data.MinCol(); } + int MaxCol() const { return m_data.MaxCol(); } + int NumRows() const { return m_data.NumRows(); } + int NumColumns() const { return m_data.NumColumns(); } + + bool IsSquare() const { return MinRow() == MinCol() && MaxRow() == MaxCol(); } + + // row ops used internally + void SwitchRows(int i, int j) { m_data.SwitchRows(i, j); } + template void GetColumn(int j, V &) const; + template void SetColumn(int j, const V &); + template void GetRow(int row, V &) const; + template void SetRow(int row, const V &); + // vector helpers used internally + template bool CheckRow(const V &v) const { return m_data.CheckRow(v); } + template bool CheckColumn(const V &v) const { return m_data.CheckColumn(v); } + bool CheckBounds(const Matrix &M) const { return m_data.CheckBounds(M.m_data); } + /// @name Extracting rows and columns //@{ - bool IsSquare() const - { - return this->MinRow() == this->MinCol() && this->MaxRow() == this->MaxCol(); - } Vector Row(int) const; Vector Column(int) const; //@} /// @name Comparison operators //@{ - bool operator==(const Matrix &) const; - bool operator!=(const Matrix &) const; + bool operator==(const Matrix &) const; + bool operator!=(const Matrix &) const; bool operator==(const T &) const; bool operator!=(const T &) const; //@} /// @name Additive operators //@{ - Matrix operator+(const Matrix &) const; - Matrix operator-(const Matrix &) const; - Matrix &operator+=(const Matrix &); - Matrix &operator-=(const Matrix &); + Matrix operator+(const Matrix &) const; + Matrix operator-(const Matrix &) const; + Matrix &operator+=(const Matrix &); + Matrix &operator-=(const Matrix &); - Matrix operator-(); + Matrix operator-(); //@} /// @name Multiplicative operators @@ -87,18 +115,18 @@ template class Matrix : public RectArray { void CMultiply(const Vector &, Vector &) const; /// "in-place" row (transposed) multiply void RMultiply(const Vector &, Vector &) const; - Matrix operator*(const Matrix &) const; + Matrix operator*(const Matrix &) const; Vector operator*(const Vector &) const; - Matrix operator*(const T &) const; - Matrix &operator*=(const T &); + Matrix operator*(const T &) const; + Matrix &operator*=(const T &); - Matrix operator/(const T &) const; - Matrix &operator/=(const T &); + Matrix operator/(const T &) const; + Matrix &operator/=(const T &); //@ /// @name Other operations //@{ - Matrix Transpose() const; + Matrix Transpose() const; /// Set matrix to identity matrix void MakeIdent(); void Pivot(int, int); @@ -110,6 +138,62 @@ template class Matrix : public RectArray { template Vector operator*(const Vector &, const Matrix &); +template template void Matrix::GetColumn(int col, Vector &v) const +{ + if (col < MinCol() || col > MaxCol()) { + throw std::out_of_range("Index out of range in Matrix::GetColumn"); + } + if (!CheckColumn(v)) { + throw DimensionException(); + } + + for (int i = MinRow(); i <= MaxRow(); ++i) { + v[i] = (*this)(i, col); + } +} + +template template void Matrix::SetColumn(int col, const Vector &v) +{ + if (col < MinCol() || col > MaxCol()) { + throw std::out_of_range("Index out of range in Matrix::SetColumn"); + } + if (!CheckColumn(v)) { + throw DimensionException(); + } + + for (int i = MinRow(); i <= MaxRow(); ++i) { + (*this)(i, col) = v[i]; + } +} + +template template void Matrix::GetRow(int row, Vector &v) const +{ + if (row < MinRow() || row > MaxRow()) { + throw std::out_of_range("Index out of range in Matrix::GetRow"); + } + if (!CheckRow(v)) { + throw DimensionException(); + } + + for (int j = MinCol(); j <= MaxCol(); ++j) { + v[j] = (*this)(row, j); + } +} + +template template void Matrix::SetRow(int row, const Vector &v) +{ + if (row < MinRow() || row > MaxRow()) { + throw std::out_of_range("Index out of range in Matrix::SetRow"); + } + if (!CheckRow(v)) { + throw DimensionException(); + } + + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(row, j) = v[j]; + } +} + } // end namespace Gambit #endif // GAMBIT_CORE_MATRIX_H diff --git a/src/core/matrix.imp b/src/core/matrix.imp index 5650583fd2..0fda5cbca8 100644 --- a/src/core/matrix.imp +++ b/src/core/matrix.imp @@ -28,37 +28,11 @@ namespace Gambit { // Matrix: Constructors, destructors, constructive operators //------------------------------------------------------------------------- -template Matrix::Matrix() = default; - -template -Matrix::Matrix(unsigned int rows, unsigned int cols) : RectArray(rows, cols) -{ -} - -template -Matrix::Matrix(unsigned int rows, unsigned int cols, int minrows) - : RectArray(minrows, minrows + rows - 1, 1, cols) -{ -} - -template Matrix::Matrix(int rl, int rh, int cl, int ch) : RectArray(rl, rh, cl, ch) -{ -} - -template Matrix::Matrix(const Matrix &M) : RectArray(M) {} - -template Matrix::~Matrix() = default; - -template Matrix &Matrix::operator=(const Matrix &M) -{ - RectArray::operator=(M); - return *this; -} template Matrix &Matrix::operator=(const T &c) { - for (int i = this->minrow; i <= this->maxrow; i++) { - for (int j = this->mincol; j <= this->maxcol; j++) { + for (int i = MinRow(); i <= MaxRow(); i++) { + for (int j = MinCol(); j <= MaxCol(); j++) { (*this)(i, j) = c; } } @@ -67,9 +41,9 @@ template Matrix &Matrix::operator=(const T &c) template Matrix Matrix::operator-() { - Matrix tmp(this->minrow, this->maxrow, this->mincol, this->maxcol); - for (int i = this->minrow; i <= this->maxrow; i++) { - for (int j = this->mincol; j <= this->maxcol; j++) { + Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); + for (int i = MinRow(); i <= MaxRow(); i++) { + for (int j = MinCol(); j <= MaxCol(); j++) { tmp(i, j) = -(*this)(i, j); } } @@ -86,17 +60,14 @@ template Matrix Matrix::operator+(const Matrix &M) const throw DimensionException(); } - const Matrix tmp(this->minrow, this->maxrow, this->mincol, this->maxcol); - for (int i = this->minrow; i <= this->maxrow; i++) { - T *src1 = this->data[i] + this->mincol; - T *src2 = M.data[i] + this->mincol; - T *dst = tmp.data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) = *(src1++) + *(src2++); + Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); + + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + tmp(i, j) = (*this)(i, j) + M(i, j); } - // assert((dst - 1) == tmp.data[i] + this->maxcol ); } + return tmp; } @@ -106,17 +77,14 @@ template Matrix Matrix::operator-(const Matrix &M) const throw DimensionException(); } - const Matrix tmp(this->minrow, this->maxrow, this->mincol, this->maxcol); - for (int i = this->minrow; i <= this->maxrow; i++) { - T *src1 = this->data[i] + this->mincol; - T *src2 = M.data[i] + this->mincol; - T *dst = tmp.data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) = *(src1++) - *(src2++); + Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); + + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + tmp(i, j) = (*this)(i, j) - M(i, j); } - // assert((dst - 1) == tmp.data[i] + this->maxcol); } + return tmp; } @@ -126,15 +94,13 @@ template Matrix &Matrix::operator+=(const Matrix &M) throw DimensionException(); } - for (int i = this->minrow; i <= this->maxrow; i++) { - T *src = M.data[i] + this->mincol; - T *dst = this->data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) += *(src++); + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(i, j) += M(i, j); } } - return (*this); + + return *this; } template Matrix &Matrix::operator-=(const Matrix &M) @@ -143,16 +109,13 @@ template Matrix &Matrix::operator-=(const Matrix &M) throw DimensionException(); } - for (int i = this->minrow; i <= this->maxrow; i++) { - T *src = M.data[i] + this->mincol; - T *dst = this->data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) -= *(src++); + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(i, j) -= M(i, j); } - // assert((dst - 1) == this->data[i] + this->maxcol); } - return (*this); + + return *this; } //------------------------------------------------------------------------- @@ -165,29 +128,28 @@ template void Matrix::CMultiply(const Vector &in, Vector &out throw DimensionException(); } - for (int i = this->minrow; i <= this->maxrow; i++) { + for (int i = MinRow(); i <= MaxRow(); ++i) { T sum = (T)0; - T *src1 = this->data[i] + this->mincol; auto src2 = in.begin(); - int j = this->maxcol - this->mincol + 1; - while (j--) { - sum += *(src1++) * *(src2++); + for (int j = MinCol(); j <= MaxCol(); ++j) { + sum += (*this)(i, j) * *(src2++); } + out[i] = sum; } } template Matrix Matrix::operator*(const Matrix &M) const { - if (this->mincol != M.minrow || this->maxcol != M.maxrow) { + if (MinCol() != M.MinRow() || MaxCol() != M.MaxRow()) { throw DimensionException(); } - Matrix tmp(this->minrow, this->maxrow, M.mincol, M.maxcol); - Vector column(M.minrow, M.maxrow); - Vector result(this->minrow, this->maxrow); - for (int j = M.mincol; j <= M.maxcol; j++) { + Matrix tmp(MinRow(), MaxRow(), M.MinCol(), M.MaxCol()); + Vector column(M.MinRow(), M.MaxRow()); + Vector result(MinRow(), MaxRow()); + for (int j = M.MinCol(); j <= M.MaxCol(); j++) { M.GetColumn(j, column); CMultiply(column, result); tmp.SetColumn(j, result); @@ -201,7 +163,7 @@ template Vector Matrix::operator*(const Vector &v) const throw DimensionException(); } - Vector tmp(this->minrow, this->maxrow); + Vector tmp(MinRow(), MaxRow()); CMultiply(v, tmp); return tmp; } @@ -213,16 +175,12 @@ template void Matrix::RMultiply(const Vector &in, Vector &out } out = (T)0; - for (int i = this->minrow; i <= this->maxrow; i++) { - T k = in[i]; - T *src = this->data[i] + this->mincol; - auto dst = out.begin(); - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) += *(src++) * k; + for (int i = MinRow(); i <= MaxRow(); ++i) { + T k = in[i]; + for (int j = MinCol(); j <= MaxCol(); ++j) { + out[j] += (*this)(i, j) * k; } - // assert(src - 1 == this->data[i] + this->maxcol); } } @@ -240,29 +198,25 @@ template Vector operator*(const Vector &v, const Matrix &M) template Matrix Matrix::operator*(const T &s) const { - const Matrix tmp(this->minrow, this->maxrow, this->mincol, this->maxcol); - for (int i = this->minrow; i <= this->maxrow; i++) { - T *src = this->data[i] + this->mincol; - T *dst = tmp.data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) = *(src++) * s; + Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); + + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + tmp(i, j) = (*this)(i, j) * s; } - // assert((src - 1) == this->data[i] + this->maxcol); } + return tmp; } template Matrix &Matrix::operator*=(const T &s) { - for (int i = this->minrow; i <= this->maxrow; i++) { - T *dst = this->data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) *= s; + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(i, j) *= s; } } - return (*this); + return *this; } template Matrix Matrix::operator/(const T &s) const @@ -271,15 +225,14 @@ template Matrix Matrix::operator/(const T &s) const throw ZeroDivideException(); } - const Matrix tmp(this->minrow, this->maxrow, this->mincol, this->maxcol); - for (int i = this->minrow; i <= this->maxrow; i++) { - T *src = this->data[i] + this->mincol; - T *dst = tmp.data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) = *(src++) / s; + Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); + + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + tmp(i, j) = (*this)(i, j) / s; } } + return tmp; } @@ -289,14 +242,13 @@ template Matrix &Matrix::operator/=(const T &s) throw ZeroDivideException(); } - for (int i = this->minrow; i <= this->maxrow; i++) { - T *dst = this->data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) /= s; + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(i, j) /= s; } } - return (*this); + + return *this; } //------------------------------------------------------------------------- @@ -305,10 +257,10 @@ template Matrix &Matrix::operator/=(const T &s) template Matrix Matrix::Transpose() const { - Matrix tmp(this->mincol, this->maxcol, this->minrow, this->maxrow); + Matrix tmp(MinCol(), MaxCol(), MinRow(), MaxRow()); - for (int i = this->minrow; i <= this->maxrow; i++) { - for (int j = this->mincol; j <= this->maxcol; j++) { + for (int i = MinRow(); i <= MaxRow(); i++) { + for (int j = MinCol(); j <= MaxCol(); j++) { tmp(j, i) = (*this)(i, j); } } @@ -326,18 +278,14 @@ template bool Matrix::operator==(const Matrix &M) const throw DimensionException(); } - for (int i = this->minrow; i <= this->maxrow; i++) { - // inner loop - T *src1 = M.data[i] + this->mincol; - T *src2 = this->data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - if (*(src1++) != *(src2++)) { + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + if ((*this)(i, j) != M(i, j)) { return false; } } - // assert(src1 - 1 == M.data[i] + this->maxcol); } + return true; } @@ -345,27 +293,25 @@ template bool Matrix::operator!=(const Matrix &M) const { return template bool Matrix::operator==(const T &s) const { - for (int i = this->minrow; i <= this->maxrow; i++) { - T *src = this->data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - if (*(src++) != s) { + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { + if ((*this)(i, j) != s) { return false; } } - // assert(src - 1 == this->data[i] + this->maxcol); } return true; } + template bool Matrix::operator!=(const T &s) const { return !(*this == s); } // Information template Vector Matrix::Row(int i) const { - Vector answer(this->mincol, this->maxcol); - for (int j = this->mincol; j <= this->maxcol; j++) { + Vector answer(MinCol(), MaxCol()); + for (int j = MinCol(); j <= MaxCol(); j++) { answer[j] = (*this)(i, j); } return answer; @@ -373,8 +319,8 @@ template Vector Matrix::Row(int i) const template Vector Matrix::Column(int j) const { - Vector answer(this->minrow, this->maxrow); - for (int i = this->minrow; i <= this->maxrow; i++) { + Vector answer(MinRow(), MaxRow()); + for (int i = MinRow(); i <= MaxRow(); i++) { answer[i] = (*this)(i, j); } return answer; @@ -387,8 +333,8 @@ template void Matrix::MakeIdent() if (!IsSquare()) { throw DimensionException(); } - for (int i = this->minrow; i <= this->maxrow; i++) { - for (int j = this->mincol; j <= this->maxcol; j++) { + for (int i = MinRow(); i <= MaxRow(); i++) { + for (int j = MinCol(); j <= MaxCol(); j++) { if (i == j) { (*this)(i, j) = (T)1; } @@ -404,27 +350,24 @@ template void Matrix::Pivot(int row, int col) if (!this->CheckRow(row) || !this->CheckColumn(col)) { throw std::out_of_range("Index out of range in Matrix::Pivot"); } - if (this->data[row][col] == (T)0) { + if ((*this)(row, col) == (T)0) { throw ZeroDivideException(); } - T mult = (T)1 / this->data[row][col]; - for (int j = this->mincol; j <= this->maxcol; j++) { - this->data[row][j] *= mult; + T inv = (T)1 / (*this)(row, col); + + // scale pivot row + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(row, j) *= inv; } - for (int i = this->minrow; i <= this->maxrow; i++) { + + // eliminate column + for (int i = MinRow(); i <= MaxRow(); ++i) { if (i != row) { - mult = this->data[i][col]; - - // inner loop - T *src = this->data[row] + this->mincol; - T *dst = this->data[i] + this->mincol; - int j = this->maxcol - this->mincol + 1; - while (j--) { - *(dst++) -= *(src++) * mult; + T mult = (*this)(i, col); + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(i, j) -= (*this)(row, j) * mult; } - // assert( dst-1 == this->data[i] + this->maxcol ); // debug - // end inner loop } } } @@ -436,12 +379,12 @@ template Matrix Matrix::Inverse() const } Matrix copy(*this); - Matrix inv(this->MaxRow(), this->MaxRow()); + Matrix inv(MaxRow(), MaxRow()); // initialize inverse matrix and prescale row vectors - for (int i = this->MinRow(); i <= this->MaxRow(); i++) { + for (int i = MinRow(); i <= MaxRow(); i++) { T max = (T)0; - for (int j = this->MinCol(); j <= this->MaxCol(); j++) { + for (int j = MinCol(); j <= MaxCol(); j++) { T abs = copy(i, j); if (abs < (T)0) { abs = -abs; @@ -456,7 +399,7 @@ template Matrix Matrix::Inverse() const } T scale = (T)1 / max; - for (int j = this->MinCol(); j <= this->MaxCol(); j++) { + for (int j = MinCol(); j <= MaxCol(); j++) { copy(i, j) *= scale; if (i == j) { inv(i, j) = scale; @@ -467,14 +410,14 @@ template Matrix Matrix::Inverse() const } } - for (int i = this->MinCol(); i <= this->MaxCol(); i++) { + for (int i = MinCol(); i <= MaxCol(); i++) { // find pivot row T max = copy(i, i); if (max < (T)0) { max = -max; } int row = i; - for (int j = i + 1; j <= this->MaxRow(); j++) { + for (int j = i + 1; j <= MaxRow(); j++) { T abs = copy(j, i); if (abs < (T)0) { abs = -abs; @@ -493,16 +436,16 @@ template Matrix Matrix::Inverse() const inv.SwitchRows(i, row); // scale pivot row T factor = (T)1 / copy(i, i); - for (int k = this->MinCol(); k <= this->MaxCol(); k++) { + for (int k = MinCol(); k <= MaxCol(); k++) { copy(i, k) *= factor; inv(i, k) *= factor; } // reduce other rows - for (int j = this->MinRow(); j <= this->MaxRow(); j++) { + for (int j = MinRow(); j <= MaxRow(); j++) { if (j != i) { T mult = copy(j, i); - for (int k = this->MinCol(); k <= this->MaxCol(); k++) { + for (int k = MinCol(); k <= MaxCol(); k++) { copy(j, k) -= copy(i, k) * mult; inv(j, k) -= inv(i, k) * mult; } @@ -522,44 +465,44 @@ template T Matrix::Determinant() const T factor = (T)1; Matrix M(*this); - for (int row = this->MinRow(); row <= this->MaxRow(); row++) { + for (int row = MinRow(); row <= MaxRow(); row++) { // Experience (as of 3/22/99) suggests that, in the interest of // numerical stability, it might be best to do Gaussian // elimination with respect to the row (of those feasible) // whose entry has the largest absolute value. int swap_row = row; - for (int i = row + 1; i <= this->MaxRow(); i++) { - if (abs(M.data[i][row]) > abs(M.data[swap_row][row])) { + for (int i = row + 1; i <= MaxRow(); i++) { + if (abs(M(i, row)) > abs(M(swap_row, row))) { swap_row = i; } } if (swap_row != row) { M.SwitchRows(row, swap_row); - for (int j = this->MinCol(); j <= this->MaxCol(); j++) { - M.data[row][j] *= (T)-1; + for (int j = MinCol(); j <= MaxCol(); j++) { + M(row, j) *= (T)-1; } } - if (M.data[row][row] == (T)0) { + if (M(row, row) == (T)0) { return (T)0; } // now do row operations to clear the row'th column // below the diagonal - for (int row1 = row + 1; row1 <= this->MaxRow(); row1++) { - factor = -M.data[row1][row] / M.data[row][row]; - for (int i = this->MinCol(); i <= this->MaxCol(); i++) { - M.data[row1][i] += M.data[row][i] * factor; + for (int row1 = row + 1; row1 <= MaxRow(); row1++) { + factor = -M(row1, row) / M(row, row); + for (int i = MinCol(); i <= MaxCol(); i++) { + M(row1, i) += M(row, i) * factor; } } } // finally we multiply the diagonal elements T det = (T)1; - for (int row = this->MinRow(); row <= this->MaxRow(); row++) { - det *= M.data[row][row]; + for (int row = MinRow(); row <= MaxRow(); row++) { + det *= M(row, row); } return det; } diff --git a/src/core/recarray.h b/src/core/recarray.h index 90540087ba..7d239d5173 100644 --- a/src/core/recarray.h +++ b/src/core/recarray.h @@ -31,8 +31,39 @@ namespace Gambit { template class RectArray { protected: int minrow, maxrow, mincol, maxcol; - T **data; + std::vector storage; + size_t row_stride() const { return static_cast(maxcol - mincol + 1); } + + size_t index(int r, int c) const + { + return static_cast(r - minrow) * row_stride() + static_cast(c - mincol); + } + +public: + /// @name Lifecycle + //@{ + RectArray() : minrow(1), maxrow(0), mincol(1), maxcol(0), storage() {} + RectArray(unsigned int nrows, unsigned int ncols); + RectArray(int minr, int maxr, int minc, int maxc) + : minrow(minr), maxrow(maxr), mincol(minc), maxcol(maxc), + storage((maxr >= minr && maxc >= minc) ? (maxr - minr + 1) * (maxc - minc + 1) : 0) + { + } + RectArray(const RectArray &) = default; + RectArray(RectArray &&) noexcept = default; + + virtual ~RectArray() = default; + + RectArray &operator=(const RectArray &) = default; + RectArray &operator=(RectArray &&) noexcept = default; + //@} + + /// check array for same row and column boundaries + bool CheckBounds(const RectArray &m) const + { + return (minrow == m.minrow && maxrow == m.maxrow && mincol == m.mincol && maxcol == m.maxcol); + } /// @name Range checking functions; returns true only if valid index/size //@{ /// check for correct row index @@ -51,25 +82,8 @@ template class RectArray { } /// check row and column indices bool Check(int row, int col) const { return CheckRow(row) && CheckColumn(col); } - /// check matrix for same row and column boundaries - bool CheckBounds(const RectArray &m) const - { - return (minrow == m.minrow && maxrow == m.maxrow && mincol == m.mincol && maxcol == m.maxcol); - } - //@} - -public: - /// @name Lifecycle - //@{ - RectArray() : minrow(1), maxrow(0), mincol(1), maxcol(0), data(nullptr) {} - RectArray(unsigned int nrows, unsigned int ncols); - RectArray(int minr, int maxr, int minc, int maxc); - RectArray(const RectArray &); - virtual ~RectArray(); - RectArray &operator=(const RectArray &); //@} - /// @name General data access //@{ size_t NumRows() const { return maxrow - minrow + 1; } @@ -87,14 +101,14 @@ template class RectArray { if (!Check(r, c)) { throw std::out_of_range("Index out of range in RectArray"); } - return data[r][c]; + return storage[index(r, c)]; } const T &operator()(int r, int c) const { if (!Check(r, c)) { throw std::out_of_range("Index out of range in RectArray"); } - return data[r][c]; + return storage[index(r, c)]; } //@} @@ -111,7 +125,17 @@ template class RectArray { if (!Check(i, j)) { throw std::out_of_range("Index out of range in RectArray"); } - std::swap(data[i], data[j]); + if (i == j) { + return; + } + + auto stride = row_stride(); + size_t ai = index(i, mincol); + size_t aj = index(j, mincol); + + for (size_t k = 0; k < stride; ++k) { + std::swap(storage[ai + k], storage[aj + k]); + } } template void GetRow(int, Vector &) const; @@ -125,77 +149,8 @@ template class RectArray { //------------------------------------------------------------------------ template -RectArray::RectArray(unsigned int rows, unsigned int cols) - : minrow(1), maxrow(rows), mincol(1), maxcol(cols), - data((rows > 0) ? new T *[maxrow] - 1 : nullptr) +RectArray::RectArray(unsigned int rows, unsigned int cols) : RectArray(1, rows, 1, cols) { - for (int i = 1; i <= maxrow; data[i++] = (cols > 0) ? new T[maxcol] - 1 : nullptr) - ; -} - -template -RectArray::RectArray(int minr, int maxr, int minc, int maxc) - : minrow(minr), maxrow(maxr), mincol(minc), maxcol(maxc), - data((maxrow >= minrow) ? new T *[maxrow - minrow + 1] - minrow : nullptr) -{ - for (int i = minrow; i <= maxrow; - data[i++] = (maxcol - mincol + 1) ? new T[maxcol - mincol + 1] - mincol : nullptr) - ; -} - -template -RectArray::RectArray(const RectArray &a) - : minrow(a.minrow), maxrow(a.maxrow), mincol(a.mincol), maxcol(a.maxcol), - data((maxrow >= minrow) ? new T *[maxrow - minrow + 1] - minrow : nullptr) -{ - for (int i = minrow; i <= maxrow; i++) { - data[i] = (maxcol >= mincol) ? new T[maxcol - mincol + 1] - mincol : nullptr; - for (int j = mincol; j <= maxcol; j++) { - data[i][j] = a.data[i][j]; - } - } -} - -template RectArray::~RectArray() -{ - for (int i = minrow; i <= maxrow; i++) { - if (data[i]) { - delete[] (data[i] + mincol); - } - } - if (data) { - delete[] (data + minrow); - } -} - -template RectArray &RectArray::operator=(const RectArray &a) -{ - if (this != &a) { - for (int i = minrow; i <= maxrow; i++) { - if (data[i]) { - delete[] (data[i] + mincol); - } - } - if (data) { - delete[] (data + minrow); - } - - minrow = a.minrow; - maxrow = a.maxrow; - mincol = a.mincol; - maxcol = a.maxcol; - - data = (maxrow >= minrow) ? new T *[maxrow - minrow + 1] - minrow : nullptr; - - for (int i = minrow; i <= maxrow; i++) { - data[i] = (maxcol >= mincol) ? new T[maxcol - mincol + 1] - mincol : nullptr; - for (int j = mincol; j <= maxcol; j++) { - data[i][j] = a.data[i][j]; - } - } - } - - return *this; } //------------------------------------------------------------------------ @@ -208,11 +163,12 @@ template void RectArray::RotateUp(int lo, int hi) throw std::out_of_range("Index out of range in RectArray"); } - T *temp = data[lo]; - for (int k = lo; k < hi; k++) { - data[k] = data[k + 1]; - } - data[hi] = temp; + auto stride = row_stride(); + + size_t first = index(lo, mincol); + size_t last = index(hi + 1, mincol); + + std::rotate(storage.begin() + first, storage.begin() + first + stride, storage.begin() + last); } template void RectArray::RotateDown(int lo, int hi) @@ -221,11 +177,12 @@ template void RectArray::RotateDown(int lo, int hi) throw std::out_of_range("Index out of range in RectArray"); } - T *temp = data[hi]; - for (int k = hi; k > lo; k--) { - data[k] = data[k - 1]; - } - data[lo] = temp; + auto stride = row_stride(); + + size_t first = index(lo, mincol); + size_t last = index(hi + 1, mincol); + + std::rotate(storage.begin() + first, storage.begin() + last - stride, storage.begin() + last); } //------------------------------------------------------------------------- @@ -240,9 +197,9 @@ template template void RectArray::GetRow(int row, Ve if (!CheckRow(v)) { throw DimensionException(); } - T *rowptr = data[row]; - for (int i = mincol; i <= maxcol; i++) { - v[i] = rowptr[i]; + size_t base = index(row, mincol); + for (int c = mincol; c <= maxcol; ++c) { + v[c] = storage[base + (c - mincol)]; } } @@ -258,8 +215,8 @@ template template void RectArray::GetColumn(int col, if (!CheckColumn(v)) { throw DimensionException(); } - for (int i = minrow; i <= maxrow; i++) { - v[i] = data[i][col]; + for (int r = minrow; r <= maxrow; ++r) { + v[r] = storage[index(r, col)]; } } @@ -271,8 +228,8 @@ template template void RectArray::SetColumn(int col, if (!CheckColumn(v)) { throw DimensionException(); } - for (int i = minrow; i <= maxrow; i++) { - data[i][col] = v[i]; + for (int r = minrow; r <= maxrow; ++r) { + storage[index(r, col)] = v[r]; } } From c7b6a7f68471c6c392817a96fbfd9371509e613a Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Wed, 17 Dec 2025 21:12:11 +0000 Subject: [PATCH 02/16] Add missing include --- src/core/recarray.h | 80 ++++++++++++++++++++++++--------------------- 1 file changed, 42 insertions(+), 38 deletions(-) diff --git a/src/core/recarray.h b/src/core/recarray.h index 7d239d5173..bc0e50bdc4 100644 --- a/src/core/recarray.h +++ b/src/core/recarray.h @@ -23,6 +23,7 @@ #ifndef GAMBIT_CORE_RECARRAY_H #define GAMBIT_CORE_RECARRAY_H +#include #include "util.h" namespace Gambit { @@ -30,24 +31,24 @@ namespace Gambit { /// This class implements a rectangular (two-dimensional) array template class RectArray { protected: - int minrow, maxrow, mincol, maxcol; - std::vector storage; + int m_minrow, m_maxrow, m_mincol, m_maxcol; + std::vector m_storage; - size_t row_stride() const { return static_cast(maxcol - mincol + 1); } + size_t row_stride() const { return static_cast(m_maxcol - m_mincol + 1); } size_t index(int r, int c) const { - return static_cast(r - minrow) * row_stride() + static_cast(c - mincol); + return static_cast(r - m_minrow) * row_stride() + static_cast(c - m_mincol); } public: /// @name Lifecycle //@{ - RectArray() : minrow(1), maxrow(0), mincol(1), maxcol(0), storage() {} + RectArray() : m_minrow(1), m_maxrow(0), m_mincol(1), m_maxcol(0), m_storage() {} RectArray(unsigned int nrows, unsigned int ncols); RectArray(int minr, int maxr, int minc, int maxc) - : minrow(minr), maxrow(maxr), mincol(minc), maxcol(maxc), - storage((maxr >= minr && maxc >= minc) ? (maxr - minr + 1) * (maxc - minc + 1) : 0) + : m_minrow(minr), m_maxrow(maxr), m_mincol(minc), m_maxcol(maxc), + m_storage((maxr >= minr && maxc >= minc) ? (maxr - minr + 1) * (maxc - minc + 1) : 0) { } RectArray(const RectArray &) = default; @@ -62,23 +63,24 @@ template class RectArray { /// check array for same row and column boundaries bool CheckBounds(const RectArray &m) const { - return (minrow == m.minrow && maxrow == m.maxrow && mincol == m.mincol && maxcol == m.maxcol); + return (m_minrow == m.m_minrow && m_maxrow == m.m_maxrow && m_mincol == m.m_mincol && + m_maxcol == m.m_maxcol); } /// @name Range checking functions; returns true only if valid index/size //@{ /// check for correct row index - bool CheckRow(int row) const { return (minrow <= row && row <= maxrow); } + bool CheckRow(int row) const { return (m_minrow <= row && row <= m_maxrow); } /// check row vector for correct column boundaries template bool CheckRow(const Vector &v) const { - return (v.front_index() == mincol && v.back_index() == maxcol); + return (v.front_index() == m_mincol && v.back_index() == m_maxcol); } /// check for correct column index - bool CheckColumn(int col) const { return (mincol <= col && col <= maxcol); } + bool CheckColumn(int col) const { return (m_mincol <= col && col <= m_maxcol); } /// check column vector for correct row boundaries template bool CheckColumn(const Vector &v) const { - return (v.front_index() == minrow && v.back_index() == maxrow); + return (v.front_index() == m_minrow && v.back_index() == m_maxrow); } /// check row and column indices bool Check(int row, int col) const { return CheckRow(row) && CheckColumn(col); } @@ -86,12 +88,12 @@ template class RectArray { //@} /// @name General data access //@{ - size_t NumRows() const { return maxrow - minrow + 1; } - size_t NumColumns() const { return maxcol - mincol + 1; } - int MinRow() const { return minrow; } - int MaxRow() const { return maxrow; } - int MinCol() const { return mincol; } - int MaxCol() const { return maxcol; } + size_t NumRows() const { return m_maxrow - m_minrow + 1; } + size_t NumColumns() const { return m_maxcol - m_mincol + 1; } + int MinRow() const { return m_minrow; } + int MaxRow() const { return m_maxrow; } + int MinCol() const { return m_mincol; } + int MaxCol() const { return m_maxcol; } //@} /// @name Indexing operations @@ -101,14 +103,14 @@ template class RectArray { if (!Check(r, c)) { throw std::out_of_range("Index out of range in RectArray"); } - return storage[index(r, c)]; + return m_storage[index(r, c)]; } const T &operator()(int r, int c) const { if (!Check(r, c)) { throw std::out_of_range("Index out of range in RectArray"); } - return storage[index(r, c)]; + return m_storage[index(r, c)]; } //@} @@ -130,11 +132,11 @@ template class RectArray { } auto stride = row_stride(); - size_t ai = index(i, mincol); - size_t aj = index(j, mincol); + size_t ai = index(i, m_mincol); + size_t aj = index(j, m_mincol); for (size_t k = 0; k < stride; ++k) { - std::swap(storage[ai + k], storage[aj + k]); + std::swap(m_storage[ai + k], m_storage[aj + k]); } } template void GetRow(int, Vector &) const; @@ -159,30 +161,32 @@ RectArray::RectArray(unsigned int rows, unsigned int cols) : RectArray(1, row template void RectArray::RotateUp(int lo, int hi) { - if (lo < minrow || hi < lo || maxrow < hi) { + if (lo < m_minrow || hi < lo || m_maxrow < hi) { throw std::out_of_range("Index out of range in RectArray"); } auto stride = row_stride(); - size_t first = index(lo, mincol); - size_t last = index(hi + 1, mincol); + size_t first = index(lo, m_mincol); + size_t last = index(hi + 1, m_mincol); - std::rotate(storage.begin() + first, storage.begin() + first + stride, storage.begin() + last); + std::rotate(m_storage.begin() + first, m_storage.begin() + first + stride, + m_storage.begin() + last); } template void RectArray::RotateDown(int lo, int hi) { - if (lo < minrow || hi < lo || maxrow < hi) { + if (lo < m_minrow || hi < lo || m_maxrow < hi) { throw std::out_of_range("Index out of range in RectArray"); } auto stride = row_stride(); - size_t first = index(lo, mincol); - size_t last = index(hi + 1, mincol); + size_t first = index(lo, m_mincol); + size_t last = index(hi + 1, m_mincol); - std::rotate(storage.begin() + first, storage.begin() + last - stride, storage.begin() + last); + std::rotate(m_storage.begin() + first, m_storage.begin() + last - stride, + m_storage.begin() + last); } //------------------------------------------------------------------------- @@ -197,9 +201,9 @@ template template void RectArray::GetRow(int row, Ve if (!CheckRow(v)) { throw DimensionException(); } - size_t base = index(row, mincol); - for (int c = mincol; c <= maxcol; ++c) { - v[c] = storage[base + (c - mincol)]; + size_t base = index(row, m_mincol); + for (int c = m_mincol; c <= m_maxcol; ++c) { + v[c] = m_storage[base + (c - m_mincol)]; } } @@ -215,8 +219,8 @@ template template void RectArray::GetColumn(int col, if (!CheckColumn(v)) { throw DimensionException(); } - for (int r = minrow; r <= maxrow; ++r) { - v[r] = storage[index(r, col)]; + for (int r = m_minrow; r <= m_maxrow; ++r) { + v[r] = m_storage[index(r, col)]; } } @@ -228,8 +232,8 @@ template template void RectArray::SetColumn(int col, if (!CheckColumn(v)) { throw DimensionException(); } - for (int r = minrow; r <= maxrow; ++r) { - storage[index(r, col)] = v[r]; + for (int r = m_minrow; r <= m_maxrow; ++r) { + m_storage[index(r, col)] = v[r]; } } From 726a5f932d198805e6cea3a4212be964451d5975 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Thu, 18 Dec 2025 09:28:24 +0000 Subject: [PATCH 03/16] Satisfy standards checks --- src/core/matrix.imp | 2 -- src/core/recarray.h | 5 ++--- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/core/matrix.imp b/src/core/matrix.imp index 0fda5cbca8..3c1e189585 100644 --- a/src/core/matrix.imp +++ b/src/core/matrix.imp @@ -28,7 +28,6 @@ namespace Gambit { // Matrix: Constructors, destructors, constructive operators //------------------------------------------------------------------------- - template Matrix &Matrix::operator=(const T &c) { for (int i = MinRow(); i <= MaxRow(); i++) { @@ -303,7 +302,6 @@ template bool Matrix::operator==(const T &s) const return true; } - template bool Matrix::operator!=(const T &s) const { return !(*this == s); } // Information diff --git a/src/core/recarray.h b/src/core/recarray.h index bc0e50bdc4..208f3e5354 100644 --- a/src/core/recarray.h +++ b/src/core/recarray.h @@ -132,9 +132,8 @@ template class RectArray { } auto stride = row_stride(); - size_t ai = index(i, m_mincol); - size_t aj = index(j, m_mincol); - + const size_t ai = index(i, m_mincol); + const size_t aj = index(j, m_mincol); for (size_t k = 0; k < stride; ++k) { std::swap(m_storage[ai + k], m_storage[aj + k]); } From 33607cc63e7aa5e422aa00582d9e9d3268e0373f Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Thu, 18 Dec 2025 09:47:34 +0000 Subject: [PATCH 04/16] More standards fixups --- src/core/matrix.h | 18 ++++---------- src/core/matrix.imp | 1 - src/core/recarray.h | 60 +++++++++++++++++---------------------------- 3 files changed, 28 insertions(+), 51 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 5d96ac7426..760732a082 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -34,11 +34,7 @@ class SingularMatrixException final : public std::runtime_error { ~SingularMatrixException() noexcept override = default; }; -template Vector operator*(const Vector &, const Matrix &); - template class Matrix { - friend Vector operator* <>(const Vector &, const Matrix &); - RectArray m_data; public: @@ -115,7 +111,7 @@ template class Matrix { void CMultiply(const Vector &, Vector &) const; /// "in-place" row (transposed) multiply void RMultiply(const Vector &, Vector &) const; - Matrix operator*(const Matrix &) const; + Matrix operator*(const Matrix &) const; Vector operator*(const Vector &) const; Matrix operator*(const T &) const; Matrix &operator*=(const T &); @@ -138,7 +134,7 @@ template class Matrix { template Vector operator*(const Vector &, const Matrix &); -template template void Matrix::GetColumn(int col, Vector &v) const +template template void Matrix::GetColumn(int col, V &v) const { if (col < MinCol() || col > MaxCol()) { throw std::out_of_range("Index out of range in Matrix::GetColumn"); @@ -146,13 +142,12 @@ template template void Matrix::GetColumn(int col, Ve if (!CheckColumn(v)) { throw DimensionException(); } - for (int i = MinRow(); i <= MaxRow(); ++i) { v[i] = (*this)(i, col); } } -template template void Matrix::SetColumn(int col, const Vector &v) +template template void Matrix::SetColumn(int col, const V &v) { if (col < MinCol() || col > MaxCol()) { throw std::out_of_range("Index out of range in Matrix::SetColumn"); @@ -160,13 +155,12 @@ template template void Matrix::SetColumn(int col, co if (!CheckColumn(v)) { throw DimensionException(); } - for (int i = MinRow(); i <= MaxRow(); ++i) { (*this)(i, col) = v[i]; } } -template template void Matrix::GetRow(int row, Vector &v) const +template template void Matrix::GetRow(int row, V &v) const { if (row < MinRow() || row > MaxRow()) { throw std::out_of_range("Index out of range in Matrix::GetRow"); @@ -174,13 +168,12 @@ template template void Matrix::GetRow(int row, Vecto if (!CheckRow(v)) { throw DimensionException(); } - for (int j = MinCol(); j <= MaxCol(); ++j) { v[j] = (*this)(row, j); } } -template template void Matrix::SetRow(int row, const Vector &v) +template template void Matrix::SetRow(int row, const V &v) { if (row < MinRow() || row > MaxRow()) { throw std::out_of_range("Index out of range in Matrix::SetRow"); @@ -188,7 +181,6 @@ template template void Matrix::SetRow(int row, const if (!CheckRow(v)) { throw DimensionException(); } - for (int j = MinCol(); j <= MaxCol(); ++j) { (*this)(row, j) = v[j]; } diff --git a/src/core/matrix.imp b/src/core/matrix.imp index 3c1e189585..c0cbe3ec1b 100644 --- a/src/core/matrix.imp +++ b/src/core/matrix.imp @@ -184,7 +184,6 @@ template void Matrix::RMultiply(const Vector &in, Vector &out } // transposed (row) vector*matrix multiplication operator -// a friend function of Matrix template Vector operator*(const Vector &v, const Matrix &M) { if (!M.CheckColumn(v)) { diff --git a/src/core/recarray.h b/src/core/recarray.h index 208f3e5354..7acf31f6da 100644 --- a/src/core/recarray.h +++ b/src/core/recarray.h @@ -45,45 +45,46 @@ template class RectArray { /// @name Lifecycle //@{ RectArray() : m_minrow(1), m_maxrow(0), m_mincol(1), m_maxcol(0), m_storage() {} - RectArray(unsigned int nrows, unsigned int ncols); - RectArray(int minr, int maxr, int minc, int maxc) - : m_minrow(minr), m_maxrow(maxr), m_mincol(minc), m_maxcol(maxc), - m_storage((maxr >= minr && maxc >= minc) ? (maxr - minr + 1) * (maxc - minc + 1) : 0) + RectArray(const size_t nrows, const size_t ncols) : RectArray(1, nrows, 1, ncols) {} + RectArray(const int minrow, const int maxrow, const int mincol, const int maxcol) + : m_minrow(minrow), m_maxrow(maxrow), m_mincol(mincol), m_maxcol(maxcol), + m_storage((maxrow >= minrow && maxcol >= mincol) + ? (maxrow - minrow + 1) * (maxcol - mincol + 1) + : 0) { } RectArray(const RectArray &) = default; RectArray(RectArray &&) noexcept = default; - - virtual ~RectArray() = default; + ~RectArray() = default; RectArray &operator=(const RectArray &) = default; RectArray &operator=(RectArray &&) noexcept = default; //@} + /// @name Range checking functions; returns true only if valid index/size + //@{ /// check array for same row and column boundaries bool CheckBounds(const RectArray &m) const { return (m_minrow == m.m_minrow && m_maxrow == m.m_maxrow && m_mincol == m.m_mincol && m_maxcol == m.m_maxcol); } - /// @name Range checking functions; returns true only if valid index/size - //@{ /// check for correct row index - bool CheckRow(int row) const { return (m_minrow <= row && row <= m_maxrow); } + bool CheckRow(const int row) const { return (m_minrow <= row && row <= m_maxrow); } /// check row vector for correct column boundaries - template bool CheckRow(const Vector &v) const + template bool CheckRow(const V &v) const { - return (v.front_index() == m_mincol && v.back_index() == m_maxcol); + return v.front_index() == m_mincol && v.back_index() == m_maxcol; } /// check for correct column index - bool CheckColumn(int col) const { return (m_mincol <= col && col <= m_maxcol); } + bool CheckColumn(const int col) const { return (m_mincol <= col && col <= m_maxcol); } /// check column vector for correct row boundaries - template bool CheckColumn(const Vector &v) const + template bool CheckColumn(const V &v) const { return (v.front_index() == m_minrow && v.back_index() == m_maxrow); } /// check row and column indices - bool Check(int row, int col) const { return CheckRow(row) && CheckColumn(col); } + bool Check(const int row, const int col) const { return CheckRow(row) && CheckColumn(col); } //@} /// @name General data access @@ -131,7 +132,7 @@ template class RectArray { return; } - auto stride = row_stride(); + const auto stride = row_stride(); const size_t ai = index(i, m_mincol); const size_t aj = index(j, m_mincol); for (size_t k = 0; k < stride; ++k) { @@ -145,15 +146,6 @@ template class RectArray { //@} }; -//------------------------------------------------------------------------ -// RectArray: Constructors, destructor, constructive operators -//------------------------------------------------------------------------ - -template -RectArray::RectArray(unsigned int rows, unsigned int cols) : RectArray(1, rows, 1, cols) -{ -} - //------------------------------------------------------------------------ // RectArray: Row and column rotation //------------------------------------------------------------------------ @@ -163,12 +155,9 @@ template void RectArray::RotateUp(int lo, int hi) if (lo < m_minrow || hi < lo || m_maxrow < hi) { throw std::out_of_range("Index out of range in RectArray"); } - - auto stride = row_stride(); - - size_t first = index(lo, m_mincol); - size_t last = index(hi + 1, m_mincol); - + const auto stride = row_stride(); + const size_t first = index(lo, m_mincol); + const size_t last = index(hi + 1, m_mincol); std::rotate(m_storage.begin() + first, m_storage.begin() + first + stride, m_storage.begin() + last); } @@ -178,12 +167,9 @@ template void RectArray::RotateDown(int lo, int hi) if (lo < m_minrow || hi < lo || m_maxrow < hi) { throw std::out_of_range("Index out of range in RectArray"); } - - auto stride = row_stride(); - - size_t first = index(lo, m_mincol); - size_t last = index(hi + 1, m_mincol); - + const auto stride = row_stride(); + const size_t first = index(lo, m_mincol); + const size_t last = index(hi + 1, m_mincol); std::rotate(m_storage.begin() + first, m_storage.begin() + last - stride, m_storage.begin() + last); } @@ -200,7 +186,7 @@ template template void RectArray::GetRow(int row, Ve if (!CheckRow(v)) { throw DimensionException(); } - size_t base = index(row, m_mincol); + const size_t base = index(row, m_mincol); for (int c = m_mincol; c <= m_maxcol; ++c) { v[c] = m_storage[base + (c - m_mincol)]; } From 5f597131105f7198d8c98446b8f92708c322eff2 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 08:51:38 +0000 Subject: [PATCH 05/16] Removed unused Row() and Column() from Matrix --- src/core/matrix.h | 6 ------ src/core/matrix.imp | 19 ------------------- 2 files changed, 25 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 760732a082..5d9291ae6d 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -81,12 +81,6 @@ template class Matrix { template bool CheckColumn(const V &v) const { return m_data.CheckColumn(v); } bool CheckBounds(const Matrix &M) const { return m_data.CheckBounds(M.m_data); } - /// @name Extracting rows and columns - //@{ - Vector Row(int) const; - Vector Column(int) const; - //@} - /// @name Comparison operators //@{ bool operator==(const Matrix &) const; diff --git a/src/core/matrix.imp b/src/core/matrix.imp index c0cbe3ec1b..c2d093c88f 100644 --- a/src/core/matrix.imp +++ b/src/core/matrix.imp @@ -303,25 +303,6 @@ template bool Matrix::operator==(const T &s) const template bool Matrix::operator!=(const T &s) const { return !(*this == s); } -// Information - -template Vector Matrix::Row(int i) const -{ - Vector answer(MinCol(), MaxCol()); - for (int j = MinCol(); j <= MaxCol(); j++) { - answer[j] = (*this)(i, j); - } - return answer; -} - -template Vector Matrix::Column(int j) const -{ - Vector answer(MinRow(), MaxRow()); - for (int i = MinRow(); i <= MaxRow(); i++) { - answer[i] = (*this)(i, j); - } - return answer; -} // more complex functions From 6630e766e0363013221629a9625b3ac58c5c6a04 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 09:33:10 +0000 Subject: [PATCH 06/16] Create concept of element iteration in RectArray, and use it to simplify a number of Matrix arithmetic operations --- src/core/matrix.h | 99 +++++++++++++++++++++++- src/core/matrix.imp | 184 -------------------------------------------- src/core/recarray.h | 13 ++++ 3 files changed, 111 insertions(+), 185 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 5d9291ae6d..82ef6f871e 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -96,7 +96,7 @@ template class Matrix { Matrix &operator+=(const Matrix &); Matrix &operator-=(const Matrix &); - Matrix operator-(); + Matrix operator-() const; //@} /// @name Multiplicative operators @@ -128,6 +128,103 @@ template class Matrix { template Vector operator*(const Vector &, const Matrix &); +template Matrix &Matrix::operator=(const T &c) +{ + std::fill(m_data.elements_begin(), m_data.elements_end(), c); + return *this; +} + +template bool Matrix::operator==(const Matrix &M) const +{ + if (!this->CheckBounds(M)) { + throw DimensionException(); + } + return std::equal(m_data.elements_begin(), m_data.elements_end(), M.m_data.elements_begin()); +} + +template bool Matrix::operator!=(const Matrix &M) const { return !(*this == M); } + +template bool Matrix::operator==(const T &c) const +{ + return std::all_of(m_data.elements_begin(), m_data.elements_end(), + [&c](const auto &v) { return v == c; }); +} + +template bool Matrix::operator!=(const T &c) const { return !(*this == c); } + +template Matrix &Matrix::operator+=(const Matrix &M) +{ + if (!this->CheckBounds(M)) { + throw DimensionException(); + } + std::transform(m_data.elements_begin(), m_data.elements_end(), M.m_data.elements_begin(), + m_data.elements_begin(), std::plus<>()); + return *this; +} + +template Matrix Matrix::operator+(const Matrix &M) const +{ + Matrix tmp(*this); + tmp += M; + return tmp; +} + +template Matrix &Matrix::operator-=(const Matrix &M) +{ + if (!this->CheckBounds(M)) { + throw DimensionException(); + } + std::transform(m_data.elements_begin(), m_data.elements_end(), M.m_data.elements_begin(), + m_data.elements_begin(), std::minus<>()); + return *this; +} + +template Matrix Matrix::operator-(const Matrix &M) const +{ + Matrix tmp(*this); + tmp -= M; + return tmp; +} + +template Matrix Matrix::operator-() const +{ + Matrix tmp(*this); + std::transform(tmp.m_data.elements_begin(), tmp.m_data.elements_end(), + tmp.m_data.elements_begin(), std::negate<>()); + return tmp; +} + +template Matrix &Matrix::operator*=(const T &c) +{ + std::transform(m_data.elements_begin(), m_data.elements_end(), m_data.elements_begin(), + [&c](const T &v) { return v * c; }); + return *this; +} + +template Matrix Matrix::operator*(const T &c) const +{ + Matrix tmp(*this); + tmp *= c; + return tmp; +} + +template Matrix &Matrix::operator/=(const T &c) +{ + if (c == static_cast(0)) { + throw ZeroDivideException(); + } + std::transform(m_data.elements_begin(), m_data.elements_end(), m_data.elements_begin(), + [&c](const T &v) { return v / c; }); + return *this; +} + +template Matrix Matrix::operator/(const T &c) const +{ + Matrix tmp(*this); + tmp /= c; + return tmp; +} + template template void Matrix::GetColumn(int col, V &v) const { if (col < MinCol() || col > MaxCol()) { diff --git a/src/core/matrix.imp b/src/core/matrix.imp index c2d093c88f..4dcd79f0d7 100644 --- a/src/core/matrix.imp +++ b/src/core/matrix.imp @@ -24,99 +24,6 @@ namespace Gambit { -//------------------------------------------------------------------------- -// Matrix: Constructors, destructors, constructive operators -//------------------------------------------------------------------------- - -template Matrix &Matrix::operator=(const T &c) -{ - for (int i = MinRow(); i <= MaxRow(); i++) { - for (int j = MinCol(); j <= MaxCol(); j++) { - (*this)(i, j) = c; - } - } - return *this; -} - -template Matrix Matrix::operator-() -{ - Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); - for (int i = MinRow(); i <= MaxRow(); i++) { - for (int j = MinCol(); j <= MaxCol(); j++) { - tmp(i, j) = -(*this)(i, j); - } - } - return tmp; -} - -//------------------------------------------------------------------------- -// Matrix: Additive operators -//------------------------------------------------------------------------- - -template Matrix Matrix::operator+(const Matrix &M) const -{ - if (!this->CheckBounds(M)) { - throw DimensionException(); - } - - Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); - - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - tmp(i, j) = (*this)(i, j) + M(i, j); - } - } - - return tmp; -} - -template Matrix Matrix::operator-(const Matrix &M) const -{ - if (!this->CheckBounds(M)) { - throw DimensionException(); - } - - Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); - - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - tmp(i, j) = (*this)(i, j) - M(i, j); - } - } - - return tmp; -} - -template Matrix &Matrix::operator+=(const Matrix &M) -{ - if (!this->CheckBounds(M)) { - throw DimensionException(); - } - - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - (*this)(i, j) += M(i, j); - } - } - - return *this; -} - -template Matrix &Matrix::operator-=(const Matrix &M) -{ - if (!this->CheckBounds(M)) { - throw DimensionException(); - } - - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - (*this)(i, j) -= M(i, j); - } - } - - return *this; -} - //------------------------------------------------------------------------- // Matrix: Multiplicative operators //------------------------------------------------------------------------- @@ -194,60 +101,6 @@ template Vector operator*(const Vector &v, const Matrix &M) return tmp; } -template Matrix Matrix::operator*(const T &s) const -{ - Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); - - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - tmp(i, j) = (*this)(i, j) * s; - } - } - - return tmp; -} - -template Matrix &Matrix::operator*=(const T &s) -{ - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - (*this)(i, j) *= s; - } - } - return *this; -} - -template Matrix Matrix::operator/(const T &s) const -{ - if (s == (T)0) { - throw ZeroDivideException(); - } - - Matrix tmp(MinRow(), MaxRow(), MinCol(), MaxCol()); - - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - tmp(i, j) = (*this)(i, j) / s; - } - } - - return tmp; -} - -template Matrix &Matrix::operator/=(const T &s) -{ - if (s == (T)0) { - throw ZeroDivideException(); - } - - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - (*this)(i, j) /= s; - } - } - - return *this; -} //------------------------------------------------------------------------- // Matrix: Transpose @@ -266,43 +119,6 @@ template Matrix Matrix::Transpose() const return tmp; } -//------------------------------------------------------------------------- -// Matrix: Comparison operators -//------------------------------------------------------------------------- - -template bool Matrix::operator==(const Matrix &M) const -{ - if (!this->CheckBounds(M)) { - throw DimensionException(); - } - - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - if ((*this)(i, j) != M(i, j)) { - return false; - } - } - } - - return true; -} - -template bool Matrix::operator!=(const Matrix &M) const { return !(*this == M); } - -template bool Matrix::operator==(const T &s) const -{ - for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - if ((*this)(i, j) != s) { - return false; - } - } - } - return true; -} - -template bool Matrix::operator!=(const T &s) const { return !(*this == s); } - // more complex functions diff --git a/src/core/recarray.h b/src/core/recarray.h index 7acf31f6da..a1beda0be1 100644 --- a/src/core/recarray.h +++ b/src/core/recarray.h @@ -144,6 +144,19 @@ template class RectArray { template void GetColumn(int, Vector &) const; template void SetColumn(int, const Vector &); //@} + + /// @name Iteration + /// @{ + using element_iterator = typename std::vector::iterator; + using const_element_iterator = typename std::vector::const_iterator; + + element_iterator elements_begin() noexcept { return m_storage.begin(); } + element_iterator elements_end() noexcept { return m_storage.end(); } + const_element_iterator elements_begin() const noexcept { return m_storage.begin(); } + const_element_iterator elements_end() const noexcept { return m_storage.end(); } + const_element_iterator elements_cbegin() const noexcept { return m_storage.cbegin(); } + const_element_iterator elements_cend() const noexcept { return m_storage.cend(); } + /// @} }; //------------------------------------------------------------------------ From 76af5bacde0b8f717ba7056855c018afcc89effd Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 09:58:11 +0000 Subject: [PATCH 07/16] Develop row/column views for RectArray; use in simplifying implementation of CMultiply --- src/core/matrix.imp | 11 +----- src/core/recarray.h | 92 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 94 insertions(+), 9 deletions(-) diff --git a/src/core/matrix.imp b/src/core/matrix.imp index 4dcd79f0d7..be67d3836d 100644 --- a/src/core/matrix.imp +++ b/src/core/matrix.imp @@ -33,16 +33,9 @@ template void Matrix::CMultiply(const Vector &in, Vector &out if (!this->CheckRow(in) || !this->CheckColumn(out)) { throw DimensionException(); } - for (int i = MinRow(); i <= MaxRow(); ++i) { - T sum = (T)0; - - auto src2 = in.begin(); - for (int j = MinCol(); j <= MaxCol(); ++j) { - sum += (*this)(i, j) * *(src2++); - } - - out[i] = sum; + auto row = m_data.GetRowView(i); + out[i] = std::inner_product(row.begin(), row.end(), in.begin(), static_cast(0)); } } diff --git a/src/core/recarray.h b/src/core/recarray.h index a1beda0be1..da36b35cb7 100644 --- a/src/core/recarray.h +++ b/src/core/recarray.h @@ -42,6 +42,94 @@ template class RectArray { } public: + class RowView { + RectArray *m_array; + int m_row; + + public: + class iterator { + RowView *m_view; + int m_col; + + public: + using iterator_category = std::forward_iterator_tag; + using value_type = T; + using difference_type = int; + using pointer = T *; + using reference = T &; + + iterator(RowView *p_view, int p_col) : m_view(p_view), m_col(p_col) {} + + reference operator*() const { return (*m_view)[m_col]; } + + iterator &operator++() + { + ++m_col; + return *this; + } + + bool operator==(const iterator &p_other) const { return m_col == p_other.m_col; } + bool operator!=(const iterator &p_other) const { return !(*this == p_other); } + }; + + RowView(RectArray &p_array, int p_row) : m_array(&p_array), m_row(p_row) + { + if (!m_array->CheckRow(m_row)) { + throw std::out_of_range("RowView"); + } + } + T &operator[](int c) { return (*m_array)(m_row, c); } + const T &operator[](int c) const { return (*m_array)(m_row, c); } + int MinIndex() const { return m_array->MinCol(); } + int MaxIndex() const { return m_array->MaxCol(); } + iterator begin() { return iterator(this, MinIndex()); } + iterator end() { return iterator(this, MaxIndex() + 1); } + }; + + class ColumnView { + RectArray *m_array; + int m_col; + + public: + class iterator { + ColumnView *m_view; + int m_row; + + public: + using iterator_category = std::forward_iterator_tag; + using value_type = T; + using difference_type = int; + using pointer = T *; + using reference = T &; + + iterator(ColumnView *p_view, int p_row) : m_view(p_view), m_row(p_row) {} + + reference operator*() const { return (*m_view)[m_row]; } + + iterator &operator++() + { + ++m_row; + return *this; + } + + bool operator==(const iterator &p_other) const { return m_row == p_other.m_row; } + bool operator!=(const iterator &p_other) const { return !(*this == p_other); } + }; + + ColumnView(RectArray &p_array, int p_col) : m_array(&p_array), m_col(p_col) + { + if (!m_array->CheckColumn(m_col)) { + throw std::out_of_range("ColumnView"); + } + } + T &operator[](int r) { return (*m_array)(r, m_col); } + const T &operator[](int r) const { return (*m_array)(r, m_col); } + int MinIndex() const { return m_array->MinRow(); } + int MaxIndex() const { return m_array->MaxRow(); } + iterator begin() { return iterator(this, MinIndex()); } + iterator end() { return iterator(this, MaxIndex() + 1); } + }; + /// @name Lifecycle //@{ RectArray() : m_minrow(1), m_maxrow(0), m_mincol(1), m_maxcol(0), m_storage() {} @@ -139,6 +227,10 @@ template class RectArray { std::swap(m_storage[ai + k], m_storage[aj + k]); } } + RowView GetRowView(int r) { return RowView(*this, r); } + RowView GetRowView(int r) const { return RowView(const_cast(*this), r); } + ColumnView GetColumnView(int c) { return ColumnView(*this, c); } + ColumnView GetColumnView(int c) const { return ColumnView(const_cast(*this), c); } template void GetRow(int, Vector &) const; template void GetColumn(int, Vector &) const; From b695a24d91211f82c38b2b63d650850c265994d3 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 10:17:33 +0000 Subject: [PATCH 08/16] Rewrite matrix multiply --- src/core/matrix.imp | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/src/core/matrix.imp b/src/core/matrix.imp index be67d3836d..62b57a2c4f 100644 --- a/src/core/matrix.imp +++ b/src/core/matrix.imp @@ -35,7 +35,7 @@ template void Matrix::CMultiply(const Vector &in, Vector &out } for (int i = MinRow(); i <= MaxRow(); ++i) { auto row = m_data.GetRowView(i); - out[i] = std::inner_product(row.begin(), row.end(), in.begin(), static_cast(0)); + out[i] = std::inner_product(row.begin(), row.end(), in.begin(), T{0}); } } @@ -44,14 +44,13 @@ template Matrix Matrix::operator*(const Matrix &M) const if (MinCol() != M.MinRow() || MaxCol() != M.MaxRow()) { throw DimensionException(); } - Matrix tmp(MinRow(), MaxRow(), M.MinCol(), M.MaxCol()); - Vector column(M.MinRow(), M.MaxRow()); - Vector result(MinRow(), MaxRow()); - for (int j = M.MinCol(); j <= M.MaxCol(); j++) { - M.GetColumn(j, column); - CMultiply(column, result); - tmp.SetColumn(j, result); + for (int i = MinRow(); i <= MaxRow(); ++i) { + auto row = m_data.GetRowView(i); + for (int j = M.MinCol(); j <= M.MaxCol(); ++j) { + auto col = M.m_data.GetColumnView(j); + tmp(i, j) = std::inner_product(row.begin(), row.end(), col.begin(), T{0}); + } } return tmp; } @@ -73,12 +72,13 @@ template void Matrix::RMultiply(const Vector &in, Vector &out throw DimensionException(); } - out = (T)0; - + out = T{0}; for (int i = MinRow(); i <= MaxRow(); ++i) { - T k = in[i]; - for (int j = MinCol(); j <= MaxCol(); ++j) { - out[j] += (*this)(i, j) * k; + auto row = m_data.GetRowView(i); + const T k = in[i]; + auto dst = out.begin(); + for (auto it = row.begin(); it != row.end(); ++it, ++dst) { + *dst += (*it) * k; } } } From 40ccbd474fc88cc4948581385171d43a5c01ed25 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 10:51:36 +0000 Subject: [PATCH 09/16] Make Matrix header-only --- Makefile.am | 1 - src/core/matrix.cc | 8 +- src/core/matrix.h | 264 +++++++++++++++++++++++++++++++++++++++ src/core/matrix.imp | 297 -------------------------------------------- 4 files changed, 265 insertions(+), 305 deletions(-) delete mode 100644 src/core/matrix.imp diff --git a/Makefile.am b/Makefile.am index f5178efe5a..380a8d6f23 100644 --- a/Makefile.am +++ b/Makefile.am @@ -249,7 +249,6 @@ core_SOURCES = \ src/core/vector.h \ src/core/recarray.h \ src/core/matrix.h \ - src/core/matrix.imp \ src/core/integer.cc \ src/core/integer.h \ src/core/rational.cc \ diff --git a/src/core/matrix.cc b/src/core/matrix.cc index 29ec79b303..e8badc350c 100644 --- a/src/core/matrix.cc +++ b/src/core/matrix.cc @@ -20,8 +20,7 @@ // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. // -#include "core.h" -#include "matrix.imp" +#include "matrix.h" namespace Gambit { @@ -30,9 +29,4 @@ template class Matrix; template class Matrix; template class Matrix; -template Vector operator*(const Vector &, const Matrix &); -template Vector operator*(const Vector &, const Matrix &); -template Vector operator*(const Vector &, const Matrix &); -template Vector operator*(const Vector &, const Matrix &); - } // end namespace Gambit diff --git a/src/core/matrix.h b/src/core/matrix.h index 82ef6f871e..7055abeadd 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -25,6 +25,7 @@ #include "recarray.h" #include "vector.h" +#include "rational.h" namespace Gambit { @@ -277,6 +278,269 @@ template template void Matrix::SetRow(int row, const V &v } } +template void Matrix::CMultiply(const Vector &in, Vector &out) const +{ + if (!this->CheckRow(in) || !this->CheckColumn(out)) { + throw DimensionException(); + } + for (int i = MinRow(); i <= MaxRow(); ++i) { + auto row = m_data.GetRowView(i); + out[i] = std::inner_product(row.begin(), row.end(), in.begin(), T{0}); + } +} + +template Matrix Matrix::operator*(const Matrix &M) const +{ + if (MinCol() != M.MinRow() || MaxCol() != M.MaxRow()) { + throw DimensionException(); + } + Matrix tmp(MinRow(), MaxRow(), M.MinCol(), M.MaxCol()); + for (int i = MinRow(); i <= MaxRow(); ++i) { + auto row = m_data.GetRowView(i); + for (int j = M.MinCol(); j <= M.MaxCol(); ++j) { + auto col = M.m_data.GetColumnView(j); + tmp(i, j) = std::inner_product(row.begin(), row.end(), col.begin(), T{0}); + } + } + return tmp; +} + +template Vector Matrix::operator*(const Vector &v) const +{ + if (!this->CheckRow(v)) { + throw DimensionException(); + } + + Vector tmp(MinRow(), MaxRow()); + CMultiply(v, tmp); + return tmp; +} + +template void Matrix::RMultiply(const Vector &in, Vector &out) const +{ + if (!this->CheckColumn(in) || !this->CheckRow(out)) { + throw DimensionException(); + } + + out = T{0}; + for (int i = MinRow(); i <= MaxRow(); ++i) { + auto row = m_data.GetRowView(i); + const T k = in[i]; + auto dst = out.begin(); + for (auto it = row.begin(); it != row.end(); ++it, ++dst) { + *dst += (*it) * k; + } + } +} + +// transposed (row) vector*matrix multiplication operator +template Vector operator*(const Vector &v, const Matrix &M) +{ + if (!M.CheckColumn(v)) { + throw DimensionException(); + } + Vector tmp(M.MinCol(), M.MaxCol()); + M.RMultiply(v, tmp); + return tmp; +} + +template Matrix Matrix::Transpose() const +{ + Matrix tmp(MinCol(), MaxCol(), MinRow(), MaxRow()); + + for (int i = MinRow(); i <= MaxRow(); i++) { + for (int j = MinCol(); j <= MaxCol(); j++) { + tmp(j, i) = (*this)(i, j); + } + } + + return tmp; +} + +template void Matrix::MakeIdent() +{ + if (!IsSquare()) { + throw DimensionException(); + } + for (int i = MinRow(); i <= MaxRow(); i++) { + for (int j = MinCol(); j <= MaxCol(); j++) { + if (i == j) { + (*this)(i, j) = (T)1; + } + else { + (*this)(i, j) = (T)0; + } + } + } +} + +template void Matrix::Pivot(int row, int col) +{ + if (!this->CheckRow(row) || !this->CheckColumn(col)) { + throw std::out_of_range("Index out of range in Matrix::Pivot"); + } + if ((*this)(row, col) == (T)0) { + throw ZeroDivideException(); + } + + T inv = (T)1 / (*this)(row, col); + + // scale pivot row + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(row, j) *= inv; + } + + // eliminate column + for (int i = MinRow(); i <= MaxRow(); ++i) { + if (i != row) { + T mult = (*this)(i, col); + for (int j = MinCol(); j <= MaxCol(); ++j) { + (*this)(i, j) -= (*this)(row, j) * mult; + } + } + } +} + +template Matrix Matrix::Inverse() const +{ + if (!IsSquare()) { + throw DimensionException(); + } + + Matrix copy(*this); + Matrix inv(MaxRow(), MaxRow()); + + // initialize inverse matrix and prescale row vectors + for (int i = MinRow(); i <= MaxRow(); i++) { + T max = (T)0; + for (int j = MinCol(); j <= MaxCol(); j++) { + T abs = copy(i, j); + if (abs < (T)0) { + abs = -abs; + } + if (abs > max) { + max = abs; + } + } + + if (max == (T)0) { + throw SingularMatrixException(); + } + + T scale = (T)1 / max; + for (int j = MinCol(); j <= MaxCol(); j++) { + copy(i, j) *= scale; + if (i == j) { + inv(i, j) = scale; + } + else { + inv(i, j) = (T)0; + } + } + } + + for (int i = MinCol(); i <= MaxCol(); i++) { + // find pivot row + T max = copy(i, i); + if (max < (T)0) { + max = -max; + } + int row = i; + for (int j = i + 1; j <= MaxRow(); j++) { + T abs = copy(j, i); + if (abs < (T)0) { + abs = -abs; + } + if (abs > max) { + max = abs; + row = j; + } + } + + if (max <= (T)0) { + throw SingularMatrixException(); + } + + copy.SwitchRows(i, row); + inv.SwitchRows(i, row); + // scale pivot row + T factor = (T)1 / copy(i, i); + for (int k = MinCol(); k <= MaxCol(); k++) { + copy(i, k) *= factor; + inv(i, k) *= factor; + } + + // reduce other rows + for (int j = MinRow(); j <= MaxRow(); j++) { + if (j != i) { + T mult = copy(j, i); + for (int k = MinCol(); k <= MaxCol(); k++) { + copy(j, k) -= copy(i, k) * mult; + inv(j, k) -= inv(i, k) * mult; + } + } + } + } + + return inv; +} + +template T Matrix::Determinant() const +{ + if (!IsSquare()) { + throw DimensionException(); + } + + T factor = (T)1; + Matrix M(*this); + + for (int row = MinRow(); row <= MaxRow(); row++) { + + // Experience (as of 3/22/99) suggests that, in the interest of + // numerical stability, it might be best to do Gaussian + // elimination with respect to the row (of those feasible) + // whose entry has the largest absolute value. + int swap_row = row; + for (int i = row + 1; i <= MaxRow(); i++) { + if (abs(M(i, row)) > abs(M(swap_row, row))) { + swap_row = i; + } + } + + if (swap_row != row) { + M.SwitchRows(row, swap_row); + for (int j = MinCol(); j <= MaxCol(); j++) { + M(row, j) *= (T)-1; + } + } + + if (M(row, row) == (T)0) { + return (T)0; + } + + // now do row operations to clear the row'th column + // below the diagonal + for (int row1 = row + 1; row1 <= MaxRow(); row1++) { + factor = -M(row1, row) / M(row, row); + for (int i = MinCol(); i <= MaxCol(); i++) { + M(row1, i) += M(row, i) * factor; + } + } + } + + // finally we multiply the diagonal elements + T det = (T)1; + for (int row = MinRow(); row <= MaxRow(); row++) { + det *= M(row, row); + } + return det; +} + +extern template class Matrix; +extern template class Matrix; +extern template class Matrix; +extern template class Matrix; + } // end namespace Gambit #endif // GAMBIT_CORE_MATRIX_H diff --git a/src/core/matrix.imp b/src/core/matrix.imp deleted file mode 100644 index 62b57a2c4f..0000000000 --- a/src/core/matrix.imp +++ /dev/null @@ -1,297 +0,0 @@ -// -// This file is part of Gambit -// Copyright (c) 1994-2025, The Gambit Project (https://www.gambit-project.org) -// -// FILE: src/core/matrix.imp -// Implementation of matrix method functions -// -// This program is free software; you can redistribute it and/or modify -// it under the terms of the GNU General Public License as published by -// the Free Software Foundation; either version 2 of the License, or -// (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, -// but WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -// GNU General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program; if not, write to the Free Software -// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. -// - -#include "matrix.h" - -namespace Gambit { - -//------------------------------------------------------------------------- -// Matrix: Multiplicative operators -//------------------------------------------------------------------------- - -template void Matrix::CMultiply(const Vector &in, Vector &out) const -{ - if (!this->CheckRow(in) || !this->CheckColumn(out)) { - throw DimensionException(); - } - for (int i = MinRow(); i <= MaxRow(); ++i) { - auto row = m_data.GetRowView(i); - out[i] = std::inner_product(row.begin(), row.end(), in.begin(), T{0}); - } -} - -template Matrix Matrix::operator*(const Matrix &M) const -{ - if (MinCol() != M.MinRow() || MaxCol() != M.MaxRow()) { - throw DimensionException(); - } - Matrix tmp(MinRow(), MaxRow(), M.MinCol(), M.MaxCol()); - for (int i = MinRow(); i <= MaxRow(); ++i) { - auto row = m_data.GetRowView(i); - for (int j = M.MinCol(); j <= M.MaxCol(); ++j) { - auto col = M.m_data.GetColumnView(j); - tmp(i, j) = std::inner_product(row.begin(), row.end(), col.begin(), T{0}); - } - } - return tmp; -} - -template Vector Matrix::operator*(const Vector &v) const -{ - if (!this->CheckRow(v)) { - throw DimensionException(); - } - - Vector tmp(MinRow(), MaxRow()); - CMultiply(v, tmp); - return tmp; -} - -template void Matrix::RMultiply(const Vector &in, Vector &out) const -{ - if (!this->CheckColumn(in) || !this->CheckRow(out)) { - throw DimensionException(); - } - - out = T{0}; - for (int i = MinRow(); i <= MaxRow(); ++i) { - auto row = m_data.GetRowView(i); - const T k = in[i]; - auto dst = out.begin(); - for (auto it = row.begin(); it != row.end(); ++it, ++dst) { - *dst += (*it) * k; - } - } -} - -// transposed (row) vector*matrix multiplication operator -template Vector operator*(const Vector &v, const Matrix &M) -{ - if (!M.CheckColumn(v)) { - throw DimensionException(); - } - Vector tmp(M.MinCol(), M.MaxCol()); - M.RMultiply(v, tmp); - return tmp; -} - - -//------------------------------------------------------------------------- -// Matrix: Transpose -//------------------------------------------------------------------------- - -template Matrix Matrix::Transpose() const -{ - Matrix tmp(MinCol(), MaxCol(), MinRow(), MaxRow()); - - for (int i = MinRow(); i <= MaxRow(); i++) { - for (int j = MinCol(); j <= MaxCol(); j++) { - tmp(j, i) = (*this)(i, j); - } - } - - return tmp; -} - - -// more complex functions - -template void Matrix::MakeIdent() -{ - if (!IsSquare()) { - throw DimensionException(); - } - for (int i = MinRow(); i <= MaxRow(); i++) { - for (int j = MinCol(); j <= MaxCol(); j++) { - if (i == j) { - (*this)(i, j) = (T)1; - } - else { - (*this)(i, j) = (T)0; - } - } - } -} - -template void Matrix::Pivot(int row, int col) -{ - if (!this->CheckRow(row) || !this->CheckColumn(col)) { - throw std::out_of_range("Index out of range in Matrix::Pivot"); - } - if ((*this)(row, col) == (T)0) { - throw ZeroDivideException(); - } - - T inv = (T)1 / (*this)(row, col); - - // scale pivot row - for (int j = MinCol(); j <= MaxCol(); ++j) { - (*this)(row, j) *= inv; - } - - // eliminate column - for (int i = MinRow(); i <= MaxRow(); ++i) { - if (i != row) { - T mult = (*this)(i, col); - for (int j = MinCol(); j <= MaxCol(); ++j) { - (*this)(i, j) -= (*this)(row, j) * mult; - } - } - } -} - -template Matrix Matrix::Inverse() const -{ - if (!IsSquare()) { - throw DimensionException(); - } - - Matrix copy(*this); - Matrix inv(MaxRow(), MaxRow()); - - // initialize inverse matrix and prescale row vectors - for (int i = MinRow(); i <= MaxRow(); i++) { - T max = (T)0; - for (int j = MinCol(); j <= MaxCol(); j++) { - T abs = copy(i, j); - if (abs < (T)0) { - abs = -abs; - } - if (abs > max) { - max = abs; - } - } - - if (max == (T)0) { - throw SingularMatrixException(); - } - - T scale = (T)1 / max; - for (int j = MinCol(); j <= MaxCol(); j++) { - copy(i, j) *= scale; - if (i == j) { - inv(i, j) = scale; - } - else { - inv(i, j) = (T)0; - } - } - } - - for (int i = MinCol(); i <= MaxCol(); i++) { - // find pivot row - T max = copy(i, i); - if (max < (T)0) { - max = -max; - } - int row = i; - for (int j = i + 1; j <= MaxRow(); j++) { - T abs = copy(j, i); - if (abs < (T)0) { - abs = -abs; - } - if (abs > max) { - max = abs; - row = j; - } - } - - if (max <= (T)0) { - throw SingularMatrixException(); - } - - copy.SwitchRows(i, row); - inv.SwitchRows(i, row); - // scale pivot row - T factor = (T)1 / copy(i, i); - for (int k = MinCol(); k <= MaxCol(); k++) { - copy(i, k) *= factor; - inv(i, k) *= factor; - } - - // reduce other rows - for (int j = MinRow(); j <= MaxRow(); j++) { - if (j != i) { - T mult = copy(j, i); - for (int k = MinCol(); k <= MaxCol(); k++) { - copy(j, k) -= copy(i, k) * mult; - inv(j, k) -= inv(i, k) * mult; - } - } - } - } - - return inv; -} - -template T Matrix::Determinant() const -{ - if (!IsSquare()) { - throw DimensionException(); - } - - T factor = (T)1; - Matrix M(*this); - - for (int row = MinRow(); row <= MaxRow(); row++) { - - // Experience (as of 3/22/99) suggests that, in the interest of - // numerical stability, it might be best to do Gaussian - // elimination with respect to the row (of those feasible) - // whose entry has the largest absolute value. - int swap_row = row; - for (int i = row + 1; i <= MaxRow(); i++) { - if (abs(M(i, row)) > abs(M(swap_row, row))) { - swap_row = i; - } - } - - if (swap_row != row) { - M.SwitchRows(row, swap_row); - for (int j = MinCol(); j <= MaxCol(); j++) { - M(row, j) *= (T)-1; - } - } - - if (M(row, row) == (T)0) { - return (T)0; - } - - // now do row operations to clear the row'th column - // below the diagonal - for (int row1 = row + 1; row1 <= MaxRow(); row1++) { - factor = -M(row1, row) / M(row, row); - for (int i = MinCol(); i <= MaxCol(); i++) { - M(row1, i) += M(row, i) * factor; - } - } - } - - // finally we multiply the diagonal elements - T det = (T)1; - for (int row = MinRow(); row <= MaxRow(); row++) { - det *= M(row, row); - } - return det; -} - -} // end namespace Gambit From b485bcd9a4800b0094c7bd28c5020368c1888a24 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 11:07:49 +0000 Subject: [PATCH 10/16] Rearrange structure of matrix header --- src/core/matrix.h | 148 ++++++++++++++++++++++++---------------------- 1 file changed, 78 insertions(+), 70 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 7055abeadd..e2d96bb92f 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -85,15 +85,25 @@ template class Matrix { /// @name Comparison operators //@{ bool operator==(const Matrix &) const; - bool operator!=(const Matrix &) const; + bool operator!=(const Matrix &M) const { return !(*this == M); } bool operator==(const T &) const; - bool operator!=(const T &) const; + bool operator!=(const T &c) const { return !(*this == c); } //@} /// @name Additive operators //@{ - Matrix operator+(const Matrix &) const; - Matrix operator-(const Matrix &) const; + Matrix operator+(const Matrix &M) const + { + Matrix tmp(*this); + tmp += M; + return tmp; + } + Matrix operator-(const Matrix &M) const + { + Matrix tmp(*this); + tmp -= M; + return tmp; + } Matrix &operator+=(const Matrix &); Matrix &operator-=(const Matrix &); @@ -108,10 +118,20 @@ template class Matrix { void RMultiply(const Vector &, Vector &) const; Matrix operator*(const Matrix &) const; Vector operator*(const Vector &) const; - Matrix operator*(const T &) const; + Matrix operator*(const T &c) const + { + Matrix tmp(*this); + tmp *= c; + return tmp; + } Matrix &operator*=(const T &); - Matrix operator/(const T &) const; + Matrix operator/(const T &c) const + { + Matrix tmp(*this); + tmp /= c; + return tmp; + } Matrix &operator/=(const T &); //@ @@ -127,14 +147,16 @@ template class Matrix { T Determinant() const; }; -template Vector operator*(const Vector &, const Matrix &); - template Matrix &Matrix::operator=(const T &c) { std::fill(m_data.elements_begin(), m_data.elements_end(), c); return *this; } +// ---------------------------------------------------------------------------- +// Implementation of element-wise operations +// ---------------------------------------------------------------------------- + template bool Matrix::operator==(const Matrix &M) const { if (!this->CheckBounds(M)) { @@ -143,16 +165,12 @@ template bool Matrix::operator==(const Matrix &M) const return std::equal(m_data.elements_begin(), m_data.elements_end(), M.m_data.elements_begin()); } -template bool Matrix::operator!=(const Matrix &M) const { return !(*this == M); } - template bool Matrix::operator==(const T &c) const { return std::all_of(m_data.elements_begin(), m_data.elements_end(), [&c](const auto &v) { return v == c; }); } -template bool Matrix::operator!=(const T &c) const { return !(*this == c); } - template Matrix &Matrix::operator+=(const Matrix &M) { if (!this->CheckBounds(M)) { @@ -163,13 +181,6 @@ template Matrix &Matrix::operator+=(const Matrix &M) return *this; } -template Matrix Matrix::operator+(const Matrix &M) const -{ - Matrix tmp(*this); - tmp += M; - return tmp; -} - template Matrix &Matrix::operator-=(const Matrix &M) { if (!this->CheckBounds(M)) { @@ -180,13 +191,6 @@ template Matrix &Matrix::operator-=(const Matrix &M) return *this; } -template Matrix Matrix::operator-(const Matrix &M) const -{ - Matrix tmp(*this); - tmp -= M; - return tmp; -} - template Matrix Matrix::operator-() const { Matrix tmp(*this); @@ -202,13 +206,6 @@ template Matrix &Matrix::operator*=(const T &c) return *this; } -template Matrix Matrix::operator*(const T &c) const -{ - Matrix tmp(*this); - tmp *= c; - return tmp; -} - template Matrix &Matrix::operator/=(const T &c) { if (c == static_cast(0)) { @@ -219,12 +216,9 @@ template Matrix &Matrix::operator/=(const T &c) return *this; } -template Matrix Matrix::operator/(const T &c) const -{ - Matrix tmp(*this); - tmp /= c; - return tmp; -} +// ---------------------------------------------------------------------------- +// Implementation of row/column access +// ---------------------------------------------------------------------------- template template void Matrix::GetColumn(int col, V &v) const { @@ -278,6 +272,10 @@ template template void Matrix::SetRow(int row, const V &v } } +// ---------------------------------------------------------------------------- +// Implementation of linear algebra concepts +// ---------------------------------------------------------------------------- + template void Matrix::CMultiply(const Vector &in, Vector &out) const { if (!this->CheckRow(in) || !this->CheckColumn(out)) { @@ -289,20 +287,21 @@ template void Matrix::CMultiply(const Vector &in, Vector &out } } -template Matrix Matrix::operator*(const Matrix &M) const +template void Matrix::RMultiply(const Vector &in, Vector &out) const { - if (MinCol() != M.MinRow() || MaxCol() != M.MaxRow()) { + if (!this->CheckColumn(in) || !this->CheckRow(out)) { throw DimensionException(); } - Matrix tmp(MinRow(), MaxRow(), M.MinCol(), M.MaxCol()); + + out = T{0}; for (int i = MinRow(); i <= MaxRow(); ++i) { auto row = m_data.GetRowView(i); - for (int j = M.MinCol(); j <= M.MaxCol(); ++j) { - auto col = M.m_data.GetColumnView(j); - tmp(i, j) = std::inner_product(row.begin(), row.end(), col.begin(), T{0}); + const T k = in[i]; + auto dst = out.begin(); + for (auto it = row.begin(); it != row.end(); ++it, ++dst) { + *dst += (*it) * k; } } - return tmp; } template Vector Matrix::operator*(const Vector &v) const @@ -310,46 +309,33 @@ template Vector Matrix::operator*(const Vector &v) const if (!this->CheckRow(v)) { throw DimensionException(); } - Vector tmp(MinRow(), MaxRow()); CMultiply(v, tmp); return tmp; } -template void Matrix::RMultiply(const Vector &in, Vector &out) const +template Matrix Matrix::operator*(const Matrix &M) const { - if (!this->CheckColumn(in) || !this->CheckRow(out)) { + if (MinCol() != M.MinRow() || MaxCol() != M.MaxRow()) { throw DimensionException(); } - - out = T{0}; + Matrix tmp(MinRow(), MaxRow(), M.MinCol(), M.MaxCol()); for (int i = MinRow(); i <= MaxRow(); ++i) { auto row = m_data.GetRowView(i); - const T k = in[i]; - auto dst = out.begin(); - for (auto it = row.begin(); it != row.end(); ++it, ++dst) { - *dst += (*it) * k; + for (int j = M.MinCol(); j <= M.MaxCol(); ++j) { + auto col = M.m_data.GetColumnView(j); + tmp(i, j) = std::inner_product(row.begin(), row.end(), col.begin(), T{0}); } } -} - -// transposed (row) vector*matrix multiplication operator -template Vector operator*(const Vector &v, const Matrix &M) -{ - if (!M.CheckColumn(v)) { - throw DimensionException(); - } - Vector tmp(M.MinCol(), M.MaxCol()); - M.RMultiply(v, tmp); return tmp; } template Matrix Matrix::Transpose() const { - Matrix tmp(MinCol(), MaxCol(), MinRow(), MaxRow()); + Matrix tmp(MinCol(), MaxCol(), MinRow(), MaxRow()); - for (int i = MinRow(); i <= MaxRow(); i++) { - for (int j = MinCol(); j <= MaxCol(); j++) { + for (int i = MinRow(); i <= MaxRow(); ++i) { + for (int j = MinCol(); j <= MaxCol(); ++j) { tmp(j, i) = (*this)(i, j); } } @@ -357,6 +343,10 @@ template Matrix Matrix::Transpose() const return tmp; } +// ---------------------------------------------------------------------------- +// Implementation of additional operations +// ---------------------------------------------------------------------------- + template void Matrix::MakeIdent() { if (!IsSquare()) { @@ -536,10 +526,28 @@ template T Matrix::Determinant() const return det; } +// ---------------------------------------------------------------------------- +// Implementation of operators +// ---------------------------------------------------------------------------- + +template Vector operator*(const Vector &v, const Matrix &M) +{ + if (!M.CheckColumn(v)) { + throw DimensionException(); + } + Vector tmp(M.MinCol(), M.MaxCol()); + M.RMultiply(v, tmp); + return tmp; +} + +// ---------------------------------------------------------------------------- +// Explicit instantiations +// ---------------------------------------------------------------------------- + +extern template class Matrix; extern template class Matrix; -extern template class Matrix; extern template class Matrix; -extern template class Matrix; +extern template class Matrix; } // end namespace Gambit From 5ef28efb832458f2a2706b492d482d78c15097e4 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 11:15:39 +0000 Subject: [PATCH 11/16] Move hapax MakeIdent into path-following --- src/core/matrix.h | 19 ------------------- src/solvers/logit/path.cc | 10 +++++++++- 2 files changed, 9 insertions(+), 20 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index e2d96bb92f..d5c0f59232 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -138,8 +138,6 @@ template class Matrix { /// @name Other operations //@{ Matrix Transpose() const; - /// Set matrix to identity matrix - void MakeIdent(); void Pivot(int, int); //@} @@ -347,23 +345,6 @@ template Matrix Matrix::Transpose() const // Implementation of additional operations // ---------------------------------------------------------------------------- -template void Matrix::MakeIdent() -{ - if (!IsSquare()) { - throw DimensionException(); - } - for (int i = MinRow(); i <= MaxRow(); i++) { - for (int j = MinCol(); j <= MaxCol(); j++) { - if (i == j) { - (*this)(i, j) = (T)1; - } - else { - (*this)(i, j) = (T)0; - } - } - } -} - template void Matrix::Pivot(int row, int col) { if (!this->CheckRow(row) || !this->CheckColumn(col)) { diff --git a/src/solvers/logit/path.cc b/src/solvers/logit/path.cc index 2bfc2ef043..03d9a7d1b7 100644 --- a/src/solvers/logit/path.cc +++ b/src/solvers/logit/path.cc @@ -70,9 +70,17 @@ void Givens(Matrix &b, Matrix &q, double &c1, double &c2, int l1 c2 = 0.0; } +void SetAsIdentity(Matrix &M) +{ + M = 0.0; + for (int i = M.MinRow(); i <= M.MaxRow(); ++i) { + M(i, i) = 1.0; + } +} + void QRDecomp(Matrix &b, Matrix &q) { - q.MakeIdent(); + SetAsIdentity(q); for (size_t m = 1; m <= b.NumColumns(); m++) { for (size_t k = m + 1; k <= b.NumRows(); k++) { Givens(b, q, b(m, m), b(k, m), m, k, m + 1); From 36af1a625d49e45da4e7a2fa9a74cc5e2efcb6c5 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 11:25:33 +0000 Subject: [PATCH 12/16] Remove unused Matrix::Pivot operation --- src/core/matrix.h | 38 ++++++++------------------------------ 1 file changed, 8 insertions(+), 30 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index d5c0f59232..4ecd227ef3 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -138,7 +138,6 @@ template class Matrix { /// @name Other operations //@{ Matrix Transpose() const; - void Pivot(int, int); //@} Matrix Inverse() const; @@ -333,8 +332,14 @@ template Matrix Matrix::Transpose() const Matrix tmp(MinCol(), MaxCol(), MinRow(), MaxRow()); for (int i = MinRow(); i <= MaxRow(); ++i) { - for (int j = MinCol(); j <= MaxCol(); ++j) { - tmp(j, i) = (*this)(i, j); + auto src_row = m_data.GetRowView(i); + auto dst_col = tmp.m_data.GetColumnView(i); + + auto src = src_row.begin(); + auto dst = dst_col.begin(); + + for (; src != src_row.end(); ++src, ++dst) { + *dst = *src; } } @@ -345,33 +350,6 @@ template Matrix Matrix::Transpose() const // Implementation of additional operations // ---------------------------------------------------------------------------- -template void Matrix::Pivot(int row, int col) -{ - if (!this->CheckRow(row) || !this->CheckColumn(col)) { - throw std::out_of_range("Index out of range in Matrix::Pivot"); - } - if ((*this)(row, col) == (T)0) { - throw ZeroDivideException(); - } - - T inv = (T)1 / (*this)(row, col); - - // scale pivot row - for (int j = MinCol(); j <= MaxCol(); ++j) { - (*this)(row, j) *= inv; - } - - // eliminate column - for (int i = MinRow(); i <= MaxRow(); ++i) { - if (i != row) { - T mult = (*this)(i, col); - for (int j = MinCol(); j <= MaxCol(); ++j) { - (*this)(i, j) -= (*this)(row, j) * mult; - } - } - } -} - template Matrix Matrix::Inverse() const { if (!IsSquare()) { From 689dbcefaf470614bdfb03fd2d955810ffde9a1a Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 12:27:15 +0000 Subject: [PATCH 13/16] Rework Matrix::Inverse --- src/core/matrix.h | 98 +++++++++++++++++++++++---------------------- src/core/recarray.h | 4 ++ 2 files changed, 55 insertions(+), 47 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 4ecd227ef3..7add578e2f 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -355,78 +355,82 @@ template Matrix Matrix::Inverse() const if (!IsSquare()) { throw DimensionException(); } + const int rmin = MinRow(); + const int rmax = MaxRow(); + const int cmin = MinCol(); + const int cmax = MaxCol(); Matrix copy(*this); - Matrix inv(MaxRow(), MaxRow()); + Matrix inv(rmin, rmax, cmin, cmax); + + inv = T{0}; // initialize inverse matrix and prescale row vectors - for (int i = MinRow(); i <= MaxRow(); i++) { - T max = (T)0; - for (int j = MinCol(); j <= MaxCol(); j++) { - T abs = copy(i, j); - if (abs < (T)0) { - abs = -abs; - } - if (abs > max) { - max = abs; - } - } + for (int i = rmin; i <= rmax; ++i) { + auto copy_row = copy.m_data.GetRowView(i); + auto inv_row = inv.m_data.GetRowView(i); - if (max == (T)0) { + T max = maximize_function(copy_row, [](const T &v) { return abs(v); }); + if (max == T{0}) { throw SingularMatrixException(); } - - T scale = (T)1 / max; - for (int j = MinCol(); j <= MaxCol(); j++) { - copy(i, j) *= scale; - if (i == j) { - inv(i, j) = scale; - } - else { - inv(i, j) = (T)0; - } + const T scale = T{1} / max; + for (auto &v : copy_row) { + v *= scale; } + inv_row[i] = scale; } - for (int i = MinCol(); i <= MaxCol(); i++) { + for (int i = cmin; i <= cmax; ++i) { // find pivot row - T max = copy(i, i); - if (max < (T)0) { - max = -max; - } + auto col_i = copy.m_data.GetColumnView(i); + T max = abs(col_i[i]); int row = i; - for (int j = i + 1; j <= MaxRow(); j++) { - T abs = copy(j, i); - if (abs < (T)0) { - abs = -abs; - } - if (abs > max) { - max = abs; + for (int j = i + 1; j <= rmax; ++j) { + const T v = abs(col_i[j]); + if (v > max) { + max = v; row = j; } } - if (max <= (T)0) { + if (max <= T{0}) { throw SingularMatrixException(); } copy.SwitchRows(i, row); inv.SwitchRows(i, row); // scale pivot row - T factor = (T)1 / copy(i, i); - for (int k = MinCol(); k <= MaxCol(); k++) { - copy(i, k) *= factor; - inv(i, k) *= factor; + T factor = T{1} / copy(i, i); + auto copy_row = copy.m_data.GetRowView(i); + auto inv_row = inv.m_data.GetRowView(i); + auto copy_it = copy_row.begin(); + auto inv_it = inv_row.begin(); + for (; copy_it != copy_row.end(); ++copy_it, ++inv_it) { + *copy_it *= factor; + *inv_it *= factor; } // reduce other rows - for (int j = MinRow(); j <= MaxRow(); j++) { - if (j != i) { - T mult = copy(j, i); - for (int k = MinCol(); k <= MaxCol(); k++) { - copy(j, k) -= copy(i, k) * mult; - inv(j, k) -= inv(i, k) * mult; - } + auto pivot_copy_row = copy.m_data.GetRowView(i); + auto pivot_inv_row = inv.m_data.GetRowView(i); + + for (int j = rmin; j <= rmax; ++j) { + if (j == i) { + continue; + } + auto row_copy = copy.m_data.GetRowView(j); + auto row_inv = inv.m_data.GetRowView(j); + const T mult = row_copy[i]; + auto pivot_copy_it = pivot_copy_row.begin(); + auto pivot_inv_it = pivot_inv_row.begin(); + auto row_copy_it = row_copy.begin(); + auto row_inv_it = row_inv.begin(); + + for (; pivot_copy_it != pivot_copy_row.end(); + ++pivot_copy_it, ++pivot_inv_it, ++row_copy_it, ++row_inv_it) { + *row_copy_it -= (*pivot_copy_it) * mult; + *row_inv_it -= (*pivot_inv_it) * mult; } } } diff --git a/src/core/recarray.h b/src/core/recarray.h index da36b35cb7..2e6e3dd885 100644 --- a/src/core/recarray.h +++ b/src/core/recarray.h @@ -83,7 +83,9 @@ template class RectArray { int MinIndex() const { return m_array->MinCol(); } int MaxIndex() const { return m_array->MaxCol(); } iterator begin() { return iterator(this, MinIndex()); } + iterator begin() const { return iterator(const_cast(this), MinIndex()); } iterator end() { return iterator(this, MaxIndex() + 1); } + iterator end() const { return iterator(const_cast(this), MaxIndex() + 1); } }; class ColumnView { @@ -127,7 +129,9 @@ template class RectArray { int MinIndex() const { return m_array->MinRow(); } int MaxIndex() const { return m_array->MaxRow(); } iterator begin() { return iterator(this, MinIndex()); } + iterator begin() const { return iterator(const_cast(this), MinIndex()); } iterator end() { return iterator(this, MaxIndex() + 1); } + iterator end() const { return iterator(const_cast(this), MaxIndex() + 1); } }; /// @name Lifecycle From 7f9bba0e04d7374d1be2bbf35944b7ece437acbb Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 12:41:54 +0000 Subject: [PATCH 14/16] Rework Matrix::Determinant --- src/core/matrix.h | 39 +++++++++++++++++++++++---------------- 1 file changed, 23 insertions(+), 16 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 7add578e2f..5ffd39531f 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -443,47 +443,54 @@ template T Matrix::Determinant() const if (!IsSquare()) { throw DimensionException(); } + const int rmin = MinRow(); + const int rmax = MaxRow(); - T factor = (T)1; Matrix M(*this); - for (int row = MinRow(); row <= MaxRow(); row++) { - + for (int row = rmin; row <= rmax; ++row) { // Experience (as of 3/22/99) suggests that, in the interest of // numerical stability, it might be best to do Gaussian // elimination with respect to the row (of those feasible) // whose entry has the largest absolute value. int swap_row = row; - for (int i = row + 1; i <= MaxRow(); i++) { - if (abs(M(i, row)) > abs(M(swap_row, row))) { + T max = abs(M(row, row)); + for (int i = row + 1; i <= rmax; ++i) { + const T v = abs(M(i, row)); + if (v > max) { + max = v; swap_row = i; } } if (swap_row != row) { M.SwitchRows(row, swap_row); - for (int j = MinCol(); j <= MaxCol(); j++) { - M(row, j) *= (T)-1; + auto pivot_row = M.m_data.GetRowView(row); + for (auto &v : pivot_row) { + v = -v; } } - if (M(row, row) == (T)0) { - return (T)0; + if (M(row, row) == T{0}) { + return T{0}; } - // now do row operations to clear the row'th column // below the diagonal - for (int row1 = row + 1; row1 <= MaxRow(); row1++) { - factor = -M(row1, row) / M(row, row); - for (int i = MinCol(); i <= MaxCol(); i++) { - M(row1, i) += M(row, i) * factor; + auto pivot_row = M.m_data.GetRowView(row); + for (int row1 = row + 1; row1 <= rmax; ++row1) { + auto elim_row = M.m_data.GetRowView(row1); + const T factor = -elim_row[row] / pivot_row[row]; + auto pivot_it = pivot_row.begin(); + auto elim_it = elim_row.begin(); + for (; pivot_it != pivot_row.end(); ++pivot_it, ++elim_it) { + *elim_it += (*pivot_it) * factor; } } } // finally we multiply the diagonal elements - T det = (T)1; - for (int row = MinRow(); row <= MaxRow(); row++) { + T det = T{1}; + for (int row = rmin; row <= rmax; ++row) { det *= M(row, row); } return det; From a45d075ad3045d3d1e3b6a6347fcffe711afef1e Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 13:17:24 +0000 Subject: [PATCH 15/16] Tidy up Matrix documentation. --- src/core/matrix.h | 137 +++++++++++++++++++++++++++++++++++----------- 1 file changed, 105 insertions(+), 32 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 5ffd39531f..71daf15662 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -35,12 +35,21 @@ class SingularMatrixException final : public std::runtime_error { ~SingularMatrixException() noexcept override = default; }; +/// @brief Dense rectangular matrix with arbitrary integer index ranges. +/// +/// @tparam T Scalar element type. Must support arithmetic operations, +/// comparison with zero, and abs(T). +/// +/// @note Inverse and Determinant use legacy Gaussian elimination algorithms. +/// They are known not to be numerically optimal for near-singular matrices. +/// Current behaviour is temporarily preserved for historical compatibility. template class Matrix { RectArray m_data; public: /// @name Lifecycle - //@{ + /// Constructors, assignment, and destruction + ///@{ Matrix() = default; Matrix(unsigned int rows, unsigned int cols) : m_data(rows, cols) {} Matrix(unsigned int rows, unsigned int cols, int minrows) @@ -55,23 +64,35 @@ template class Matrix { Matrix &operator=(const Matrix &) = default; Matrix &operator=(Matrix &&) noexcept = default; Matrix &operator=(const T &); - //@} + ///@} - // element access + /// @brief Access matrix element at (row, column) with bounds checking + /// + /// @param r Row index + /// @param c Column index + /// @throws std::out_of_range on invalid index. T &operator()(int r, int c) { return m_data(r, c); } + /// @copydoc operator(int, int) const T &operator()(int r, int c) const { return m_data(r, c); } - // bounds / dimensions + /// @brief Lowest valid row index int MinRow() const { return m_data.MinRow(); } + /// @brief Highest valid row index int MaxRow() const { return m_data.MaxRow(); } + /// @brief Lowest valid column index int MinCol() const { return m_data.MinCol(); } + /// @brief Highest valid column index int MaxCol() const { return m_data.MaxCol(); } + /// @brief Number of rows in the matrix int NumRows() const { return m_data.NumRows(); } + /// @brief Number of columns in the matrix int NumColumns() const { return m_data.NumColumns(); } + /// @brief Check if the matrix is a square matrix (num rows == num columns) bool IsSquare() const { return MinRow() == MinCol() && MaxRow() == MaxCol(); } - // row ops used internally + /// @name Row and column helpers + ///@{ void SwitchRows(int i, int j) { m_data.SwitchRows(i, j); } template void GetColumn(int j, V &) const; template void SetColumn(int j, const V &); @@ -81,17 +102,24 @@ template class Matrix { template bool CheckRow(const V &v) const { return m_data.CheckRow(v); } template bool CheckColumn(const V &v) const { return m_data.CheckColumn(v); } bool CheckBounds(const Matrix &M) const { return m_data.CheckBounds(M.m_data); } + ///@} /// @name Comparison operators - //@{ + /// Element-wise comparisons + ///@{ bool operator==(const Matrix &) const; bool operator!=(const Matrix &M) const { return !(*this == M); } bool operator==(const T &) const; bool operator!=(const T &c) const { return !(*this == c); } - //@} - - /// @name Additive operators - //@{ + ///@} + + /// @name Arithmetic operators + /// Element-wise and scalar arithmetic + /// + /// All matrix-matrix operations require identical row and column bounds and + /// operate component-by-component. Scalar operations apply uniformly to all + /// elements. + ///@{ Matrix operator+(const Matrix &M) const { Matrix tmp(*this); @@ -108,16 +136,7 @@ template class Matrix { Matrix &operator-=(const Matrix &); Matrix operator-() const; - //@} - - /// @name Multiplicative operators - //@{ - /// "in-place" column multiply - void CMultiply(const Vector &, Vector &) const; - /// "in-place" row (transposed) multiply - void RMultiply(const Vector &, Vector &) const; - Matrix operator*(const Matrix &) const; - Vector operator*(const Vector &) const; + Matrix operator*(const T &c) const { Matrix tmp(*this); @@ -133,15 +152,59 @@ template class Matrix { return tmp; } Matrix &operator/=(const T &); - //@ + ///@} + + /// @name Linear algebra operations + /// Some primitives for doing linear algebra. + ///@{ + + /// @brief Multiply a matrix by a column vector + /// + /// Computes p_output = (*this) * p_input, where @p p_input is interpreted + /// as a column vector and @p p_output receives the resulting column vector. + /// + /// @param p_input Input column vector + /// @param p_output Output column vector + /// @throws DimensionException if dimensions are incompatible + void CMultiply(const Vector &p_input, Vector &p_output) const; + + /// @brief Multiply a row vector by the matrix + /// + /// Computes p_output = p_input * (*this), where @p p_input is interpreted + /// as a row vector and @p p_output receives the resulting row vector. + /// + /// @param p_input Input row vector + /// @param p_output Output row vector + /// @throws DimensionException if dimensions are incompatible + void RMultiply(const Vector &p_input, Vector &p_output) const; + + /// @brief Multiply matrix by a column vector. + /// + /// Equivalent to CMultiply(v, result) + /// + /// @param v Input column vector + /// @return Resulting column vector + /// @throws DimensionException if dimensions are incompatible + Vector operator*(const Vector &v) const; + + /// @brief Matrix-matrix multiplication + /// + /// Computes the product (*this) * M; the number of columns of this matrix + /// must equal the number of rows in @param M . + /// + /// @param M Matrix to multiply this with + /// @return Resulting matrix + /// @throws DimensionException if dimensions are incompatible + Matrix operator*(const Matrix &M) const; + + ///@ /// @name Other operations - //@{ + ///@{ Matrix Transpose() const; - //@} - Matrix Inverse() const; T Determinant() const; + ///@} }; template Matrix &Matrix::operator=(const T &c) @@ -273,28 +336,28 @@ template template void Matrix::SetRow(int row, const V &v // Implementation of linear algebra concepts // ---------------------------------------------------------------------------- -template void Matrix::CMultiply(const Vector &in, Vector &out) const +template void Matrix::CMultiply(const Vector &p_input, Vector &p_output) const { - if (!this->CheckRow(in) || !this->CheckColumn(out)) { + if (!this->CheckRow(p_input) || !this->CheckColumn(p_output)) { throw DimensionException(); } for (int i = MinRow(); i <= MaxRow(); ++i) { auto row = m_data.GetRowView(i); - out[i] = std::inner_product(row.begin(), row.end(), in.begin(), T{0}); + p_output[i] = std::inner_product(row.begin(), row.end(), p_input.begin(), T{0}); } } -template void Matrix::RMultiply(const Vector &in, Vector &out) const +template void Matrix::RMultiply(const Vector &p_input, Vector &p_output) const { - if (!this->CheckColumn(in) || !this->CheckRow(out)) { + if (!this->CheckColumn(p_input) || !this->CheckRow(p_output)) { throw DimensionException(); } - out = T{0}; + p_output = T{0}; for (int i = MinRow(); i <= MaxRow(); ++i) { auto row = m_data.GetRowView(i); - const T k = in[i]; - auto dst = out.begin(); + const T k = p_input[i]; + auto dst = p_output.begin(); for (auto it = row.begin(); it != row.end(); ++it, ++dst) { *dst += (*it) * k; } @@ -359,6 +422,7 @@ template Matrix Matrix::Inverse() const const int rmax = MaxRow(); const int cmin = MinCol(); const int cmax = MaxCol(); + using Gambit::abs; Matrix copy(*this); Matrix inv(rmin, rmax, cmin, cmax); @@ -445,6 +509,7 @@ template T Matrix::Determinant() const } const int rmin = MinRow(); const int rmax = MaxRow(); + using Gambit::abs; Matrix M(*this); @@ -500,6 +565,14 @@ template T Matrix::Determinant() const // Implementation of operators // ---------------------------------------------------------------------------- +/// @brief Multiple a row vector by a matrix +/// +/// Computes v * M, where @param v is interpreted as a row vector +/// +/// @param v The row vector +/// @param M The matrix to multiply with +/// @throws DimensionException if dimensions are incompatible +/// @sa Matrix::RMultiply template Vector operator*(const Vector &v, const Matrix &M) { if (!M.CheckColumn(v)) { From df4de2f333491202878c6a097b389b18b560bb4b Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Mon, 29 Dec 2025 13:43:07 +0000 Subject: [PATCH 16/16] Add missing header for std::negate --- src/core/matrix.h | 37 +++++++++++++++++++++---------------- src/core/recarray.h | 12 ++++++------ 2 files changed, 27 insertions(+), 22 deletions(-) diff --git a/src/core/matrix.h b/src/core/matrix.h index 71daf15662..955a7e41a5 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -23,6 +23,8 @@ #ifndef GAMBIT_CORE_MATRIX_H #define GAMBIT_CORE_MATRIX_H +#include + #include "recarray.h" #include "vector.h" #include "rational.h" @@ -98,10 +100,13 @@ template class Matrix { template void SetColumn(int j, const V &); template void GetRow(int row, V &) const; template void SetRow(int row, const V &); - // vector helpers used internally - template bool CheckRow(const V &v) const { return m_data.CheckRow(v); } - template bool CheckColumn(const V &v) const { return m_data.CheckColumn(v); } - bool CheckBounds(const Matrix &M) const { return m_data.CheckBounds(M.m_data); } + + /// @brief Test whether a vector conforms to the matrix row shape + template bool ConformsToRow(const V &v) const { return m_data.ConformsToRow(v); } + /// @brief Test whether a vector conforms to the matrix column shape + template bool ConformsToColumn(const V &v) const { return m_data.ConformsToColumn(v); } + /// @brief Test whether another matrix conforms to the shape of this matrix + bool ConformsTo(const Matrix &M) const { return m_data.ConformsTo(M.m_data); } ///@} /// @name Comparison operators @@ -219,7 +224,7 @@ template Matrix &Matrix::operator=(const T &c) template bool Matrix::operator==(const Matrix &M) const { - if (!this->CheckBounds(M)) { + if (!this->ConformsTo(M)) { throw DimensionException(); } return std::equal(m_data.elements_begin(), m_data.elements_end(), M.m_data.elements_begin()); @@ -233,7 +238,7 @@ template bool Matrix::operator==(const T &c) const template Matrix &Matrix::operator+=(const Matrix &M) { - if (!this->CheckBounds(M)) { + if (!this->ConformsTo(M)) { throw DimensionException(); } std::transform(m_data.elements_begin(), m_data.elements_end(), M.m_data.elements_begin(), @@ -243,7 +248,7 @@ template Matrix &Matrix::operator+=(const Matrix &M) template Matrix &Matrix::operator-=(const Matrix &M) { - if (!this->CheckBounds(M)) { + if (!this->ConformsTo(M)) { throw DimensionException(); } std::transform(m_data.elements_begin(), m_data.elements_end(), M.m_data.elements_begin(), @@ -268,7 +273,7 @@ template Matrix &Matrix::operator*=(const T &c) template Matrix &Matrix::operator/=(const T &c) { - if (c == static_cast(0)) { + if (c == T{0}) { throw ZeroDivideException(); } std::transform(m_data.elements_begin(), m_data.elements_end(), m_data.elements_begin(), @@ -285,7 +290,7 @@ template template void Matrix::GetColumn(int col, V &v) c if (col < MinCol() || col > MaxCol()) { throw std::out_of_range("Index out of range in Matrix::GetColumn"); } - if (!CheckColumn(v)) { + if (!ConformsToColumn(v)) { throw DimensionException(); } for (int i = MinRow(); i <= MaxRow(); ++i) { @@ -298,7 +303,7 @@ template template void Matrix::SetColumn(int col, const V if (col < MinCol() || col > MaxCol()) { throw std::out_of_range("Index out of range in Matrix::SetColumn"); } - if (!CheckColumn(v)) { + if (!ConformsToColumn(v)) { throw DimensionException(); } for (int i = MinRow(); i <= MaxRow(); ++i) { @@ -311,7 +316,7 @@ template template void Matrix::GetRow(int row, V &v) cons if (row < MinRow() || row > MaxRow()) { throw std::out_of_range("Index out of range in Matrix::GetRow"); } - if (!CheckRow(v)) { + if (!ConformsToRow(v)) { throw DimensionException(); } for (int j = MinCol(); j <= MaxCol(); ++j) { @@ -324,7 +329,7 @@ template template void Matrix::SetRow(int row, const V &v if (row < MinRow() || row > MaxRow()) { throw std::out_of_range("Index out of range in Matrix::SetRow"); } - if (!CheckRow(v)) { + if (!ConformsToRow(v)) { throw DimensionException(); } for (int j = MinCol(); j <= MaxCol(); ++j) { @@ -338,7 +343,7 @@ template template void Matrix::SetRow(int row, const V &v template void Matrix::CMultiply(const Vector &p_input, Vector &p_output) const { - if (!this->CheckRow(p_input) || !this->CheckColumn(p_output)) { + if (!this->ConformsToRow(p_input) || !this->ConformsToColumn(p_output)) { throw DimensionException(); } for (int i = MinRow(); i <= MaxRow(); ++i) { @@ -349,7 +354,7 @@ template void Matrix::CMultiply(const Vector &p_input, Vector template void Matrix::RMultiply(const Vector &p_input, Vector &p_output) const { - if (!this->CheckColumn(p_input) || !this->CheckRow(p_output)) { + if (!this->ConformsToColumn(p_input) || !this->ConformsToRow(p_output)) { throw DimensionException(); } @@ -366,7 +371,7 @@ template void Matrix::RMultiply(const Vector &p_input, Vector template Vector Matrix::operator*(const Vector &v) const { - if (!this->CheckRow(v)) { + if (!this->ConformsToRow(v)) { throw DimensionException(); } Vector tmp(MinRow(), MaxRow()); @@ -575,7 +580,7 @@ template T Matrix::Determinant() const /// @sa Matrix::RMultiply template Vector operator*(const Vector &v, const Matrix &M) { - if (!M.CheckColumn(v)) { + if (!M.ConformsToColumn(v)) { throw DimensionException(); } Vector tmp(M.MinCol(), M.MaxCol()); diff --git a/src/core/recarray.h b/src/core/recarray.h index 2e6e3dd885..ada5c89edd 100644 --- a/src/core/recarray.h +++ b/src/core/recarray.h @@ -156,7 +156,7 @@ template class RectArray { /// @name Range checking functions; returns true only if valid index/size //@{ /// check array for same row and column boundaries - bool CheckBounds(const RectArray &m) const + bool ConformsTo(const RectArray &m) const { return (m_minrow == m.m_minrow && m_maxrow == m.m_maxrow && m_mincol == m.m_mincol && m_maxcol == m.m_maxcol); @@ -164,14 +164,14 @@ template class RectArray { /// check for correct row index bool CheckRow(const int row) const { return (m_minrow <= row && row <= m_maxrow); } /// check row vector for correct column boundaries - template bool CheckRow(const V &v) const + template bool ConformsToRow(const V &v) const { return v.front_index() == m_mincol && v.back_index() == m_maxcol; } /// check for correct column index bool CheckColumn(const int col) const { return (m_mincol <= col && col <= m_maxcol); } /// check column vector for correct row boundaries - template bool CheckColumn(const V &v) const + template bool ConformsToColumn(const V &v) const { return (v.front_index() == m_minrow && v.back_index() == m_maxrow); } @@ -292,7 +292,7 @@ template template void RectArray::GetRow(int row, Ve if (!CheckRow(row)) { throw std::out_of_range("Index out of range in RectArray"); } - if (!CheckRow(v)) { + if (!ConformsToRow(v)) { throw DimensionException(); } const size_t base = index(row, m_mincol); @@ -310,7 +310,7 @@ template template void RectArray::GetColumn(int col, if (!CheckColumn(col)) { throw std::out_of_range("Index out of range in RectArray"); } - if (!CheckColumn(v)) { + if (!ConformsToColumn(v)) { throw DimensionException(); } for (int r = m_minrow; r <= m_maxrow; ++r) { @@ -323,7 +323,7 @@ template template void RectArray::SetColumn(int col, if (!CheckColumn(col)) { throw std::out_of_range("Index out of range in RectArray"); } - if (!CheckColumn(v)) { + if (!ConformsToColumn(v)) { throw DimensionException(); } for (int r = m_minrow; r <= m_maxrow; ++r) {