diff --git a/src/Imath/ImathMatrix.h b/src/Imath/ImathMatrix.h
index a596030d..a4e2590d 100644
--- a/src/Imath/ImathMatrix.h
+++ b/src/Imath/ImathMatrix.h
@@ -292,6 +292,9 @@ template <class T> class IMATH_EXPORT_TEMPLATE_TYPE Matrix22
     /// Determinant
     IMATH_HOSTDEVICE constexpr T determinant () const IMATH_NOEXCEPT;
 
+    /// Trace
+    IMATH_HOSTDEVICE constexpr T trace() const IMATH_NOEXCEPT;
+
     /// Set matrix to rotation by r (in radians)
     /// @return const referenced to this
     template <class S>
@@ -678,6 +681,9 @@ template <class T> class IMATH_EXPORT_TEMPLATE_TYPE Matrix33
     /// Determinant
     IMATH_HOSTDEVICE constexpr T determinant () const IMATH_NOEXCEPT;
 
+    /// Trace
+    IMATH_HOSTDEVICE constexpr T trace() const IMATH_NOEXCEPT;
+
     /// Set matrix to rotation by r (in radians)
     /// @return const referenced to this
     template <class S>
@@ -1162,6 +1168,9 @@ template <class T> class IMATH_EXPORT_TEMPLATE_TYPE Matrix44
     /// Determinant
     IMATH_HOSTDEVICE IMATH_CONSTEXPR14 T determinant () const IMATH_NOEXCEPT;
 
+    /// Trace
+    IMATH_HOSTDEVICE constexpr T trace() const IMATH_NOEXCEPT;
+
     /// Set matrix to rotation by XYZ euler angles (in radians)
     /// @return const referenced to this
     template <class S>
@@ -1915,6 +1924,13 @@ Matrix22<T>::determinant () const IMATH_NOEXCEPT
     return x[0][0] * x[1][1] - x[1][0] * x[0][1];
 }
 
+template <class T>
+IMATH_HOSTDEVICE constexpr inline T
+Matrix22<T>::trace () const IMATH_NOEXCEPT
+{
+    return x[0][0] + x[1][1];
+}
+
 template <class T>
 template <class S>
 IMATH_HOSTDEVICE inline const Matrix22<T>&
@@ -3086,6 +3102,13 @@ Matrix33<T>::determinant () const IMATH_NOEXCEPT
            x[0][2] * (x[1][0] * x[2][1] - x[1][1] * x[2][0]);
 }
 
+template <class T>
+IMATH_HOSTDEVICE constexpr inline T
+Matrix33<T>::trace () const IMATH_NOEXCEPT
+{
+    return x[0][0] + x[1][1] + x[2][2];
+}
+
 template <class T>
 template <class S>
 IMATH_HOSTDEVICE inline const Matrix33<T>&
@@ -4592,6 +4615,13 @@ Matrix44<T>::determinant () const IMATH_NOEXCEPT
     return sum;
 }
 
+template <class T>
+IMATH_HOSTDEVICE constexpr inline T
+Matrix44<T>::trace () const IMATH_NOEXCEPT
+{
+    return x[0][0] + x[1][1] + x[2][2] + x[3][3];
+}
+
 template <class T>
 template <class S>
 IMATH_HOSTDEVICE inline const Matrix44<T>&
diff --git a/src/ImathTest/testMatrix.cpp b/src/ImathTest/testMatrix.cpp
index 8d37fee7..f248097c 100644
--- a/src/ImathTest/testMatrix.cpp
+++ b/src/ImathTest/testMatrix.cpp
@@ -640,6 +640,116 @@ testMatrix ()
             u.baseTypeEpsilon ());
     }
 
+    // Trace
+    {
+        cout << "2x2 trace" << endl;
+
+        IMATH_INTERNAL_NAMESPACE::Rand32 random;
+
+        IMATH_INTERNAL_NAMESPACE::M22f u;
+        float                          trace = 0;
+        for (int i = 0; i < 2; i++)
+        {
+            for (int j = 0; j < 2; j++)
+            {
+                const float randomNum = random.nextf ();
+                u[i][j]               = randomNum;
+                if (i == j) { trace += randomNum; }
+            }
+        }
+
+        assert (fabsf (u.trace () - trace) <= u.baseTypeEpsilon ());
+    }
+    {
+        IMATH_INTERNAL_NAMESPACE::Rand32 random;
+
+        IMATH_INTERNAL_NAMESPACE::M22d u;
+        double                         trace = 0;
+        for (int i = 0; i < 2; i++)
+        {
+            for (int j = 0; j < 2; j++)
+            {
+                const double randomNum = random.nextf ();
+                u[i][j]                = randomNum;
+                if (i == j) { trace += randomNum; }
+            }
+        }
+
+        assert (fabsf (u.trace () - trace) <= u.baseTypeEpsilon ());
+    }
+    {
+        cout << "3x3 trace" << endl;
+
+        IMATH_INTERNAL_NAMESPACE::Rand32 random;
+
+        IMATH_INTERNAL_NAMESPACE::M33f u;
+        float                          trace = 0;
+        for (int i = 0; i < 3; i++)
+        {
+            for (int j = 0; j < 3; j++)
+            {
+                const float randomNum = random.nextf ();
+                u[i][j]               = randomNum;
+                if (i == j) { trace += randomNum; }
+            }
+        }
+
+        assert (fabsf (u.trace () - trace) <= u.baseTypeEpsilon ());
+    }
+    {
+        IMATH_INTERNAL_NAMESPACE::Rand32 random;
+
+        IMATH_INTERNAL_NAMESPACE::M33d u;
+        double                         trace = 0;
+        for (int i = 0; i < 3; i++)
+        {
+            for (int j = 0; j < 3; j++)
+            {
+                const double randomNum = random.nextf ();
+                u[i][j]                = randomNum;
+                if (i == j) { trace += randomNum; }
+            }
+        }
+
+        assert (fabsf (u.trace () - trace) <= u.baseTypeEpsilon ());
+    }
+    {
+        cout << "4x4 trace" << endl;
+
+        IMATH_INTERNAL_NAMESPACE::Rand32 random;
+
+        IMATH_INTERNAL_NAMESPACE::M44f u;
+        float                          trace = 0;
+        for (int i = 0; i < 4; i++)
+        {
+            for (int j = 0; j < 4; j++)
+            {
+                const float randomNum = random.nextf ();
+                u[i][j]               = randomNum;
+                if (i == j) { trace += randomNum; }
+            }
+        }
+
+        assert (fabsf (u.trace () - trace) <= u.baseTypeEpsilon ());
+    }
+    {
+        IMATH_INTERNAL_NAMESPACE::Rand32 random;
+
+        IMATH_INTERNAL_NAMESPACE::M44d u;
+        double                         trace = 0;
+        for (int i = 0; i < 4; i++)
+        {
+            for (int j = 0; j < 4; j++)
+            {
+                const double randomNum = random.nextf ();
+                u[i][j]                = randomNum;
+                if (i == j) { trace += randomNum; }
+            }
+        }
+
+        assert (fabsf (u.trace () - trace) <= u.baseTypeEpsilon ());
+    }
+
     // Matrix minors
     {
         cout << "4x4 matrix minors" << endl;