diff --git a/Demo_LDLT_Solver.f90 b/Demo_LDLT_Solver.f90
new file mode 100644
index 0000000..26c8aa4
--- /dev/null
+++ b/Demo_LDLT_Solver.f90
@@ -0,0 +1,73 @@
+!==============================================================================!
+  subroutine Demo_LDLT_Solver
+!------------------------------------------------------------------------------!
+!----------------------------------[Modules]-----------------------------------!
+  use Matrix_Mod
+!------------------------------------------------------------------------------!
+  implicit none
+!---------------------------------[Interfaces]---------------------------------!
+  include "Print_Matrix.int"               
+  include "Print_Matrix_Compressed.int"               
+  include "Print_Vector.int"               
+  include "Create_Matrix_Compressed.int"               
+  include "LDLT_Factorization.int"  
+  include "LDLT_Solution.int"
+  include "Load_Linear_System.int"
+  include "Matrix_Vector_Multiply.int"
+  include "Vector_Vector_Dot_Product.int"
+  include "Expand_Matrix.int"
+!-----------------------------------[Locals]-----------------------------------!
+  integer           :: n
+  real, allocatable :: a_matrix(:,:), p_matrix(:,:)
+  real, allocatable :: b(:), x(:), y(:), r(:)
+  type(Matrix)      :: c_matrix
+  real              :: error
+!==============================================================================!
+
+  ! Create compressed system matrix
+  call Create_Matrix_Compressed(c_matrix, 10, 10, 10, 0)
+  n = c_matrix % n
+  if(n<=64) call Print_Matrix_Compressed("Compressed c_matrix:", c_matrix)
+ 
+  ! Create two full matrices from the compressed one
+  call Expand_Matrix(a_matrix, c_matrix)
+  call Expand_Matrix(p_matrix, c_matrix)
+  p_matrix = 0
+
+  ! Finish memory allocation
+  allocate (b(n))
+  allocate (x(n))
+  allocate (y(n))
+  allocate (r(n))
+
+  ! Fill the right hand side
+  b = 0.1
+
+  ! Just print original matrix
+  if(n<=64) call Print_Matrix("a_matrix:", a_matrix)
+
+  ! Perform LDLT factorization on the matrix to fin the lower one
+  call LDLT_Factorization(p_matrix, a_matrix)
+  if(n<=64) call Print_Matrix("p_matrix after Cholesky factorization", p_matrix)
+
+  ! Compute x
+  call LDLT_Solution(x, p_matrix, b)
+  call Print_Vector("Solution x:", x) 
+
+  ! Check result
+  call Matrix_Vector_Multiply(y, a_matrix, x)
+  call Print_Vector("Vector y should recover the source term:", y) 
+  r = b - y
+  call Vector_Vector_Dot_Product(error, r, r)
+  write(*,*) "Error: ", sqrt(error)  
+
+  ! Free memory
+  deallocate(a_matrix)
+  deallocate(p_matrix)
+  deallocate(b)
+  deallocate(x)
+  deallocate(y)
+  deallocate(r)
+  call deallocate_Matrix(c_matrix)
+
+  end subroutine Demo_LDLT_Solver
diff --git a/Demo_Incomplete_Cholesky_T_Solver.f90 b/Demo_LDLT_Solver_From_Tflows.f90
similarity index 96%
rename from Demo_Incomplete_Cholesky_T_Solver.f90
rename to Demo_LDLT_Solver_From_Tflows.f90
index a7b6167..a71e31b 100644
--- a/Demo_Incomplete_Cholesky_T_Solver.f90
+++ b/Demo_LDLT_Solver_From_Tflows.f90
@@ -1,5 +1,5 @@
 !==============================================================================!
-  subroutine Demo_Incomplete_Cholesky_T_Solver
+  subroutine Demo_LDLT_Solver_From_Tflows
 !------------------------------------------------------------------------------!
 !----------------------------------[Modules]-----------------------------------!
   use Matrix_Mod
@@ -71,4 +71,4 @@ subroutine Demo_Incomplete_Cholesky_T_Solver
   call deallocate_Matrix(a_matrix)
   call deallocate_Matrix(p_matrix)
 
-  end subroutine Demo_Incomplete_Cholesky_T_Solver
+  end subroutine Demo_LDLT_Solver_From_Tflows
diff --git a/LDLT_Factorization.f90 b/LDLT_Factorization.f90
index 1189f20..5ac592a 100644
--- a/LDLT_Factorization.f90
+++ b/LDLT_Factorization.f90
@@ -1,7 +1,7 @@
 !==============================================================================!
-  subroutine Cholesky_Factorization(f, a)
+  subroutine LDLT_Factorization(f, a)
 !------------------------------------------------------------------------------!
-!   Computes Cholesky decomposition on full matrices.                          !
+!   Computes LDL^T decomposition on full matrices.                             !
 !------------------------------------------------------------------------------!
   implicit none
 !---------------------------------[Arguments]----------------------------------!
@@ -30,4 +30,4 @@ subroutine Cholesky_Factorization(f, a)
     end do
   end do
 
-  end subroutine Cholesky_Factorization
+  end subroutine LDLT_Factorization
diff --git a/LDLT_Factorization.int b/LDLT_Factorization.int
new file mode 100644
index 0000000..52c1ab0
--- /dev/null
+++ b/LDLT_Factorization.int
@@ -0,0 +1,7 @@
+interface
+  subroutine LDLT_Factorization(f, a)
+    implicit none
+    real, dimension(:,:) :: f
+    real, dimension(:,:) :: a
+  end subroutine LDLT_Factorization
+end interface
diff --git a/LDLT_Solution.int b/LDLT_Solution.int
new file mode 100644
index 0000000..d922e41
--- /dev/null
+++ b/LDLT_Solution.int
@@ -0,0 +1,8 @@
+interface
+  subroutine LDLT_Solution(x, f, b)
+    implicit none
+    real, dimension(:)   :: x
+    real, dimension(:,:) :: f
+    real, dimension(:)   :: b
+  end subroutine LDLT_Solution
+end interface