diff --git a/base/serial/impl/psb_c_coo_impl.F90 b/base/serial/impl/psb_c_coo_impl.F90
index 5b681d22..a56a79a4 100644
--- a/base/serial/impl/psb_c_coo_impl.F90
+++ b/base/serial/impl/psb_c_coo_impl.F90
@@ -595,19 +595,13 @@ 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 = .true.
-    else
-      ! Always keep the diagonal, even if numerically zero
-      cpy = (a%ia(i) == a%ja(i))
-    end if
-    if (cpy) then
+    ! Always keep the diagonal, even if numerically zero
+    if ((a%val(i) /= czero).or.(a%ia(i) == a%ja(i))) then
       j = j + 1
       a%val(j) = a%val(i)
       a%ia(j)  = a%ia(i)
@@ -5933,19 +5927,13 @@ subroutine  psb_lc_coo_clean_zeros(a, info)
   integer(psb_ipk_), intent(out) :: info
   !
   integer(psb_lpk_) :: 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 = .true.
-    else
-      ! Always keep the diagonal, even if numerically zero
-      cpy = (a%ia(i) == a%ja(i))
-    end if
-    if (cpy) then 
+    ! Always keep the diagonal, even if numerically zero
+    if ((a%val(i) /= czero).or.(a%ia(i) == a%ja(i))) 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 190a4d5b..7916d954 100644
--- a/base/serial/impl/psb_c_csc_impl.F90
+++ b/base/serial/impl/psb_c_csc_impl.F90
@@ -2412,7 +2412,6 @@ 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()
@@ -2422,13 +2421,8 @@ 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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ia(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= czero).or.(i == a%ia(k))) then
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
         j = j + 1
@@ -4320,7 +4314,6 @@ subroutine  psb_lc_csc_clean_zeros(a, info)
   !
   integer(psb_lpk_) :: i, j, k, nc
   integer(psb_lpk_), allocatable :: ilcp(:)
-  logical :: cpy
   
   info = 0
   call a%sync()
@@ -4330,13 +4323,8 @@ subroutine  psb_lc_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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ia(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= czero).or.(i == a%ia(k))) 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 0db9f3fa..6a31e522 100644
--- a/base/serial/impl/psb_c_csr_impl.F90
+++ b/base/serial/impl/psb_c_csr_impl.F90
@@ -3633,7 +3633,6 @@ 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()
@@ -3643,13 +3642,8 @@ 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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ja(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= czero).or.(i == a%ja(k))) then
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1
@@ -6559,7 +6553,6 @@ subroutine  psb_lc_csr_clean_zeros(a, info)
   !
   integer(psb_lpk_) :: i, j, k, nr
   integer(psb_lpk_), allocatable :: ilrp(:)
-  logical :: cpy
   
   info = 0
   call a%sync()
@@ -6569,13 +6562,8 @@ subroutine  psb_lc_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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ja(k))
-      end if
-      if (cpy) then 
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= czero).or.(i == a%ja(k))) 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 a802775a..e3e7b42c 100644
--- a/base/serial/impl/psb_d_coo_impl.F90
+++ b/base/serial/impl/psb_d_coo_impl.F90
@@ -595,19 +595,13 @@ 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 = .true.
-    else
-      ! Always keep the diagonal, even if numerically zero
-      cpy = (a%ia(i) == a%ja(i))
-    end if
-    if (cpy) then
+    ! Always keep the diagonal, even if numerically zero
+    if ((a%val(i) /= dzero).or.(a%ia(i) == a%ja(i))) then
       j = j + 1
       a%val(j) = a%val(i)
       a%ia(j)  = a%ia(i)
@@ -5933,19 +5927,13 @@ subroutine  psb_ld_coo_clean_zeros(a, info)
   integer(psb_ipk_), intent(out) :: info
   !
   integer(psb_lpk_) :: 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 = .true.
-    else
-      ! Always keep the diagonal, even if numerically zero
-      cpy = (a%ia(i) == a%ja(i))
-    end if
-    if (cpy) then 
+    ! Always keep the diagonal, even if numerically zero
+    if ((a%val(i) /= dzero).or.(a%ia(i) == a%ja(i))) 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 61e2ad14..886add04 100644
--- a/base/serial/impl/psb_d_csc_impl.F90
+++ b/base/serial/impl/psb_d_csc_impl.F90
@@ -2412,7 +2412,6 @@ 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()
@@ -2422,13 +2421,8 @@ 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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ia(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= dzero).or.(i == a%ia(k))) then
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
         j = j + 1
@@ -4320,7 +4314,6 @@ subroutine  psb_ld_csc_clean_zeros(a, info)
   !
   integer(psb_lpk_) :: i, j, k, nc
   integer(psb_lpk_), allocatable :: ilcp(:)
-  logical :: cpy
   
   info = 0
   call a%sync()
@@ -4330,13 +4323,8 @@ subroutine  psb_ld_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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ia(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= dzero).or.(i == a%ia(k))) 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 56ba8c63..f5891870 100644
--- a/base/serial/impl/psb_d_csr_impl.F90
+++ b/base/serial/impl/psb_d_csr_impl.F90
@@ -3633,7 +3633,6 @@ 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()
@@ -3643,13 +3642,8 @@ 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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ja(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= dzero).or.(i == a%ja(k))) then
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1
@@ -6559,7 +6553,6 @@ subroutine  psb_ld_csr_clean_zeros(a, info)
   !
   integer(psb_lpk_) :: i, j, k, nr
   integer(psb_lpk_), allocatable :: ilrp(:)
-  logical :: cpy
   
   info = 0
   call a%sync()
@@ -6569,13 +6562,8 @@ subroutine  psb_ld_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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ja(k))
-      end if
-      if (cpy) then 
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= dzero).or.(i == a%ja(k))) 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 0979ff85..023cde51 100644
--- a/base/serial/impl/psb_s_coo_impl.F90
+++ b/base/serial/impl/psb_s_coo_impl.F90
@@ -595,19 +595,13 @@ 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 = .true.
-    else
-      ! Always keep the diagonal, even if numerically zero
-      cpy = (a%ia(i) == a%ja(i))
-    end if
-    if (cpy) then
+    ! Always keep the diagonal, even if numerically zero
+    if ((a%val(i) /= szero).or.(a%ia(i) == a%ja(i))) then
       j = j + 1
       a%val(j) = a%val(i)
       a%ia(j)  = a%ia(i)
@@ -5933,19 +5927,13 @@ subroutine  psb_ls_coo_clean_zeros(a, info)
   integer(psb_ipk_), intent(out) :: info
   !
   integer(psb_lpk_) :: 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 = .true.
-    else
-      ! Always keep the diagonal, even if numerically zero
-      cpy = (a%ia(i) == a%ja(i))
-    end if
-    if (cpy) then 
+    ! Always keep the diagonal, even if numerically zero
+    if ((a%val(i) /= szero).or.(a%ia(i) == a%ja(i))) 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 ca41d705..3bb47d95 100644
--- a/base/serial/impl/psb_s_csc_impl.F90
+++ b/base/serial/impl/psb_s_csc_impl.F90
@@ -2412,7 +2412,6 @@ 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()
@@ -2422,13 +2421,8 @@ 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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ia(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= szero).or.(i == a%ia(k))) then
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
         j = j + 1
@@ -4320,7 +4314,6 @@ subroutine  psb_ls_csc_clean_zeros(a, info)
   !
   integer(psb_lpk_) :: i, j, k, nc
   integer(psb_lpk_), allocatable :: ilcp(:)
-  logical :: cpy
   
   info = 0
   call a%sync()
@@ -4330,13 +4323,8 @@ subroutine  psb_ls_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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ia(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= szero).or.(i == a%ia(k))) 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 323601ca..0a166b0c 100644
--- a/base/serial/impl/psb_s_csr_impl.F90
+++ b/base/serial/impl/psb_s_csr_impl.F90
@@ -3633,7 +3633,6 @@ 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()
@@ -3643,13 +3642,8 @@ 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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ja(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= szero).or.(i == a%ja(k))) then
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1
@@ -6559,7 +6553,6 @@ subroutine  psb_ls_csr_clean_zeros(a, info)
   !
   integer(psb_lpk_) :: i, j, k, nr
   integer(psb_lpk_), allocatable :: ilrp(:)
-  logical :: cpy
   
   info = 0
   call a%sync()
@@ -6569,13 +6562,8 @@ subroutine  psb_ls_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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ja(k))
-      end if
-      if (cpy) then 
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= szero).or.(i == a%ja(k))) 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 74c3f2cb..7dfceb06 100644
--- a/base/serial/impl/psb_z_coo_impl.F90
+++ b/base/serial/impl/psb_z_coo_impl.F90
@@ -595,19 +595,13 @@ 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 = .true.
-    else
-      ! Always keep the diagonal, even if numerically zero
-      cpy = (a%ia(i) == a%ja(i))
-    end if
-    if (cpy) then
+    ! Always keep the diagonal, even if numerically zero
+    if ((a%val(i) /= zzero).or.(a%ia(i) == a%ja(i))) then
       j = j + 1
       a%val(j) = a%val(i)
       a%ia(j)  = a%ia(i)
@@ -5933,19 +5927,13 @@ subroutine  psb_lz_coo_clean_zeros(a, info)
   integer(psb_ipk_), intent(out) :: info
   !
   integer(psb_lpk_) :: 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 = .true.
-    else
-      ! Always keep the diagonal, even if numerically zero
-      cpy = (a%ia(i) == a%ja(i))
-    end if
-    if (cpy) then 
+    ! Always keep the diagonal, even if numerically zero
+    if ((a%val(i) /= zzero).or.(a%ia(i) == a%ja(i))) 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 7ceff47f..32be36af 100644
--- a/base/serial/impl/psb_z_csc_impl.F90
+++ b/base/serial/impl/psb_z_csc_impl.F90
@@ -2412,7 +2412,6 @@ 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()
@@ -2422,13 +2421,8 @@ 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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ia(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= zzero).or.(i == a%ia(k))) then
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
         j = j + 1
@@ -4320,7 +4314,6 @@ subroutine  psb_lz_csc_clean_zeros(a, info)
   !
   integer(psb_lpk_) :: i, j, k, nc
   integer(psb_lpk_), allocatable :: ilcp(:)
-  logical :: cpy
   
   info = 0
   call a%sync()
@@ -4330,13 +4323,8 @@ subroutine  psb_lz_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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ia(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= zzero).or.(i == a%ia(k))) 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 54659def..e2ddf0d7 100644
--- a/base/serial/impl/psb_z_csr_impl.F90
+++ b/base/serial/impl/psb_z_csr_impl.F90
@@ -3633,7 +3633,6 @@ 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()
@@ -3643,13 +3642,8 @@ 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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ja(k))
-      end if
-      if (cpy) then
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= zzero).or.(i == a%ja(k))) then
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1
@@ -6559,7 +6553,6 @@ subroutine  psb_lz_csr_clean_zeros(a, info)
   !
   integer(psb_lpk_) :: i, j, k, nr
   integer(psb_lpk_), allocatable :: ilrp(:)
-  logical :: cpy
   
   info = 0
   call a%sync()
@@ -6569,13 +6562,8 @@ subroutine  psb_lz_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 = .true.
-      else
-        ! Always keep the diagonal, even if numerically zero
-        cpy = (i == a%ja(k))
-      end if
-      if (cpy) then 
+      ! Always keep the diagonal, even if numerically zero
+      if ((a%val(k) /= zzero).or.(i == a%ja(k))) then 
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1