Fix clean_zeros to always keep the diagonal

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