diff --git a/base/serial/impl/psb_c_coo_impl.F90 b/base/serial/impl/psb_c_coo_impl.F90
index ccf7e42d..2a6aaf13 100644
--- a/base/serial/impl/psb_c_coo_impl.F90
+++ b/base/serial/impl/psb_c_coo_impl.F90
@@ -595,12 +595,16 @@ subroutine  psb_c_coo_clean_zeros(a, info)
   integer(psb_ipk_), intent(out) :: info
   !
   integer(psb_ipk_) :: i,j,k, nzin
-
+  logical           :: cpy
+  
   info = 0
   nzin = a%get_nzeros()
   j = 0
   do i=1, nzin
-    if (a%val(i) /= czero) then
+    cpy = (a%val(i) /= czero)
+    ! Always keep the diagonal, even if numerically zero
+    if (.not.cpy) cpy = (a%ia(i) == a%ja(i))
+    if (cpy) then
       j = j + 1
       a%val(j) = a%val(i)
       a%ia(j)  = a%ia(i)
diff --git a/base/serial/impl/psb_c_csc_impl.F90 b/base/serial/impl/psb_c_csc_impl.F90
index fe7227dd..f8ab1025 100644
--- a/base/serial/impl/psb_c_csc_impl.F90
+++ b/base/serial/impl/psb_c_csc_impl.F90
@@ -2412,6 +2412,7 @@ subroutine  psb_c_csc_clean_zeros(a, info)
   !
   integer(psb_ipk_) :: i, j, k, nc
   integer(psb_ipk_), allocatable :: ilcp(:)
+  logical           :: cpy
 
   info = 0
   call a%sync()
@@ -2421,7 +2422,10 @@ subroutine  psb_c_csc_clean_zeros(a, info)
   j        = a%icp(1)
   do i=1, nc
     do k = ilcp(i), ilcp(i+1) -1
-      if (a%val(k) /= czero) then
+      cpy = (a%val(k) /= czero)
+      ! Always keep the diagonal, even if numerically zero
+      if (.not.cpy) cpy = (i == a%ia(k))
+      if (cpy) then
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
         j = j + 1
diff --git a/base/serial/impl/psb_c_csr_impl.F90 b/base/serial/impl/psb_c_csr_impl.F90
index 276d3d1c..d3494f46 100644
--- a/base/serial/impl/psb_c_csr_impl.F90
+++ b/base/serial/impl/psb_c_csr_impl.F90
@@ -3633,6 +3633,7 @@ subroutine  psb_c_csr_clean_zeros(a, info)
   !
   integer(psb_ipk_) :: i, j, k, nr
   integer(psb_ipk_), allocatable :: ilrp(:)
+  logical           :: cpy
 
   info = 0
   call a%sync()
@@ -3642,7 +3643,10 @@ subroutine  psb_c_csr_clean_zeros(a, info)
   j        = a%irp(1)
   do i=1, nr
     do k = ilrp(i), ilrp(i+1) -1
-      if (a%val(k) /= czero) then
+      cpy = (a%val(k) /= czero)
+      ! Always keep the diagonal, even if numerically zero
+      if (.not.cpy) cpy = (i == a%ja(k))
+      if (cpy) then
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1
diff --git a/base/serial/impl/psb_d_coo_impl.F90 b/base/serial/impl/psb_d_coo_impl.F90
index beb438aa..15d4c4ad 100644
--- a/base/serial/impl/psb_d_coo_impl.F90
+++ b/base/serial/impl/psb_d_coo_impl.F90
@@ -595,12 +595,16 @@ subroutine  psb_d_coo_clean_zeros(a, info)
   integer(psb_ipk_), intent(out) :: info
   !
   integer(psb_ipk_) :: i,j,k, nzin
-
+  logical           :: cpy
+  
   info = 0
   nzin = a%get_nzeros()
   j = 0
   do i=1, nzin
-    if (a%val(i) /= dzero) then
+    cpy = (a%val(i) /= dzero)
+    ! Always keep the diagonal, even if numerically zero
+    if (.not.cpy) cpy = (a%ia(i) == a%ja(i))
+    if (cpy) then
       j = j + 1
       a%val(j) = a%val(i)
       a%ia(j)  = a%ia(i)
diff --git a/base/serial/impl/psb_d_csc_impl.F90 b/base/serial/impl/psb_d_csc_impl.F90
index d9fa2874..92a4c555 100644
--- a/base/serial/impl/psb_d_csc_impl.F90
+++ b/base/serial/impl/psb_d_csc_impl.F90
@@ -2412,6 +2412,7 @@ subroutine  psb_d_csc_clean_zeros(a, info)
   !
   integer(psb_ipk_) :: i, j, k, nc
   integer(psb_ipk_), allocatable :: ilcp(:)
+  logical           :: cpy
 
   info = 0
   call a%sync()
@@ -2421,7 +2422,10 @@ subroutine  psb_d_csc_clean_zeros(a, info)
   j        = a%icp(1)
   do i=1, nc
     do k = ilcp(i), ilcp(i+1) -1
-      if (a%val(k) /= dzero) then
+      cpy = (a%val(k) /= dzero)
+      ! Always keep the diagonal, even if numerically zero
+      if (.not.cpy) cpy = (i == a%ia(k))
+      if (cpy) then
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
         j = j + 1
diff --git a/base/serial/impl/psb_d_csr_impl.F90 b/base/serial/impl/psb_d_csr_impl.F90
index a0aaeeee..09dffe0e 100644
--- a/base/serial/impl/psb_d_csr_impl.F90
+++ b/base/serial/impl/psb_d_csr_impl.F90
@@ -3633,6 +3633,7 @@ subroutine  psb_d_csr_clean_zeros(a, info)
   !
   integer(psb_ipk_) :: i, j, k, nr
   integer(psb_ipk_), allocatable :: ilrp(:)
+  logical           :: cpy
 
   info = 0
   call a%sync()
@@ -3642,7 +3643,10 @@ subroutine  psb_d_csr_clean_zeros(a, info)
   j        = a%irp(1)
   do i=1, nr
     do k = ilrp(i), ilrp(i+1) -1
-      if (a%val(k) /= dzero) then
+      cpy = (a%val(k) /= dzero)
+      ! Always keep the diagonal, even if numerically zero
+      if (.not.cpy) cpy = (i == a%ja(k))
+      if (cpy) then
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1
diff --git a/base/serial/impl/psb_s_coo_impl.F90 b/base/serial/impl/psb_s_coo_impl.F90
index d833cf96..217b39c2 100644
--- a/base/serial/impl/psb_s_coo_impl.F90
+++ b/base/serial/impl/psb_s_coo_impl.F90
@@ -595,12 +595,16 @@ subroutine  psb_s_coo_clean_zeros(a, info)
   integer(psb_ipk_), intent(out) :: info
   !
   integer(psb_ipk_) :: i,j,k, nzin
-
+  logical           :: cpy
+  
   info = 0
   nzin = a%get_nzeros()
   j = 0
   do i=1, nzin
-    if (a%val(i) /= szero) then
+    cpy = (a%val(i) /= szero)
+    ! Always keep the diagonal, even if numerically zero
+    if (.not.cpy) cpy = (a%ia(i) == a%ja(i))
+    if (cpy) then
       j = j + 1
       a%val(j) = a%val(i)
       a%ia(j)  = a%ia(i)
diff --git a/base/serial/impl/psb_s_csc_impl.F90 b/base/serial/impl/psb_s_csc_impl.F90
index c135c9ef..ffa32f8f 100644
--- a/base/serial/impl/psb_s_csc_impl.F90
+++ b/base/serial/impl/psb_s_csc_impl.F90
@@ -2412,6 +2412,7 @@ subroutine  psb_s_csc_clean_zeros(a, info)
   !
   integer(psb_ipk_) :: i, j, k, nc
   integer(psb_ipk_), allocatable :: ilcp(:)
+  logical           :: cpy
 
   info = 0
   call a%sync()
@@ -2421,7 +2422,10 @@ subroutine  psb_s_csc_clean_zeros(a, info)
   j        = a%icp(1)
   do i=1, nc
     do k = ilcp(i), ilcp(i+1) -1
-      if (a%val(k) /= szero) then
+      cpy = (a%val(k) /= szero)
+      ! Always keep the diagonal, even if numerically zero
+      if (.not.cpy) cpy = (i == a%ia(k))
+      if (cpy) then
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
         j = j + 1
diff --git a/base/serial/impl/psb_s_csr_impl.F90 b/base/serial/impl/psb_s_csr_impl.F90
index 9d5dc145..5e18de7b 100644
--- a/base/serial/impl/psb_s_csr_impl.F90
+++ b/base/serial/impl/psb_s_csr_impl.F90
@@ -3633,6 +3633,7 @@ subroutine  psb_s_csr_clean_zeros(a, info)
   !
   integer(psb_ipk_) :: i, j, k, nr
   integer(psb_ipk_), allocatable :: ilrp(:)
+  logical           :: cpy
 
   info = 0
   call a%sync()
@@ -3642,7 +3643,10 @@ subroutine  psb_s_csr_clean_zeros(a, info)
   j        = a%irp(1)
   do i=1, nr
     do k = ilrp(i), ilrp(i+1) -1
-      if (a%val(k) /= szero) then
+      cpy = (a%val(k) /= szero)
+      ! Always keep the diagonal, even if numerically zero
+      if (.not.cpy) cpy = (i == a%ja(k))
+      if (cpy) then
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1
diff --git a/base/serial/impl/psb_z_coo_impl.F90 b/base/serial/impl/psb_z_coo_impl.F90
index ea6b2492..89950c4c 100644
--- a/base/serial/impl/psb_z_coo_impl.F90
+++ b/base/serial/impl/psb_z_coo_impl.F90
@@ -595,12 +595,16 @@ subroutine  psb_z_coo_clean_zeros(a, info)
   integer(psb_ipk_), intent(out) :: info
   !
   integer(psb_ipk_) :: i,j,k, nzin
-
+  logical           :: cpy
+  
   info = 0
   nzin = a%get_nzeros()
   j = 0
   do i=1, nzin
-    if (a%val(i) /= zzero) then
+    cpy = (a%val(i) /= zzero)
+    ! Always keep the diagonal, even if numerically zero
+    if (.not.cpy) cpy = (a%ia(i) == a%ja(i))
+    if (cpy) then
       j = j + 1
       a%val(j) = a%val(i)
       a%ia(j)  = a%ia(i)
diff --git a/base/serial/impl/psb_z_csc_impl.F90 b/base/serial/impl/psb_z_csc_impl.F90
index becb7003..a37e04e5 100644
--- a/base/serial/impl/psb_z_csc_impl.F90
+++ b/base/serial/impl/psb_z_csc_impl.F90
@@ -2412,6 +2412,7 @@ subroutine  psb_z_csc_clean_zeros(a, info)
   !
   integer(psb_ipk_) :: i, j, k, nc
   integer(psb_ipk_), allocatable :: ilcp(:)
+  logical           :: cpy
 
   info = 0
   call a%sync()
@@ -2421,7 +2422,10 @@ subroutine  psb_z_csc_clean_zeros(a, info)
   j        = a%icp(1)
   do i=1, nc
     do k = ilcp(i), ilcp(i+1) -1
-      if (a%val(k) /= zzero) then
+      cpy = (a%val(k) /= zzero)
+      ! Always keep the diagonal, even if numerically zero
+      if (.not.cpy) cpy = (i == a%ia(k))
+      if (cpy) then
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
         j = j + 1
diff --git a/base/serial/impl/psb_z_csr_impl.F90 b/base/serial/impl/psb_z_csr_impl.F90
index 7f11c3bd..57d31f80 100644
--- a/base/serial/impl/psb_z_csr_impl.F90
+++ b/base/serial/impl/psb_z_csr_impl.F90
@@ -3633,6 +3633,7 @@ subroutine  psb_z_csr_clean_zeros(a, info)
   !
   integer(psb_ipk_) :: i, j, k, nr
   integer(psb_ipk_), allocatable :: ilrp(:)
+  logical           :: cpy
 
   info = 0
   call a%sync()
@@ -3642,7 +3643,10 @@ subroutine  psb_z_csr_clean_zeros(a, info)
   j        = a%irp(1)
   do i=1, nr
     do k = ilrp(i), ilrp(i+1) -1
-      if (a%val(k) /= zzero) then
+      cpy = (a%val(k) /= zzero)
+      ! Always keep the diagonal, even if numerically zero
+      if (.not.cpy) cpy = (i == a%ja(k))
+      if (cpy) then
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1