diff --git a/.gitignore b/.gitignore
index 303ac4d..11328a7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,4 +3,5 @@
 *.vs*
 */bin/*
 */obj/*
-*/packages/*
\ No newline at end of file
+*/packages/*
+GenericTensor.sln.DotSettings.user
\ No newline at end of file
diff --git a/GenericTensor.sln.DotSettings.user b/GenericTensor.sln.DotSettings.user
deleted file mode 100644
index e0814be..0000000
--- a/GenericTensor.sln.DotSettings.user
+++ /dev/null
@@ -1,2 +0,0 @@
-<wpf:ResourceDictionary xml:space="preserve" xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml" xmlns:s="clr-namespace:System;assembly=mscorlib" xmlns:ss="urn:shemas-jetbrains-com:settings-storage-xaml" xmlns:wpf="http://schemas.microsoft.com/winfx/2006/xaml/presentation">
-	<s:String x:Key="/Default/Profiling/Configurations/=0/@EntryIndexedValue">&lt;data&gt;&lt;HostParameters type="LocalHostParameters" /&gt;&lt;Argument type="AttachArgument"&gt;&lt;DisplayName&gt;D:\main\vs_prj\GenericTensor\GenericTensor\Sample\bin\Release\netcoreapp3.1\Sample.exe&lt;/DisplayName&gt;&lt;ProcessId&gt;12684&lt;/ProcessId&gt;&lt;RuntimeType&gt;CoreClr&lt;/RuntimeType&gt;&lt;/Argument&gt;&lt;Info type="TimelineInfo" /&gt;&lt;CoreOptions type="CoreOptions"&gt;&lt;CoreTempPath IsNull="False"&gt;&lt;/CoreTempPath&gt;&lt;RemoteEndPoint IsNull="False"&gt;&lt;/RemoteEndPoint&gt;&lt;AdditionalEnvironmentVariables /&gt;&lt;/CoreOptions&gt;&lt;HostOptions type="HostOptions"&gt;&lt;HostTempPath IsNull="False"&gt;&lt;/HostTempPath&gt;&lt;/HostOptions&gt;&lt;/data&gt;</s:String></wpf:ResourceDictionary>
\ No newline at end of file
diff --git a/GenericTensor/Core/Exceptions.cs b/GenericTensor/Core/Exceptions.cs
index 17d8563..ab4d4d3 100644
--- a/GenericTensor/Core/Exceptions.cs
+++ b/GenericTensor/Core/Exceptions.cs
@@ -55,4 +55,13 @@ public class InvalidDeterminantException : DataException
         internal InvalidDeterminantException(string msg) : base(msg) {}
         internal InvalidDeterminantException() : base() {}
     }
+    
+    /// <summary>
+    /// Thrown when there is no decomposition for provided matrix
+    /// </summary>
+    public class ImpossibleDecomposition : DataException
+    {
+        internal ImpossibleDecomposition(string msg) : base(msg) {}
+        internal ImpossibleDecomposition() : base() {}
+    }
 }
diff --git a/GenericTensor/Declaration.cs b/GenericTensor/Declaration.cs
index 5dec790..bd98273 100644
--- a/GenericTensor/Declaration.cs
+++ b/GenericTensor/Declaration.cs
@@ -1,4 +1,5 @@
 #region copyright
+
 /*
  * MIT License
  * 
@@ -22,6 +23,7 @@
  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  * SOFTWARE.
  */
+
 #endregion
 
 
@@ -40,6 +42,7 @@ public sealed partial class GenTensor<T, TWrapper>
         : ICloneable where TWrapper : struct, IOperations<T>
     {
         #region Composition
+
         /// <summary>
         /// Creates a new axis that is put backward
         /// and then sets all elements as children
@@ -67,6 +70,7 @@ public static GenTensor<T, TWrapper> Concat(GenTensor<T, TWrapper> a, GenTensor<
         #endregion
 
         #region Constructors
+
         /// <summary>
         /// Creates a tensor whose all matrices are identity matrices
         /// <para>1 is achieved with <see cref="IOperations{T}.CreateOne"/></para>
@@ -145,36 +149,31 @@ public static GenTensor<T, TWrapper> CreateMatrix(int width, int height)
         /// (its only argument is an array of integers which are indices of the tensor)
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> CreateTensor(TensorShape shape, Func<int[], T> operation, Threading threading = Threading.Single)
-            => Constructors<T, TWrapper>.CreateTensor(shape, operation, threading);
+        public static GenTensor<T, TWrapper> CreateTensor(TensorShape shape, Func<int[], T> operation, Threading threading = Threading.Single) => Constructors<T, TWrapper>.CreateTensor(shape, operation, threading);
 
         /// <summary>
         /// Creates a tensor from an array
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> CreateTensor(T[] data)
-            => Constructors<T, TWrapper>.CreateTensor(data);
+        public static GenTensor<T, TWrapper> CreateTensor(T[] data) => Constructors<T, TWrapper>.CreateTensor(data);
 
         /// <summary>
         /// Creates a tensor from a two-dimensional array
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> CreateTensor(T[,] data)
-            => Constructors<T, TWrapper>.CreateTensor(data);
+        public static GenTensor<T, TWrapper> CreateTensor(T[,] data) => Constructors<T, TWrapper>.CreateTensor(data);
 
         /// <summary>
         /// Creates a tensor from a three-dimensional array
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> CreateTensor(T[,,] data)
-            => Constructors<T, TWrapper>.CreateTensor(data);
+        public static GenTensor<T, TWrapper> CreateTensor(T[,,] data) => Constructors<T, TWrapper>.CreateTensor(data);
 
         /// <summary>
         /// Creates a tensor from an n-dimensional array
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> CreateTensor(Array data)
-            => Constructors<T, TWrapper>.CreateTensor(data);
+        public static GenTensor<T, TWrapper> CreateTensor(Array data) => Constructors<T, TWrapper>.CreateTensor(data);
 
         #endregion
 
@@ -186,6 +185,42 @@ public static GenTensor<T, TWrapper> CreateTensor(Array data)
         /// </summary>
         public GenTensor<T, TWrapper> RowEchelonFormSimple()
             => EchelonForm<T, TWrapper>.RowEchelonFormSimple(this);
+        
+        /// <summary>
+        /// Decomposes a matrix into a triangular one.
+        /// Is of the Row Echelon Form (leading elements might be differ from ones).
+        /// 
+        /// In addition, returns a permutation that algorithm performs with rows.
+        /// Permutation array is size of rows there are in matrix.
+        ///
+        /// Initial state of that array is:
+        /// 1 2 3 ... numberOfRows
+        /// 
+        /// For example, algorithm swaps first and third rows then:
+        /// 3 2 1 ... numberOfRows
+        ///
+        /// It can be useful performing decompositions
+        /// </summary>
+        public (GenTensor<T, TWrapper>, int[]) RowEchelonFormPermuteSimple()
+            => EchelonForm<T, TWrapper>.RowEchelonFormPermuteSimple(this);
+        
+        /// <summary>
+        /// Decomposes a matrix into a triangular one.
+        /// Is of the Row Echelon Form (leading elements might be differ from ones).
+        /// 
+        /// In addition, returns a permutation that algorithm performs with rows.
+        /// Permutation array is size of rows there are in matrix.
+        ///
+        /// Initial state of that array is:
+        /// 1 2 3 ... numberOfRows
+        /// 
+        /// For example, algorithm swaps first and third rows then:
+        /// 3 2 1 ... numberOfRows
+        ///
+        /// It can be useful performing decompositions
+        /// </summary>
+        public (GenTensor<T, TWrapper>, int[]) RowEchelonFormPermuteSafeDivision()
+            => EchelonForm<T, TWrapper>.RowEchelonFormPermuteSafeDivision(this);
 
         /// <summary>
         /// Decomposes a matrix into a triangular one.
@@ -196,7 +231,6 @@ public GenTensor<T, TWrapper> RowEchelonFormSafeDivision()
             => EchelonForm<T, TWrapper>.RowEchelonFormSafeDivision(this);
 
 
-
         /// <summary>
         /// Decomposes a matrix into a triangular one.
         /// Is of the Row Echelon Form (leading elements are ones).
@@ -215,7 +249,6 @@ public GenTensor<T, TWrapper> RowEchelonFormLeadingOnesSafeDivision()
             => EchelonForm<T, TWrapper>.RowEchelonFormLeadingOnesSafeDivision(this);
 
 
-
         /// <summary>
         /// Finds the reduced echelon form of a matrix.
         /// </summary>
@@ -228,6 +261,16 @@ public GenTensor<T, TWrapper> ReducedRowEchelonFormSimple()
         /// </summary>
         public GenTensor<T, TWrapper> ReducedRowEchelonFormSafeDivision()
             => EchelonForm<T, TWrapper>.ReducedRowEchelonFormSafeDivision(this);
+        
+        /// <summary>
+        /// Finds the reduced echelon form of a matrix.
+        /// Uses safe division, i. e. perform division only when computing the final result.
+        ///
+        /// Additionally returns row permutation
+        /// </summary>
+        public (GenTensor<T, TWrapper>, int[]) ReducedRowEchelonFormPermuteSafeDivision()
+            => EchelonForm<T, TWrapper>.ReducedRowEchelonFormPermuteSafeDivision(this);
+
         #endregion
 
         #region Determinant
@@ -307,6 +350,7 @@ public void InvertMatrix()
         /// </summary>
         public void TensorMatrixInvert()
             => Inversion<T, TWrapper>.TensorMatrixInvert(this);
+
         #endregion
 
         #region Elementary matrix operations
@@ -398,70 +442,80 @@ public static GenTensor<T, TWrapper> TensorMatrixMultiply(GenTensor<T, TWrapper>
         /// T1 + T2
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseAdd(a, b, threading);
 
         /// <summary>
         /// T1 - T2
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseSubtract(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseSubtract(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseSubtract(a, b, threading);
 
         /// <summary>
         /// T1 * T2
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseMultiply(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseMultiply(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseMultiply(a, b, threading);
 
         /// <summary>
         /// T1 / T2
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseDivide(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseDivide(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseDivide(a, b, threading);
 
         /// <summary>
         /// T1 + const
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, T b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseAdd(GenTensor<T, TWrapper> a, T b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseAdd(a, b, threading);
 
         /// <summary>
         /// T1 - const
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseSubtract(GenTensor<T, TWrapper> a, T b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseSubtract(GenTensor<T, TWrapper> a, T b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseSubtract(a, b, threading);
 
         /// <summary>
         /// const - T1
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseSubtract(T a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseSubtract(T a, GenTensor<T, TWrapper> b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseSubtract(a, b, threading);
 
         /// <summary>
         /// T1 * const
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseMultiply(GenTensor<T, TWrapper> a, T b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseMultiply(GenTensor<T, TWrapper> a, T b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseMultiply(a, b, threading);
 
         /// <summary>
         /// T1 / const
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseDivide(GenTensor<T, TWrapper> a, T b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseDivide(GenTensor<T, TWrapper> a, T b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseDivide(a, b, threading);
 
         /// <summary>
         /// const / T1
         /// </summary>
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
-        public static GenTensor<T, TWrapper> PiecewiseDivide(T a, GenTensor<T, TWrapper> b, Threading threading = Threading.Single)
+        public static GenTensor<T, TWrapper> PiecewiseDivide(T a, GenTensor<T, TWrapper> b,
+            Threading threading = Threading.Single)
             => PiecewiseArithmetics<T, TWrapper>.PiecewiseDivide(a, b, threading);
 
         #endregion
@@ -493,6 +547,40 @@ public GenTensor<T, TWrapper> TensorMatrixPower(int power, Threading threading =
 
         #endregion
 
+        #region Decompositions
+
+        /// <summary>
+        /// https://www.geeksforgeeks.org/l-u-decomposition-system-linear-equations/
+        /// </summary>
+        /// <returns>
+        /// LU decomposition
+        /// </returns>
+        public (GenTensor<T, TWrapper> l, GenTensor<T, TWrapper> u) LuDecomposition()
+            => LuDecomposition<T, TWrapper>.Decompose(this);
+        
+        /// <summary>
+        /// Find PLU decomposition: matrices P, L, U such that for original matrix A: PA = LU.
+        /// 
+        /// P stands for permutation matrix, permutations are made during the Gauss elimination process
+        /// L stands for LIBERTY lower triangle matrix
+        /// U stands for upper triangle matrix
+        ///
+        /// Algorithm, given matrix A:
+        /// 1. Form an adjacent matrix (A|E)
+        /// 2. Find row echelon form of that matrix (U|L_0) and permutation of the rows
+        /// 3. Form permutation matrix P such that P_ij = \delta_{}
+        /// 4. Compute L = P * L_0^{-1}
+        ///
+        /// Results are: P, L, U
+        /// </summary>
+        /// <returns>
+        /// LUP decomposition of given matrix
+        /// </returns>
+        public (GenTensor<T, TWrapper> p, GenTensor<T, TWrapper> l, GenTensor<T, TWrapper> u) PluDecomposition()
+            => PluDecomposition<T, TWrapper>.Decompose(this);
+
+        #endregion
+
         #region ToString & GetHashCode
 
         /// <inheritdoc/>
@@ -517,7 +605,8 @@ public static GenTensor<T, TWrapper> VectorCrossProduct(GenTensor<T, TWrapper> a
         /// <summary>
         /// Calls VectorCrossProduct for every vector in the tensor
         /// </summary>
-        public static GenTensor<T, TWrapper> TensorVectorCrossProduct(GenTensor<T, TWrapper> a, GenTensor<T, TWrapper> b)
+        public static GenTensor<T, TWrapper> TensorVectorCrossProduct(GenTensor<T, TWrapper> a,
+            GenTensor<T, TWrapper> b)
             => VectorProduct<T, TWrapper>.TensorVectorCrossProduct(a, b);
 
         /// <summary>
@@ -562,6 +651,7 @@ public GenTensor<T, TWrapper> Forward()
         #endregion
 
         #region Serialization
+
         /*
          * Serialization protocol:
          *
@@ -595,6 +685,7 @@ public byte[] Serialize()
         /// </returns>
         public static GenTensor<T, TWrapper> Deserialize(byte[] data)
             => Serializer<T, TWrapper>.Deserialize(data);
+
         #endregion
     }
-}
+}
\ No newline at end of file
diff --git a/GenericTensor/Functions/Determinant.cs b/GenericTensor/Functions/Determinant.cs
index c3dce7a..d67448b 100644
--- a/GenericTensor/Functions/Determinant.cs
+++ b/GenericTensor/Functions/Determinant.cs
@@ -87,7 +87,7 @@ internal static T DeterminantGaussianSafeDivision(GenTensor<T, TWrapper> t, int
                 return t.GetValueNoCheck(0, 0);
 
             var n = diagLength;
-            var elemMatrix = EchelonForm<T, TWrapper>.InnerGaussianEliminationSafeDivision(t, n, n, out var swapCount);
+            var elemMatrix = EchelonForm<T, TWrapper>.InnerGaussianEliminationSafeDivision(t, n, n, null, out var swapCount);
 
             var det = default(EchelonForm<T, TWrapper>.WrapperSafeDivisionWrapper<T, TWrapper>).CreateOne();
             for (int i = 0; i < n; i++)
diff --git a/GenericTensor/Functions/EchelonForm.cs b/GenericTensor/Functions/EchelonForm.cs
index b782158..7432d0b 100644
--- a/GenericTensor/Functions/EchelonForm.cs
+++ b/GenericTensor/Functions/EchelonForm.cs
@@ -1,4 +1,5 @@
 #region copyright
+
 /*
  * MIT License
  * 
@@ -22,6 +23,7 @@
  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  * SOFTWARE.
  */
+
 #endregion
 
 
@@ -34,11 +36,18 @@ namespace GenericTensor.Functions
 {
     internal static class EchelonFormExtensions
     {
-        internal static GenTensor<T, TWrapper> SafeDivisionToSimple<T, TWrapper>(this GenTensor<EchelonForm<T, TWrapper>.SafeDivisionWrapper<T, TWrapper>, EchelonForm<T, TWrapper>.WrapperSafeDivisionWrapper<T, TWrapper>> t) where TWrapper : struct, IOperations<T>
+        internal static GenTensor<T, TWrapper> SafeDivisionToSimple<T, TWrapper>(
+            this GenTensor<EchelonForm<T, TWrapper>.SafeDivisionWrapper<T, TWrapper>,
+                EchelonForm<T, TWrapper>.WrapperSafeDivisionWrapper<T, TWrapper>> t)
+            where TWrapper : struct, IOperations<T>
             => GenTensor<T, TWrapper>.CreateMatrix(t.Shape[0], t.Shape[1], (x, y) => t.GetValueNoCheck(x, y).Count());
 
-        internal static GenTensor<EchelonForm<T, TWrapper>.SafeDivisionWrapper<T, TWrapper>, EchelonForm<T, TWrapper>.WrapperSafeDivisionWrapper<T, TWrapper>> SimpleToSafeDivision<T, TWrapper>(this GenTensor<T, TWrapper> t) where TWrapper : struct, IOperations<T>
-            => GenTensor<EchelonForm<T, TWrapper>.SafeDivisionWrapper<T, TWrapper>, EchelonForm<T, TWrapper>.WrapperSafeDivisionWrapper<T, TWrapper>>
+        internal static
+            GenTensor<EchelonForm<T, TWrapper>.SafeDivisionWrapper<T, TWrapper>,
+                EchelonForm<T, TWrapper>.WrapperSafeDivisionWrapper<T, TWrapper>>
+            SimpleToSafeDivision<T, TWrapper>(this GenTensor<T, TWrapper> t) where TWrapper : struct, IOperations<T>
+            => GenTensor<EchelonForm<T, TWrapper>.SafeDivisionWrapper<T, TWrapper>,
+                    EchelonForm<T, TWrapper>.WrapperSafeDivisionWrapper<T, TWrapper>>
                 .CreateMatrix(t.Shape[0], t.Shape[1],
                     (x, y) => new EchelonForm<T, TWrapper>.SafeDivisionWrapper<T, TWrapper>(t.GetValueNoCheck(x, y))
                 );
@@ -47,6 +56,7 @@ internal static GenTensor<EchelonForm<T, TWrapper>.SafeDivisionWrapper<T, TWrapp
     internal static class EchelonForm<T, TWrapper> where TWrapper : struct, IOperations<T>
     {
         #region Gaussian elimination safe division
+
         internal struct SafeDivisionWrapper<W, TW> where TW : struct, IOperations<W>
         {
             internal W num;
@@ -63,28 +73,31 @@ public SafeDivisionWrapper(W num, W den)
                 this.num = num;
                 this.den = den;
             }
-            
+
             public W Count() => default(TW).Divide(num, den);
         }
 
-        internal struct WrapperSafeDivisionWrapper<W, TW> : IOperations<SafeDivisionWrapper<W, TW>> where TW : struct, IOperations<W>
+        internal struct WrapperSafeDivisionWrapper<W, TW> : IOperations<SafeDivisionWrapper<W, TW>>
+            where TW : struct, IOperations<W>
         {
-
             public SafeDivisionWrapper<W, TW> Add(SafeDivisionWrapper<W, TW> a, SafeDivisionWrapper<W, TW> b)
             {
-                return new SafeDivisionWrapper<W, TW>(default(TW).Add(default(TW).Multiply(a.num, b.den), default(TW).Multiply(b.num, a.den)),
+                return new SafeDivisionWrapper<W, TW>(
+                    default(TW).Add(default(TW).Multiply(a.num, b.den), default(TW).Multiply(b.num, a.den)),
                     default(TW).Multiply(a.den, b.den));
             }
 
             public SafeDivisionWrapper<W, TW> Subtract(SafeDivisionWrapper<W, TW> a, SafeDivisionWrapper<W, TW> b)
             {
-                return new SafeDivisionWrapper<W, TW>(default(TW).Subtract(default(TW).Multiply(a.num, b.den), default(TW).Multiply(b.num, a.den)),
+                return new SafeDivisionWrapper<W, TW>(
+                    default(TW).Subtract(default(TW).Multiply(a.num, b.den), default(TW).Multiply(b.num, a.den)),
                     default(TW).Multiply(a.den, b.den));
             }
 
             public SafeDivisionWrapper<W, TW> Multiply(SafeDivisionWrapper<W, TW> a, SafeDivisionWrapper<W, TW> b)
             {
-                return new SafeDivisionWrapper<W, TW>(default(TW).Multiply(a.num, b.num), default(TW).Multiply(a.den, b.den));
+                return new SafeDivisionWrapper<W, TW>(default(TW).Multiply(a.num, b.num),
+                    default(TW).Multiply(a.den, b.den));
             }
 
             public SafeDivisionWrapper<W, TW> Negate(SafeDivisionWrapper<W, TW> a)
@@ -94,7 +107,8 @@ public SafeDivisionWrapper<W, TW> Negate(SafeDivisionWrapper<W, TW> a)
 
             public SafeDivisionWrapper<W, TW> Divide(SafeDivisionWrapper<W, TW> a, SafeDivisionWrapper<W, TW> b)
             {
-                return new SafeDivisionWrapper<W, TW>(default(TW).Multiply(a.num, b.den), default(TW).Multiply(a.den, b.num));
+                return new SafeDivisionWrapper<W, TW>(default(TW).Multiply(a.num, b.den),
+                    default(TW).Multiply(a.den, b.num));
             }
 
             public SafeDivisionWrapper<W, TW> CreateOne()
@@ -143,31 +157,33 @@ public SafeDivisionWrapper<W, TW> Deserialize(byte[] data)
             }
         }
 
-        internal static GenTensor<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>> InnerGaussianEliminationSafeDivision(GenTensor<T, TWrapper> t, int m, int n, out int swapCount)
+        internal static GenTensor<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>>
+            InnerGaussianEliminationSafeDivision(GenTensor<T, TWrapper> t, int m, int n, int[]? permutations, out int swapCount)
         {
             var elemMatrix = t.SimpleToSafeDivision();
             swapCount = 0;
-            EchelonForm<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>>.InnerGaussianEliminationSimple(elemMatrix, 0, ref swapCount);
+            EchelonForm<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>>
+                .InnerGaussianEliminationSimple(elemMatrix, 0, permutations, ref swapCount);
             return elemMatrix;
         }
 
         public static GenTensor<T, TWrapper> RowEchelonFormSafeDivision(GenTensor<T, TWrapper> t)
         {
-            #if ALLOW_EXCEPTIONS
+#if ALLOW_EXCEPTIONS
             if (!t.IsMatrix)
                 throw new InvalidShapeException("this should be matrix");
-            #endif
-            var wrp = InnerGaussianEliminationSafeDivision(t, t.Shape[0], t.Shape[1], out _);
+#endif
+            var wrp = InnerGaussianEliminationSafeDivision(t, t.Shape[0], t.Shape[1], null, out _);
             return wrp.SafeDivisionToSimple();
         }
 
         internal static void InnerGaussianEliminationSimpleDiscardSwapCount(GenTensor<T, TWrapper> t, int off)
         {
             var intoNowhere = 0;
-            InnerGaussianEliminationSimple(t, off, ref intoNowhere);
+            InnerGaussianEliminationSimple(t, off, null, ref intoNowhere);
         }
 
-        internal static void InnerGaussianEliminationSimple(GenTensor<T, TWrapper> t, int off, ref int swapCount)
+        internal static void InnerGaussianEliminationSimple(GenTensor<T, TWrapper> t, int off, int[]? permutations, ref int swapCount)
         {
             // Here we are sticking to the algorithm,
             // provided here: https://www.math.purdue.edu/~shao92/documents/Algorithm%20REF.pdf
@@ -185,6 +201,9 @@ internal static void InnerGaussianEliminationSimple(GenTensor<T, TWrapper> t, in
             {
                 t.RowSwap(off, pivotId);
                 swapCount++;
+
+                if (permutations is not null)
+                    (permutations[off], permutations[pivotId]) = (permutations[pivotId], permutations[off]);
             }
 
 
@@ -200,7 +219,7 @@ internal static void InnerGaussianEliminationSimple(GenTensor<T, TWrapper> t, in
 
 
             // VI. Let us apply the algorithm for the inner matrix
-            InnerGaussianEliminationSimple(t, off + 1, ref swapCount);
+            InnerGaussianEliminationSimple(t, off + 1, permutations, ref swapCount);
 
 
             static int? NonZeroColumn(GenTensor<T, TWrapper> t, int c, int off)
@@ -223,18 +242,47 @@ internal static void InnerGaussianEliminationSimple(GenTensor<T, TWrapper> t, in
 
         public static GenTensor<T, TWrapper> RowEchelonFormSimple(GenTensor<T, TWrapper> t)
         {
-            #if ALLOW_EXCEPTIONS
+#if ALLOW_EXCEPTIONS
             if (!t.IsMatrix)
                 throw new InvalidShapeException("this should be matrix");
-            #endif
+#endif
             var res = t.Copy(copyElements: false);
             InnerGaussianEliminationSimpleDiscardSwapCount(res, 0);
             return res;
         }
 
+        public static (GenTensor<T, TWrapper>, int[]) RowEchelonFormPermuteSimple(GenTensor<T, TWrapper> t)
+        {
+#if ALLOW_EXCEPTIONS
+            if (!t.IsMatrix)
+                throw new InvalidShapeException("this should be matrix");
+#endif
+            var res = t.Copy(copyElements: false);
+            var permute = new int[t.Shape[0]];
+            for (var i = 0; i < permute.Length; i++) permute[i] = i + 1;
+
+            var _ = 0;
+            InnerGaussianEliminationSimple(res, 0, permute, ref _);
+            return (res, permute);
+        }
+
+        public static (GenTensor<T, TWrapper>, int[]) RowEchelonFormPermuteSafeDivision(GenTensor<T, TWrapper> t)
+        {
+#if ALLOW_EXCEPTIONS
+            if (!t.IsMatrix)
+                throw new InvalidShapeException("this should be matrix");
+#endif
+            var permutations = new int[t.Shape[0]];
+            for (var i = 0; i < permutations.Length; i++) permutations[i] = i + 1;
+
+            var wrp = InnerGaussianEliminationSafeDivision(t, t.Shape[0], t.Shape[1], permutations, out _);
+            return (wrp.SafeDivisionToSimple(), permutations);
+        }
+
         #endregion
 
         #region Row echelon form leading ones
+
         private static GenTensor<T, TWrapper> InnerRowEchelonFormLeadingOnes(GenTensor<T, TWrapper> t)
         {
             var rowForm = t.Copy(copyElements: false);
@@ -247,37 +295,39 @@ private static GenTensor<T, TWrapper> InnerRowEchelonFormLeadingOnes(GenTensor<T
 
         public static GenTensor<T, TWrapper> RowEchelonFormLeadingOnesSimple(GenTensor<T, TWrapper> t)
         {
-            #if ALLOW_EXCEPTIONS
+#if ALLOW_EXCEPTIONS
             if (!t.IsMatrix)
                 throw new InvalidShapeException("this should be matrix");
-            #endif
+#endif
             return InnerRowEchelonFormLeadingOnes(t);
         }
 
         public static GenTensor<T, TWrapper> RowEchelonFormLeadingOnesSafeDivision(GenTensor<T, TWrapper> t)
         {
-            #if ALLOW_EXCEPTIONS
+#if ALLOW_EXCEPTIONS
             if (!t.IsMatrix)
                 throw new InvalidShapeException("this should be matrix");
-            #endif
-            return EchelonForm<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>>.InnerRowEchelonFormLeadingOnes(t.SimpleToSafeDivision()).SafeDivisionToSimple();
+#endif
+            return EchelonForm<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>>
+                .InnerRowEchelonFormLeadingOnes(t.SimpleToSafeDivision()).SafeDivisionToSimple();
         }
 
         #endregion
 
         #region Reduced row echelon form
 
-        private static GenTensor<T, TWrapper> InnerReducedRowEchelonForm(GenTensor<T, TWrapper> t, out int swapCount)
+        private static GenTensor<T, TWrapper> InnerReducedRowEchelonForm(GenTensor<T, TWrapper> t, int[]? permutations, out int swapCount)
         {
             var upper = t.Copy(copyElements: false);
             swapCount = 0;
-            InnerGaussianEliminationSimple(upper, 0, ref swapCount);
+            InnerGaussianEliminationSimple(upper, 0, permutations, ref swapCount);
             for (int r = t.Shape[0] - 1; r >= 0; r--)
             {
                 if (upper.RowGetLeadingElement(r) is not { } leading)
                     continue;
                 for (int i = 0; i < r; i++)
-                    upper.RowSubtract(i, r, default(TWrapper).Divide(upper.GetValueNoCheck(i, leading.index), leading.value));
+                    upper.RowSubtract(i, r,
+                        default(TWrapper).Divide(upper.GetValueNoCheck(i, leading.index), leading.value));
 
                 upper.RowMultiply(r, default(TWrapper).Divide(default(TWrapper).CreateOne(), leading.value));
             }
@@ -287,22 +337,38 @@ private static GenTensor<T, TWrapper> InnerReducedRowEchelonForm(GenTensor<T, TW
 
         public static GenTensor<T, TWrapper> ReducedRowEchelonFormSimple(GenTensor<T, TWrapper> t)
         {
-            #if ALLOW_EXCEPTIONS
+#if ALLOW_EXCEPTIONS
             if (!t.IsMatrix)
                 throw new InvalidShapeException("this should be matrix");
-            #endif
-            return InnerReducedRowEchelonForm(t, out _);
+#endif
+            return InnerReducedRowEchelonForm(t, null, out _);
         }
 
         public static GenTensor<T, TWrapper> ReducedRowEchelonFormSafeDivision(GenTensor<T, TWrapper> t)
         {
+#if ALLOW_EXCEPTIONS
+            if (!t.IsMatrix)
+                throw new InvalidShapeException("this should be matrix");
+#endif
+            return EchelonForm<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>>
+                .InnerReducedRowEchelonForm(t.SimpleToSafeDivision(), null, out var _).SafeDivisionToSimple();
+        }
+
+        public static (GenTensor<T, TWrapper> result, int[] permutations) ReducedRowEchelonFormPermuteSafeDivision(GenTensor<T, TWrapper> t)
+        {
             #if ALLOW_EXCEPTIONS
             if (!t.IsMatrix)
                 throw new InvalidShapeException("this should be matrix");
             #endif
-            return EchelonForm<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>>.InnerReducedRowEchelonForm(t.SimpleToSafeDivision(), out var _).SafeDivisionToSimple();
+            
+            var permutation = new int[t.Shape[0]];
+            for (var i = 0; i < permutation.Length; i++) permutation[i] = i + 1;
+
+            return (EchelonForm<SafeDivisionWrapper<T, TWrapper>, WrapperSafeDivisionWrapper<T, TWrapper>>
+                    .InnerReducedRowEchelonForm(t.SimpleToSafeDivision(), permutation, out var _).SafeDivisionToSimple(),
+                permutation);
         }
 
         #endregion
     }
-}
+}
\ No newline at end of file
diff --git a/GenericTensor/Functions/LuDecomposition.cs b/GenericTensor/Functions/LuDecomposition.cs
new file mode 100644
index 0000000..5899cc2
--- /dev/null
+++ b/GenericTensor/Functions/LuDecomposition.cs
@@ -0,0 +1,90 @@
+#region copyright
+
+/*
+ * MIT License
+ * 
+ * Copyright (c) 2020-2021 WhiteBlackGoose
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ */
+
+#endregion
+
+
+using GenericTensor.Core;
+
+namespace GenericTensor.Functions
+{
+    internal static class LuDecomposition<T, TWrapper> where TWrapper : struct, IOperations<T>
+    {
+        public static (GenTensor<T, TWrapper>, GenTensor<T, TWrapper>) Decompose(GenTensor<T, TWrapper> t)
+        {
+#if ALLOW_EXCEPTIONS
+            if (!t.IsSquareMatrix)
+                throw new InvalidShapeException("this should be a square matrix");
+#endif
+
+            var n = t.Shape[0];
+            var lower = GenTensor<T, TWrapper>.CreateMatrix(n, n);
+            var upper = GenTensor<T, TWrapper>.CreateMatrix(n, n);
+
+            var tw = default(TWrapper);
+            var zero = tw.CreateZero();
+            var one = tw.CreateOne();
+
+            for (var i = 0; i < n; i++)
+            {
+                // Upper triangular
+                for (var k = i; k < n; k++)
+                {
+                    var sum = zero;
+                    for (var j = 0; j < i; j++)
+                        sum = tw.Add(sum, tw.Multiply(lower[i, j], upper[j, k]));
+
+                    upper[i, k] = tw.Subtract(t[i, k], sum);
+                }
+
+                // Lower triangular
+                for (var k = i; k < n; k++)
+                {
+                    if (i == k)
+                        lower[i, i] = one; // 1 at diagonals
+                    else
+                    {
+                        var sum = zero;
+                        for (var j = 0; j < i; j++)
+                            sum = tw.Add(sum, tw.Multiply(lower[k, j], upper[j, i]));
+
+                        #if ALLOW_EXCEPTIONS
+                        if (Equals(upper[i, i], zero))
+                        {
+                            throw new ImpossibleDecomposition("There is no LU decomposition for given matrix");
+                        }
+                        #endif
+                        
+                        lower[k, i]
+                            = tw.Divide(tw.Subtract(t[k, i], sum), upper[i, i]);
+                    }
+                }
+            }
+
+            return (lower, upper);
+        }
+    }
+}
\ No newline at end of file
diff --git a/GenericTensor/Functions/PluDecomposition.cs b/GenericTensor/Functions/PluDecomposition.cs
new file mode 100644
index 0000000..d7356a7
--- /dev/null
+++ b/GenericTensor/Functions/PluDecomposition.cs
@@ -0,0 +1,68 @@
+#region copyright
+
+/*
+ * MIT License
+ * 
+ * Copyright (c) 2020-2021 WhiteBlackGoose
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ */
+
+#endregion
+
+
+using GenericTensor.Core;
+
+namespace GenericTensor.Functions
+{
+    internal static class PluDecomposition<T, TWrapper> where TWrapper : struct, IOperations<T>
+    {
+        public static (GenTensor<T, TWrapper>, GenTensor<T, TWrapper>, GenTensor<T, TWrapper>) Decompose(GenTensor<T, TWrapper> original)
+        {
+            var t = original.Copy(copyElements: true);
+
+            #if ALLOW_EXCEPTIONS
+            if (!t.IsSquareMatrix)
+                throw new InvalidShapeException("this should be a square matrix");
+            #endif
+
+            var n = t.Shape[0];
+            var m = t.Shape[1];
+            var tw = default(TWrapper);
+
+            var identity = GenTensor<T, TWrapper>.CreateIdentityMatrix(m);
+
+            t.TransposeMatrix();
+            var adj = GenTensor<T, TWrapper>.Concat(t, identity);
+            adj.TransposeMatrix();
+
+            var (echelon, permute) = adj.RowEchelonFormPermuteSafeDivision();
+            var upper = GenTensor<T, TWrapper>.CreateMatrix(n, m, (i, j) => echelon[i, j]);
+            var lowerZero = GenTensor<T, TWrapper>.CreateMatrix(m, m, (i, j) => echelon[i, m + j]);
+
+            lowerZero.InvertMatrix();
+
+            var permuteMatrix =
+                GenTensor<T, TWrapper>.CreateMatrix(m, m, (i, j) => j == permute[i] - 1 ? tw.CreateOne() : tw.CreateZero());
+
+            var lower = GenTensor<T, TWrapper>.MatrixMultiply(permuteMatrix, lowerZero);
+            return (permuteMatrix, lower, upper);
+        }
+    }
+}
\ No newline at end of file
diff --git a/GenericTensor/GenericTensor.csproj b/GenericTensor/GenericTensor.csproj
index bb9cd57..c168d09 100644
--- a/GenericTensor/GenericTensor.csproj
+++ b/GenericTensor/GenericTensor.csproj
@@ -18,6 +18,7 @@
     <PackageTags>tensor, generic, matrix, vector, performance</PackageTags>
     <SignAssembly>true</SignAssembly>
     <AssemblyOriginatorKeyFile>keypair.snk</AssemblyOriginatorKeyFile>
+    <Nullable>enable</Nullable>
   </PropertyGroup>
 
   <PropertyGroup Condition="'$(Configuration)|$(Platform)'=='Release|AnyCPU'">
diff --git a/GenericTensor/GenericTensor.csproj.DotSettings b/GenericTensor/GenericTensor.csproj.DotSettings
deleted file mode 100644
index b9fd6ee..0000000
--- a/GenericTensor/GenericTensor.csproj.DotSettings
+++ /dev/null
@@ -1,2 +0,0 @@
-<wpf:ResourceDictionary xml:space="preserve" xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml" xmlns:s="clr-namespace:System;assembly=mscorlib" xmlns:ss="urn:shemas-jetbrains-com:settings-storage-xaml" xmlns:wpf="http://schemas.microsoft.com/winfx/2006/xaml/presentation">
-	<s:String x:Key="/Default/CodeInspection/CSharpLanguageProject/LanguageLevel/@EntryValue">CSharp80</s:String></wpf:ResourceDictionary>
\ No newline at end of file
diff --git a/README.md b/README.md
index 6efd483..eda315f 100644
--- a/README.md
+++ b/README.md
@@ -352,6 +352,32 @@ Finds the power of a matrix.
 Works for O(log(N) * N^3)
 </p></details>
 
+#### Decompositions
+
+<details><summary><strong>PLU decomposision</strong></summary><p>
+
+```cs     
+public (GenTensor<T, TWrapper> p, GenTensor<T, TWrapper> l, GenTensor<T, TWrapper> u) PluDecomposition()
+```
+
+Find PLU decomposition: matrices P, L, U such that for original matrix A: PA = LU.
+
+Works for O(n^3)
+
+</details>
+
+<details><summary><strong>LU decomposision</strong></summary><p>
+
+```cs
+public (GenTensor<T, TWrapper>, GenTensor<T, TWrapper>) LuDecomposition()   
+```
+
+Find LU decomposition: matrices L, U such that for original matrix A: A = LU.
+
+Works for O(n^3)
+
+</details>
+
 #### Piecewise arithmetics
 
 <details><summary><strong>Tensor and Tensor</strong></summary><p>
diff --git a/UnitTests/LuDecompositions.cs b/UnitTests/LuDecompositions.cs
new file mode 100644
index 0000000..9118bfd
--- /dev/null
+++ b/UnitTests/LuDecompositions.cs
@@ -0,0 +1,136 @@
+#region copyright
+
+/*
+ * MIT License
+ * 
+ * Copyright (c) 2020-2021 WhiteBlackGoose
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ * 
+ * The above copyright notice and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ */
+
+#endregion
+
+using System;
+using GenericTensor.Core;
+using GenericTensor.Functions;
+using Microsoft.VisualStudio.TestTools.UnitTesting;
+
+namespace UnitTests
+{
+    [TestClass]
+    public class LuDecompositions
+    {
+        private void AssertMatrixIsLowerTriangular<T, TW>(GenTensor<T, TW> t) where TW : struct, IOperations<T>
+        {
+            Assert.IsTrue(t.IsSquareMatrix);
+            t.IterateOver2((i, j) =>
+            {
+                if (i < j) Assert.AreEqual(t[i, j], default(TW).CreateZero());
+            });
+        }
+
+        private void AssertMatrixIsUpperTriangular<T, TW>(GenTensor<T, TW> t) where TW : struct, IOperations<T>
+        {
+            Assert.IsTrue(t.IsSquareMatrix);
+            t.IterateOver2((i, j) =>
+            {
+                if (i > j) Assert.AreEqual(t[i, j], default(TW).CreateZero());
+            });
+        }
+
+        private void TestProvidedLu<T, TW>(T[,] data) where TW : struct, IOperations<T>
+        {
+            var m = GenTensor<T, TW>.CreateMatrix(data);
+            var (lower, upper) = m.LuDecomposition();
+            
+            AssertMatrixIsLowerTriangular(lower);
+            AssertMatrixIsUpperTriangular(upper);
+
+            Assert.AreEqual(GenTensor<T, TW>.MatrixMultiply(lower, upper), m);
+        }
+        
+        private void TestProvidedPlu<T, TW>(T[,] data) where TW : struct, IOperations<T>
+        {
+            var m = GenTensor<T, TW>.CreateMatrix(data);
+            var (p, l, u) = m.PluDecomposition();
+            
+            AssertMatrixIsLowerTriangular(l);
+            AssertMatrixIsUpperTriangular(u);
+
+            Assert.AreEqual(GenTensor<T, TW>.MatrixMultiply(p, m), GenTensor<T, TW>.MatrixMultiply(l, u));
+        }
+
+        private static void TestNoLu<T, TW>(T[,] data) where TW : struct, IOperations<T> =>
+            Assert.ThrowsException<ImpossibleDecomposition>(() =>
+                GenTensor<T, TW>.CreateMatrix(data).LuDecomposition());
+
+        [TestMethod]
+        public void LuArbitrary_1() =>
+            TestProvidedLu<int, IntWrapper>(new[,]
+            {
+                { 1, 2, 4 },
+                { 3, 8, 14 },
+                { 2, 6, 13 },
+            });
+
+        [TestMethod]
+        public void LuArbitrary_2() =>
+            TestProvidedLu<int, IntWrapper>(new[,]
+            {
+                { -6, -4 },
+                { -6, -4 }
+            });
+
+        [TestMethod]
+        public void LuArbitrary_3() =>
+            TestNoLu<int, IntWrapper>(new[,]
+            {
+                { 0, 0, 0 },
+                { 0, 0, 0 },
+                { 0, 0, 0 }
+            });
+
+        [TestMethod]
+        public void PluArbitrary_1() =>
+            TestProvidedPlu<int, IntWrapper>(new[,]
+            {
+                { 1, 2, 4 },
+                { 3, 8, 14 },
+                { 2, 6, 13 }
+            });
+        
+        [TestMethod]
+        public void PluArbitrary_2() =>
+            TestProvidedPlu<double, DoubleWrapper>(new[,]
+            {
+                { 0D, 1, 0 },
+                { -8, 8, 1 },
+                { 2, -2, 0 }
+            });
+        
+        [TestMethod]
+        public void PluArbitrary_3() =>
+            TestProvidedPlu<int, IntWrapper>(new[,]
+            {
+                {1, 0, 4},
+                {0, 0, 1},
+                {3, 2, 13}
+            });
+    }
+}
\ No newline at end of file
diff --git a/UnitTests/UnitTests.csproj b/UnitTests/UnitTests.csproj
index 9dc2d9a..3954fdb 100644
--- a/UnitTests/UnitTests.csproj
+++ b/UnitTests/UnitTests.csproj
@@ -1,9 +1,11 @@
 <Project Sdk="Microsoft.NET.Sdk">
 
   <PropertyGroup>
-    <TargetFramework>net5.0</TargetFramework>
+    <TargetFramework>netcoreapp3.1</TargetFramework>
 
     <IsPackable>false</IsPackable>
+
+    <LangVersion>9</LangVersion>
   </PropertyGroup>
   
   <ItemGroup>