diff --git a/README.md b/README.md
index 5e9c020364..c4363ab4c2 100644
--- a/README.md
+++ b/README.md
@@ -268,7 +268,7 @@ True
 
 Galois field arrays also contain matrix decomposition routines not included in NumPy. These include:
 
-- `FieldArray.row_reduce`, `FieldArray.lu_decompose`, `FieldArray.lup_decompose`
+- `FieldArray.row_reduce`, `FieldArray.lu_decompose`, `FieldArray.plu_decompose`
 
 #### NumPy ufunc methods
 
diff --git a/docs/basic-usage.rst b/docs/basic-usage.rst
index d3e79790c7..f466b40249 100644
--- a/docs/basic-usage.rst
+++ b/docs/basic-usage.rst
@@ -159,7 +159,7 @@ operation in :math:`\mathrm{GF}(p^m)` not in :math:`\mathbb{R}`. Some of these f
 
 Galois field arrays also contain matrix decomposition routines not included in NumPy. These include:
 
-- :func:`galois.FieldArray.row_reduce`, :func:`galois.FieldArray.lu_decompose`, :func:`galois.FieldArray.lup_decompose`
+- :func:`galois.FieldArray.row_reduce`, :func:`galois.FieldArray.lu_decompose`, :func:`galois.FieldArray.plu_decompose`
 
 NumPy ufunc methods
 ...................
diff --git a/galois/_fields/_linalg.py b/galois/_fields/_linalg.py
index b5abef6431..23ed908056 100644
--- a/galois/_fields/_linalg.py
+++ b/galois/_fields/_linalg.py
@@ -188,7 +188,7 @@ def lu_decompose(A):
     return L, U
 
 
-def lup_decompose(A):
+def plu_decompose(A):
     if not (A.ndim == 2 and A.shape[0] == A.shape[1]):
         raise ValueError(f"Argument `A` must be a square matrix, not {A.shape}.")
 
@@ -196,7 +196,7 @@ def lup_decompose(A):
     n = A.shape[0]
     Ai = A.copy()
     L = field.Zeros((n,n))
-    P = field.Identity(n)
+    P = field.Identity(n)  # Row permutation matrix
 
     for i in range(0, n-1):
         if Ai[i,i] == 0:
@@ -219,7 +219,8 @@ def lup_decompose(A):
     L[-1,-1] = 1  # Set the final diagonal to 1
     U = Ai
 
-    return L, U, P
+    # NOTE: Return column permutation matrix
+    return P.T, L, U
 
 
 ###############################################################################
@@ -271,7 +272,8 @@ def det(A):
     elif n == 3:
         return A[0,0]*(A[1,1]*A[2,2] - A[1,2]*A[2,1]) - A[0,1]*(A[1,0]*A[2,2] - A[1,2]*A[2,0]) + A[0,2]*(A[1,0]*A[2,1] - A[1,1]*A[2,0])
     else:
-        L, U, P = lup_decompose(A)
+        P, L, U = plu_decompose(A)
+        P = P.T  # Convert row permutation matrix into column permutation matrix
         idxs = np.arange(0, n)
         nrows = n - np.count_nonzero(P[idxs,idxs]) # The number of moved rows
         S = max(nrows - 1, 0)  # The number of row swaps
diff --git a/galois/_fields/_main.py b/galois/_fields/_main.py
index de32e18dc7..c4f74882cc 100644
--- a/galois/_fields/_main.py
+++ b/galois/_fields/_main.py
@@ -17,7 +17,7 @@
 from .._poly_conversion import integer_to_poly, poly_to_integer, str_to_integer, poly_to_str, sparse_poly_to_integer, sparse_poly_to_str, str_to_sparse_poly
 
 from ._dtypes import DTYPES
-from ._linalg import dot, row_reduce, lu_decompose, lup_decompose
+from ._linalg import dot, row_reduce, lu_decompose, plu_decompose
 from ._functions import FunctionMeta
 from ._ufuncs import UfuncMeta
 
@@ -1917,18 +1917,23 @@ def lu_decompose(self) -> "FieldArray":
         """
         return lu_decompose(self)
 
-    def lup_decompose(self) -> "FieldArray":
+    def plu_decompose(self) -> "FieldArray":
         r"""
         Decomposes the input array into the product of lower and upper triangular matrices using partial pivoting.
 
         Returns
         -------
+        galois.FieldArray
+            The column permutation matrix.
         galois.FieldArray
             The lower triangular matrix.
         galois.FieldArray
             The upper triangular matrix.
-        galois.FieldArray
-            The permutation matrix.
+
+        Notes
+        -----
+        The PLU decomposition of :math:`\mathbf{A}` is defined as :math:`\mathbf{A} = \mathbf{P} \mathbf{L} \mathbf{U}`. This is equivalent to
+        :math:`\mathbf{P}^T \mathbf{A} = \mathbf{L} \mathbf{U}`.
 
         Examples
         --------
@@ -1936,15 +1941,17 @@ def lup_decompose(self) -> "FieldArray":
 
             GF = galois.GF(5)
             A = GF([[1, 3, 2, 0], [3, 4, 2, 3], [0, 2, 1, 4], [4, 3, 3, 1]])
-            L, U, P = A.lup_decompose()
+            P, L, U = A.plu_decompose()
+            P
             L
             U
-            P
 
-            # P A = L U
-            np.array_equal(P @ A, L @ U)
+            # A = P L U
+            np.array_equal(A, P @ L @ U)
+            # P.T A = L U
+            np.array_equal(P.T @ A, L @ U)
         """
-        return lup_decompose(self)
+        return plu_decompose(self)
 
     def field_trace(self) -> "FieldArray":
         r"""
diff --git a/scripts/generate_field_test_vectors.py b/scripts/generate_field_test_vectors.py
index 3990c07c18..a525646736 100644
--- a/scripts/generate_field_test_vectors.py
+++ b/scripts/generate_field_test_vectors.py
@@ -385,6 +385,27 @@ def make_luts(field, sub_folder, seed, sparse=False):
     d = {"X": X, "L": L, "U": U}
     save_pickle(d, folder, "lu_decompose.pkl")
 
+    set_seed(seed + 204)
+    shapes = [(2,2), (3,3), (4,4), (5,5), (6,6)]
+    X = []
+    L = []
+    U = []
+    P = []
+    for i in range(len(shapes)):
+        x = randint_matrix(0, order, shapes[i])
+        dtype = x.dtype
+        X.append(x)
+        x = matrix(FIELD, [[F(e) for e in row] for row in x])
+        p, l, u = x.LU()
+        p = np.array([[I(e) for e in row] for row in p], dtype)
+        l = np.array([[I(e) for e in row] for row in l], dtype)
+        u = np.array([[I(e) for e in row] for row in u], dtype)
+        P.append(p)
+        L.append(l)
+        U.append(u)
+    d = {"X": X, "P": P, "L": L, "U": U}
+    save_pickle(d, folder, "plu_decompose.pkl")
+
     ###############################################################################
     # Polynomial arithmetic
     ###############################################################################
diff --git a/tests/fields/conftest.py b/tests/fields/conftest.py
index eefdb06121..c31a432490 100644
--- a/tests/fields/conftest.py
+++ b/tests/fields/conftest.py
@@ -307,6 +307,20 @@ def field_lu_decompose(field_folder):
     return d
 
 
+@pytest.fixture(scope="session")
+def field_plu_decompose(field_folder):
+    GF, folder = field_folder
+    with open(os.path.join(folder, "plu_decompose.pkl"), "rb") as f:
+        print(f"Loading {f}...")
+        d = pickle.load(f)
+    d["GF"] = GF
+    d["X"] = [GF(x) for x in d["X"]]
+    d["P"] = [GF(p) for p in d["P"]]
+    d["L"] = [GF(l) for l in d["L"]]
+    d["U"] = [GF(u) for u in d["U"]]
+    return d
+
+
 ###############################################################################
 # Fixtures for arithmetic methods over finite fields
 ###############################################################################
diff --git a/tests/fields/data/GF(109987^4)/plu_decompose.pkl b/tests/fields/data/GF(109987^4)/plu_decompose.pkl
new file mode 100644
index 0000000000..02fb043a7d
Binary files /dev/null and b/tests/fields/data/GF(109987^4)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(2)/plu_decompose.pkl b/tests/fields/data/GF(2)/plu_decompose.pkl
new file mode 100644
index 0000000000..f2217c8a4b
Binary files /dev/null and b/tests/fields/data/GF(2)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(2147483647)/plu_decompose.pkl b/tests/fields/data/GF(2147483647)/plu_decompose.pkl
new file mode 100644
index 0000000000..5c98455c31
Binary files /dev/null and b/tests/fields/data/GF(2147483647)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(2^100)/plu_decompose.pkl b/tests/fields/data/GF(2^100)/plu_decompose.pkl
new file mode 100644
index 0000000000..11a4f5d06c
Binary files /dev/null and b/tests/fields/data/GF(2^100)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(2^2)/plu_decompose.pkl b/tests/fields/data/GF(2^2)/plu_decompose.pkl
new file mode 100644
index 0000000000..4ef2469fe4
Binary files /dev/null and b/tests/fields/data/GF(2^2)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(2^3)/plu_decompose.pkl b/tests/fields/data/GF(2^3)/plu_decompose.pkl
new file mode 100644
index 0000000000..da5603ad6b
Binary files /dev/null and b/tests/fields/data/GF(2^3)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(2^32)/plu_decompose.pkl b/tests/fields/data/GF(2^32)/plu_decompose.pkl
new file mode 100644
index 0000000000..f8676c3584
Binary files /dev/null and b/tests/fields/data/GF(2^32)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(2^8)/plu_decompose.pkl b/tests/fields/data/GF(2^8)/plu_decompose.pkl
new file mode 100644
index 0000000000..a6b8e5ce38
Binary files /dev/null and b/tests/fields/data/GF(2^8)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(2^8, 283, 19)/plu_decompose.pkl b/tests/fields/data/GF(2^8, 283, 19)/plu_decompose.pkl
new file mode 100644
index 0000000000..c022b303ba
Binary files /dev/null and b/tests/fields/data/GF(2^8, 283, 19)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(31)/plu_decompose.pkl b/tests/fields/data/GF(31)/plu_decompose.pkl
new file mode 100644
index 0000000000..a7af0f067a
Binary files /dev/null and b/tests/fields/data/GF(31)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(3191)/plu_decompose.pkl b/tests/fields/data/GF(3191)/plu_decompose.pkl
new file mode 100644
index 0000000000..fd8dd1cc27
Binary files /dev/null and b/tests/fields/data/GF(3191)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(36893488147419103183)/plu_decompose.pkl b/tests/fields/data/GF(36893488147419103183)/plu_decompose.pkl
new file mode 100644
index 0000000000..378a391416
Binary files /dev/null and b/tests/fields/data/GF(36893488147419103183)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(5)/plu_decompose.pkl b/tests/fields/data/GF(5)/plu_decompose.pkl
new file mode 100644
index 0000000000..8a49c80cb3
Binary files /dev/null and b/tests/fields/data/GF(5)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(7)/plu_decompose.pkl b/tests/fields/data/GF(7)/plu_decompose.pkl
new file mode 100644
index 0000000000..eb6b49a999
Binary files /dev/null and b/tests/fields/data/GF(7)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(7^3)/plu_decompose.pkl b/tests/fields/data/GF(7^3)/plu_decompose.pkl
new file mode 100644
index 0000000000..0217c562d6
Binary files /dev/null and b/tests/fields/data/GF(7^3)/plu_decompose.pkl differ
diff --git a/tests/fields/data/GF(7^3, 643, 244)/plu_decompose.pkl b/tests/fields/data/GF(7^3, 643, 244)/plu_decompose.pkl
new file mode 100644
index 0000000000..09fb688fa2
Binary files /dev/null and b/tests/fields/data/GF(7^3, 643, 244)/plu_decompose.pkl differ
diff --git a/tests/fields/test_linalg.py b/tests/fields/test_linalg.py
index f7748d5664..645eab6a10 100644
--- a/tests/fields/test_linalg.py
+++ b/tests/fields/test_linalg.py
@@ -223,15 +223,6 @@ def test_matmul_2d_2d(field):
 #     assert array_equal(A @ B, np.matmul(A, B))
 
 
-def test_lup_decomposition():
-    GF = galois.GF(3)
-    I = GF.Identity(3)
-    for trial in range(100):
-        A = GF.Random((3,3))
-        L, U, P = A.lup_decompose()
-        assert array_equal(L @ U, P @ A)
-
-
 def test_det_2x2(field):
     dtype = random.choice(field.dtypes)
     a = field.Random(dtype=dtype)
@@ -342,3 +333,18 @@ def test_lu_decompose(field_lu_decompose):
         assert np.array_equal(u, U[i])
         assert type(l) is GF
         assert type(u) is GF
+
+
+def test_plu_decompose(field_plu_decompose):
+    GF, X, P, L, U = field_plu_decompose["GF"], field_plu_decompose["X"], field_plu_decompose["P"], field_plu_decompose["L"], field_plu_decompose["U"]
+
+    for i in range(len(X)):
+        dtype = random.choice(GF.dtypes)
+        x = X[i].astype(dtype)
+        p, l, u = x.plu_decompose()
+        assert np.array_equal(p, P[i])
+        assert np.array_equal(l, L[i])
+        assert np.array_equal(u, U[i])
+        assert type(p) is GF
+        assert type(l) is GF
+        assert type(u) is GF