diff --git a/Makefile.am b/Makefile.am
index 654344066..276d533b0 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -251,14 +251,11 @@ core_SOURCES = \
 	src/core/recarray.h \
 	src/core/matrix.h \
 	src/core/matrix.imp \
-	src/core/sqmatrix.h \
-	src/core/sqmatrix.imp \
 	src/core/integer.cc \
 	src/core/integer.h \
 	src/core/rational.cc \
 	src/core/rational.h \
 	src/core/matrix.cc \
-	src/core/sqmatrix.cc \
 	src/core/function.cc \
 	src/core/function.h \
 	src/core/tinyxml.cc \
diff --git a/src/core/matrix.h b/src/core/matrix.h
index fecad1d55..04920da69 100644
--- a/src/core/matrix.h
+++ b/src/core/matrix.h
@@ -28,6 +28,12 @@
 
 namespace Gambit {
 
+class SingularMatrixException final : public std::runtime_error {
+public:
+  SingularMatrixException() : std::runtime_error("Attempted to invert a singular matrix") {}
+  ~SingularMatrixException() noexcept override = default;
+};
+
 template <class T> Vector<T> operator*(const Vector<T> &, const Matrix<T> &);
 
 template <class T> class Matrix : public RectArray<T> {
@@ -49,6 +55,10 @@ template <class T> class Matrix : public RectArray<T> {
 
   /// @name Extracting rows and columns
   //@{
+  bool IsSquare() const
+  {
+    return this->MinRow() == this->MinCol() && this->MaxRow() == this->MaxCol();
+  }
   Vector<T> Row(int) const;
   Vector<T> Column(int) const;
   //@}
@@ -93,6 +103,9 @@ template <class T> class Matrix : public RectArray<T> {
   void MakeIdent();
   void Pivot(int, int);
   //@}
+
+  Matrix Inverse() const;
+  T Determinant() const;
 };
 
 template <class T> Vector<T> operator*(const Vector<T> &, const Matrix<T> &);
diff --git a/src/core/matrix.imp b/src/core/matrix.imp
index ae87aba82..5650583fd 100644
--- a/src/core/matrix.imp
+++ b/src/core/matrix.imp
@@ -384,6 +384,9 @@ template <class T> Vector<T> Matrix<T>::Column(int j) const
 
 template <class T> void Matrix<T>::MakeIdent()
 {
+  if (!IsSquare()) {
+    throw DimensionException();
+  }
   for (int i = this->minrow; i <= this->maxrow; i++) {
     for (int j = this->mincol; j <= this->maxcol; j++) {
       if (i == j) {
@@ -426,4 +429,139 @@ template <class T> void Matrix<T>::Pivot(int row, int col)
   }
 }
 
+template <class T> Matrix<T> Matrix<T>::Inverse() const
+{
+  if (!IsSquare()) {
+    throw DimensionException();
+  }
+
+  Matrix copy(*this);
+  Matrix inv(this->MaxRow(), this->MaxRow());
+
+  // initialize inverse matrix and prescale row vectors
+  for (int i = this->MinRow(); i <= this->MaxRow(); i++) {
+    T max = (T)0;
+    for (int j = this->MinCol(); j <= this->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 = this->MinCol(); j <= this->MaxCol(); j++) {
+      copy(i, j) *= scale;
+      if (i == j) {
+        inv(i, j) = scale;
+      }
+      else {
+        inv(i, j) = (T)0;
+      }
+    }
+  }
+
+  for (int i = this->MinCol(); i <= this->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++) {
+      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 = this->MinCol(); k <= this->MaxCol(); k++) {
+      copy(i, k) *= factor;
+      inv(i, k) *= factor;
+    }
+
+    // reduce other rows
+    for (int j = this->MinRow(); j <= this->MaxRow(); j++) {
+      if (j != i) {
+        T mult = copy(j, i);
+        for (int k = this->MinCol(); k <= this->MaxCol(); k++) {
+          copy(j, k) -= copy(i, k) * mult;
+          inv(j, k) -= inv(i, k) * mult;
+        }
+      }
+    }
+  }
+
+  return inv;
+}
+
+template <class T> T Matrix<T>::Determinant() const
+{
+  if (!IsSquare()) {
+    throw DimensionException();
+  }
+
+  T factor = (T)1;
+  Matrix M(*this);
+
+  for (int row = this->MinRow(); row <= this->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])) {
+        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;
+      }
+    }
+
+    if (M.data[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;
+      }
+    }
+  }
+
+  // finally we multiply the diagonal elements
+  T det = (T)1;
+  for (int row = this->MinRow(); row <= this->MaxRow(); row++) {
+    det *= M.data[row][row];
+  }
+  return det;
+}
+
 } // end namespace Gambit
diff --git a/src/core/sqmatrix.cc b/src/core/sqmatrix.cc
deleted file mode 100644
index 20f0295c9..000000000
--- a/src/core/sqmatrix.cc
+++ /dev/null
@@ -1,27 +0,0 @@
-//
-// This file is part of Gambit
-// Copyright (c) 1994-2025, The Gambit Project (https://www.gambit-project.org)
-//
-// FILE: src/core/sqmatrix.cc
-// Instantiation of common square matrix types
-//
-// 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 "sqmatrix.imp"
-#include "rational.h"
-
-template class Gambit::SquareMatrix<double>;
-template class Gambit::SquareMatrix<Gambit::Rational>;
diff --git a/src/core/sqmatrix.h b/src/core/sqmatrix.h
deleted file mode 100644
index 0a9c5152e..000000000
--- a/src/core/sqmatrix.h
+++ /dev/null
@@ -1,52 +0,0 @@
-//
-// This file is part of Gambit
-// Copyright (c) 1994-2025, The Gambit Project (https://www.gambit-project.org)
-//
-// FILE: src/core/sqmatrix.h
-// Implementation of square matrices
-//
-// 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.
-//
-
-#ifndef GAMBIT_CORE_SQMATRIX_H
-#define GAMBIT_CORE_SQMATRIX_H
-
-#include "matrix.h"
-
-namespace Gambit {
-
-class SingularMatrixException final : public std::runtime_error {
-public:
-  SingularMatrixException() : std::runtime_error("Attempted to invert a singular matrix") {}
-  ~SingularMatrixException() noexcept override = default;
-};
-
-template <class T> class SquareMatrix : public Matrix<T> {
-public:
-  SquareMatrix() = default;
-  explicit SquareMatrix(int size);
-  explicit SquareMatrix(const Matrix<T> &);
-  SquareMatrix(const SquareMatrix<T> &);
-  ~SquareMatrix() override = default;
-
-  SquareMatrix<T> &operator=(const SquareMatrix<T> &);
-
-  SquareMatrix<T> Inverse() const;
-  T Determinant() const;
-};
-
-} // end namespace Gambit
-
-#endif // GAMBIT_CORE_SQMATRIX_H
diff --git a/src/core/sqmatrix.imp b/src/core/sqmatrix.imp
deleted file mode 100644
index e9109b70a..000000000
--- a/src/core/sqmatrix.imp
+++ /dev/null
@@ -1,178 +0,0 @@
-//
-// This file is part of Gambit
-// Copyright (c) 1994-2025, The Gambit Project (https://www.gambit-project.org)
-//
-// FILE: src/core/sqmatrix.imp
-// Implementation of square matrices
-//
-// 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 "sqmatrix.h"
-
-namespace Gambit {
-
-//-----------------------------------------------------------------------------
-//     SquareMatrix<T>: Constructors, destructor, constructive operators
-//-----------------------------------------------------------------------------
-
-template <class T> SquareMatrix<T>::SquareMatrix(int size) : Matrix<T>(size, size) {}
-
-template <class T> SquareMatrix<T>::SquareMatrix(const Matrix<T> &M) : Matrix<T>(M) {}
-
-template <class T> SquareMatrix<T>::SquareMatrix(const SquareMatrix<T> &M) : Matrix<T>(M) {}
-
-template <class T> SquareMatrix<T> &SquareMatrix<T>::operator=(const SquareMatrix<T> &M)
-{
-  Matrix<T>::operator=(M);
-  return *this;
-}
-
-//-----------------------------------------------------------------------------
-//                SquareMatrix<T>: Public member functions
-//-----------------------------------------------------------------------------
-
-template <class T> SquareMatrix<T> SquareMatrix<T>::Inverse() const
-{
-  if (this->mincol != this->minrow || this->maxcol != this->maxrow) {
-    throw DimensionException();
-  }
-
-  SquareMatrix<T> copy(*this);
-  SquareMatrix<T> inv(this->maxrow);
-
-  // initialize inverse matrix and prescale row vectors
-  for (int i = this->minrow; i <= this->maxrow; i++) {
-    T max = (T)0;
-    for (int j = this->mincol; j <= this->maxcol; j++) {
-      T abs = copy.data[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 = this->mincol; j <= this->maxcol; j++) {
-      copy.data[i][j] *= scale;
-      if (i == j) {
-        inv.data[i][j] = scale;
-      }
-      else {
-        inv.data[i][j] = (T)0;
-      }
-    }
-  }
-
-  for (int i = this->mincol; i <= this->maxcol; i++) {
-    // find pivot row
-    T max = copy.data[i][i];
-    if (max < (T)0) {
-      max = -max;
-    }
-    int row = i;
-    for (int j = i + 1; j <= this->maxrow; j++) {
-      T abs = copy.data[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.data[i][i];
-    for (int k = this->mincol; k <= this->maxcol; k++) {
-      copy.data[i][k] *= factor;
-      inv.data[i][k] *= factor;
-    }
-
-    // reduce other rows
-    for (int j = this->minrow; j <= this->maxrow; j++) {
-      if (j != i) {
-        T mult = copy.data[j][i];
-        for (int k = this->mincol; k <= this->maxcol; k++) {
-          copy.data[j][k] -= copy.data[i][k] * mult;
-          inv.data[j][k] -= inv.data[i][k] * mult;
-        }
-      }
-    }
-  }
-
-  return inv;
-}
-
-template <class T> T SquareMatrix<T>::Determinant() const
-{
-  T factor = (T)1;
-  SquareMatrix<T> M(*this);
-
-  for (int row = this->minrow; row <= this->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])) {
-        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;
-      }
-    }
-
-    if (M.data[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;
-      }
-    }
-  }
-
-  // finally we multiply the diagonal elements
-  T det = (T)1;
-  for (int row = this->minrow; row <= this->maxrow; row++) {
-    det *= M.data[row][row];
-  }
-  return det;
-}
-
-} // end namespace Gambit
diff --git a/src/solvers/enumpoly/poly.h b/src/solvers/enumpoly/poly.h
index 6e0830ee6..4bd2550e7 100644
--- a/src/solvers/enumpoly/poly.h
+++ b/src/solvers/enumpoly/poly.h
@@ -24,7 +24,6 @@
 #define POLY_H
 
 #include "gambit.h"
-#include "core/sqmatrix.h"
 
 namespace Gambit {
 
@@ -193,7 +192,7 @@ template <class T> class Polynomial {
   Polynomial<T> DivideByPolynomial(const Polynomial<T> &den) const;
 
   Polynomial<T> TranslateOfMono(const Monomial<T> &, const Vector<T> &) const;
-  Polynomial<T> MonoInNewCoordinates(const Monomial<T> &, const SquareMatrix<T> &) const;
+  Polynomial<T> MonoInNewCoordinates(const Monomial<T> &, const Matrix<T> &) const;
 
 public:
   Polynomial(std::shared_ptr<VariableSpace> p) : m_space(p) {}
@@ -329,7 +328,7 @@ template <class T> class Polynomial {
   const List<Monomial<T>> &MonomialList() const { return m_terms; }
 
   Polynomial<T> TranslateOfPoly(const Vector<T> &) const;
-  Polynomial<T> PolyInNewCoordinates(const SquareMatrix<T> &) const;
+  Polynomial<T> PolyInNewCoordinates(const Matrix<T> &) const;
   T MaximalValueOfNonlinearPart(const T &) const;
   Polynomial<T> Normalize() const;
 };
diff --git a/src/solvers/enumpoly/poly.imp b/src/solvers/enumpoly/poly.imp
index 854189da4..9be7d8feb 100644
--- a/src/solvers/enumpoly/poly.imp
+++ b/src/solvers/enumpoly/poly.imp
@@ -238,8 +238,7 @@ template <class T> Polynomial<T> Polynomial<T>::TranslateOfPoly(const Vector<T>
 }
 
 template <class T>
-Polynomial<T> Polynomial<T>::MonoInNewCoordinates(const Monomial<T> &m,
-                                                  const SquareMatrix<T> &M) const
+Polynomial<T> Polynomial<T>::MonoInNewCoordinates(const Monomial<T> &m, const Matrix<T> &M) const
 {
   Polynomial<T> answer(m_space, static_cast<T>(1));
 
@@ -258,8 +257,7 @@ Polynomial<T> Polynomial<T>::MonoInNewCoordinates(const Monomial<T> &m,
   return answer;
 }
 
-template <class T>
-Polynomial<T> Polynomial<T>::PolyInNewCoordinates(const SquareMatrix<T> &M) const
+template <class T> Polynomial<T> Polynomial<T>::PolyInNewCoordinates(const Matrix<T> &M) const
 {
   Polynomial<T> answer(m_space);
   for (const auto &term : m_terms) {
diff --git a/src/solvers/enumpoly/polypartial.h b/src/solvers/enumpoly/polypartial.h
index d7084bd46..4a2eca425 100644
--- a/src/solvers/enumpoly/polypartial.h
+++ b/src/solvers/enumpoly/polypartial.h
@@ -98,7 +98,7 @@ template <class T> class PolynomialSystemDerivatives {
   }
   int GetDimension() const { return m_system.front().GetDimension(); }
   Matrix<T> DerivativeMatrix(const Vector<T> &, int) const;
-  SquareMatrix<T> SquareDerivativeMatrix(const Vector<T> &) const;
+  Matrix<T> SquareDerivativeMatrix(const Vector<T> &) const;
   Vector<T> ValuesOfRootPolys(const Vector<T> &, int) const;
 };
 
diff --git a/src/solvers/enumpoly/polypartial.imp b/src/solvers/enumpoly/polypartial.imp
index 9248c5251..41854f6fd 100644
--- a/src/solvers/enumpoly/polypartial.imp
+++ b/src/solvers/enumpoly/polypartial.imp
@@ -172,9 +172,9 @@ Matrix<T> PolynomialSystemDerivatives<T>::DerivativeMatrix(const Vector<T> &p,
 }
 
 template <class T>
-SquareMatrix<T> PolynomialSystemDerivatives<T>::SquareDerivativeMatrix(const Vector<T> &p) const
+Matrix<T> PolynomialSystemDerivatives<T>::SquareDerivativeMatrix(const Vector<T> &p) const
 {
-  SquareMatrix<T> answer(GetDimension());
+  Matrix<T> answer(GetDimension(), GetDimension());
   auto poly = m_system.begin();
   for (int i = 1; i <= GetDimension(); i++, poly++) {
     for (int j = 1; j <= GetDimension(); j++) {
diff --git a/src/solvers/enumpoly/polysolver.cc b/src/solvers/enumpoly/polysolver.cc
index 7a59fd5c9..5268ae1b2 100644
--- a/src/solvers/enumpoly/polysolver.cc
+++ b/src/solvers/enumpoly/polysolver.cc
@@ -110,7 +110,7 @@ bool PolynomialSystemSolver::NewtonRootInRectangle(const Rectangle<double> &r,
 }
 
 double PolynomialSystemSolver::MaxDistanceFromPointToVertexAfterTransformation(
-    const Rectangle<double> &r, const Vector<double> &p, const SquareMatrix<double> &M) const
+    const Rectangle<double> &r, const Vector<double> &p, const Matrix<double> &M) const
 {
   // A very early implementation of this method used a type gDouble which
   // implemented fuzzy comparisons.  Adding the epsilon parameter here is
@@ -152,7 +152,7 @@ Vector<double> PolynomialSystemSolver::NewtonStep(const Vector<double> &point) c
   const Vector<double> evals = m_derivatives.ValuesOfRootPolys(point, NumEquations());
   const Matrix<double> deriv = m_derivatives.DerivativeMatrix(point, NumEquations());
   auto transpose = deriv.Transpose();
-  auto sqmat = SquareMatrix<double>(deriv * transpose);
+  auto sqmat = Matrix<double>(deriv * transpose);
   if (std::abs(sqmat.Determinant()) <= 1.0e-9) {
     throw SingularMatrixException();
   }
@@ -164,7 +164,7 @@ Vector<double> PolynomialSystemSolver::ImprovingNewtonStep(const Vector<double>
   const Vector<double> evals = m_derivatives.ValuesOfRootPolys(point, NumEquations());
   const Matrix<double> deriv = m_derivatives.DerivativeMatrix(point, NumEquations());
   auto transpose = deriv.Transpose();
-  auto sqmat = SquareMatrix<double>(deriv * transpose);
+  auto sqmat = Matrix<double>(deriv * transpose);
   if (std::abs(sqmat.Determinant()) <= 1.0e-9) {
     throw SingularMatrixException();
   }
diff --git a/src/solvers/enumpoly/polysolver.h b/src/solvers/enumpoly/polysolver.h
index 7e9d5ecb1..38915a2a1 100644
--- a/src/solvers/enumpoly/polysolver.h
+++ b/src/solvers/enumpoly/polysolver.h
@@ -66,7 +66,7 @@ class PolynomialSystemSolver {
   // the one produced by the last step
   double MaxDistanceFromPointToVertexAfterTransformation(const Rectangle<double> &,
                                                          const Vector<double> &,
-                                                         const SquareMatrix<double> &) const;
+                                                         const Matrix<double> &) const;
 
   bool HasNoOtherRootsIn(const Rectangle<double> &, const Vector<double> &) const;
 
diff --git a/src/solvers/enumpoly/polysystem.h b/src/solvers/enumpoly/polysystem.h
index 98a469f71..919b1fafc 100644
--- a/src/solvers/enumpoly/polysystem.h
+++ b/src/solvers/enumpoly/polysystem.h
@@ -89,7 +89,7 @@ template <class T> class PolynomialSystem {
     return ret;
   }
   /// Transform system to new coordinates
-  PolynomialSystem<T> TransformCoords(const SquareMatrix<T> &M) const
+  PolynomialSystem<T> TransformCoords(const Matrix<T> &M) const
   {
     PolynomialSystem<T> ret(m_space);
     for (const auto &v : m_system) {
diff --git a/src/solvers/logit/path.cc b/src/solvers/logit/path.cc
index 00743be76..2bfc2ef04 100644
--- a/src/solvers/logit/path.cc
+++ b/src/solvers/logit/path.cc
@@ -24,7 +24,6 @@
 #include <algorithm> // for std::max
 
 #include "gambit.h"
-#include "core/sqmatrix.h"
 #include "path.h"
 
 namespace Gambit {
@@ -143,7 +142,7 @@ void PathTracer::TracePath(
   Vector<double> t(x.size()), newT(x.size());
   Vector<double> y(x.size() - 1);
   Matrix<double> b(x.size(), x.size() - 1);
-  SquareMatrix<double> q(x.size());
+  Matrix<double> q(x.size(), x.size());
 
   p_jacobian(x, b);
   QRDecomp(b, q);