From 1c235f928189169a9318d9027fd78f33740149ff Mon Sep 17 00:00:00 2001
From: sfilippone <filippone.salvatore@gmail.com>
Date: Fri, 4 Oct 2024 10:28:52 +0200
Subject: [PATCH] Improve clean_zeros

---
 base/serial/impl/psb_c_coo_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_c_csc_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_c_csr_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_d_coo_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_d_csc_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_d_csr_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_s_coo_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_s_csc_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_s_csr_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_z_coo_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_z_csc_impl.F90 | 20 +++++++++++++++-----
 base/serial/impl/psb_z_csr_impl.F90 | 20 +++++++++++++++-----
 12 files changed, 180 insertions(+), 60 deletions(-)

diff --git a/base/serial/impl/psb_c_coo_impl.F90 b/base/serial/impl/psb_c_coo_impl.F90
index 2a6aaf13..5b681d22 100644
--- a/base/serial/impl/psb_c_coo_impl.F90
+++ b/base/serial/impl/psb_c_coo_impl.F90
@@ -601,9 +601,12 @@ subroutine  psb_c_coo_clean_zeros(a, info)
   nzin = a%get_nzeros()
   j = 0
   do i=1, nzin
-    cpy = (a%val(i) /= czero)
-    ! Always keep the diagonal, even if numerically zero
-    if (.not.cpy) cpy = (a%ia(i) == a%ja(i))
+    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
       j = j + 1
       a%val(j) = a%val(i)
@@ -5930,12 +5933,19 @@ 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
+    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 
       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 f8ab1025..190a4d5b 100644
--- a/base/serial/impl/psb_c_csc_impl.F90
+++ b/base/serial/impl/psb_c_csc_impl.F90
@@ -2422,9 +2422,12 @@ subroutine  psb_c_csc_clean_zeros(a, info)
   j        = a%icp(1)
   do i=1, nc
     do k = ilcp(i), ilcp(i+1) -1
-      cpy = (a%val(k) /= czero)
-      ! Always keep the diagonal, even if numerically zero
-      if (.not.cpy) cpy = (i == a%ia(k))
+      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
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
@@ -4317,7 +4320,8 @@ 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()
   nc   = a%get_ncols()
@@ -4326,7 +4330,13 @@ 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
+      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
         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 d3494f46..0db9f3fa 100644
--- a/base/serial/impl/psb_c_csr_impl.F90
+++ b/base/serial/impl/psb_c_csr_impl.F90
@@ -3643,9 +3643,12 @@ subroutine  psb_c_csr_clean_zeros(a, info)
   j        = a%irp(1)
   do i=1, nr
     do k = ilrp(i), ilrp(i+1) -1
-      cpy = (a%val(k) /= czero)
-      ! Always keep the diagonal, even if numerically zero
-      if (.not.cpy) cpy = (i == a%ja(k))
+      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
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
@@ -6556,7 +6559,8 @@ 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()
   nr   = a%get_nrows()
@@ -6565,7 +6569,13 @@ 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
+      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 
         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 15d4c4ad..a802775a 100644
--- a/base/serial/impl/psb_d_coo_impl.F90
+++ b/base/serial/impl/psb_d_coo_impl.F90
@@ -601,9 +601,12 @@ subroutine  psb_d_coo_clean_zeros(a, info)
   nzin = a%get_nzeros()
   j = 0
   do i=1, nzin
-    cpy = (a%val(i) /= dzero)
-    ! Always keep the diagonal, even if numerically zero
-    if (.not.cpy) cpy = (a%ia(i) == a%ja(i))
+    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
       j = j + 1
       a%val(j) = a%val(i)
@@ -5930,12 +5933,19 @@ 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
+    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 
       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 92a4c555..61e2ad14 100644
--- a/base/serial/impl/psb_d_csc_impl.F90
+++ b/base/serial/impl/psb_d_csc_impl.F90
@@ -2422,9 +2422,12 @@ subroutine  psb_d_csc_clean_zeros(a, info)
   j        = a%icp(1)
   do i=1, nc
     do k = ilcp(i), ilcp(i+1) -1
-      cpy = (a%val(k) /= dzero)
-      ! Always keep the diagonal, even if numerically zero
-      if (.not.cpy) cpy = (i == a%ia(k))
+      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
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
@@ -4317,7 +4320,8 @@ 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()
   nc   = a%get_ncols()
@@ -4326,7 +4330,13 @@ 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
+      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
         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 09dffe0e..56ba8c63 100644
--- a/base/serial/impl/psb_d_csr_impl.F90
+++ b/base/serial/impl/psb_d_csr_impl.F90
@@ -3643,9 +3643,12 @@ subroutine  psb_d_csr_clean_zeros(a, info)
   j        = a%irp(1)
   do i=1, nr
     do k = ilrp(i), ilrp(i+1) -1
-      cpy = (a%val(k) /= dzero)
-      ! Always keep the diagonal, even if numerically zero
-      if (.not.cpy) cpy = (i == a%ja(k))
+      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
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
@@ -6556,7 +6559,8 @@ 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()
   nr   = a%get_nrows()
@@ -6565,7 +6569,13 @@ 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
+      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 
         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 217b39c2..0979ff85 100644
--- a/base/serial/impl/psb_s_coo_impl.F90
+++ b/base/serial/impl/psb_s_coo_impl.F90
@@ -601,9 +601,12 @@ subroutine  psb_s_coo_clean_zeros(a, info)
   nzin = a%get_nzeros()
   j = 0
   do i=1, nzin
-    cpy = (a%val(i) /= szero)
-    ! Always keep the diagonal, even if numerically zero
-    if (.not.cpy) cpy = (a%ia(i) == a%ja(i))
+    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
       j = j + 1
       a%val(j) = a%val(i)
@@ -5930,12 +5933,19 @@ 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
+    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 
       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 ffa32f8f..ca41d705 100644
--- a/base/serial/impl/psb_s_csc_impl.F90
+++ b/base/serial/impl/psb_s_csc_impl.F90
@@ -2422,9 +2422,12 @@ subroutine  psb_s_csc_clean_zeros(a, info)
   j        = a%icp(1)
   do i=1, nc
     do k = ilcp(i), ilcp(i+1) -1
-      cpy = (a%val(k) /= szero)
-      ! Always keep the diagonal, even if numerically zero
-      if (.not.cpy) cpy = (i == a%ia(k))
+      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
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
@@ -4317,7 +4320,8 @@ 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()
   nc   = a%get_ncols()
@@ -4326,7 +4330,13 @@ 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
+      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
         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 5e18de7b..323601ca 100644
--- a/base/serial/impl/psb_s_csr_impl.F90
+++ b/base/serial/impl/psb_s_csr_impl.F90
@@ -3643,9 +3643,12 @@ subroutine  psb_s_csr_clean_zeros(a, info)
   j        = a%irp(1)
   do i=1, nr
     do k = ilrp(i), ilrp(i+1) -1
-      cpy = (a%val(k) /= szero)
-      ! Always keep the diagonal, even if numerically zero
-      if (.not.cpy) cpy = (i == a%ja(k))
+      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
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
@@ -6556,7 +6559,8 @@ 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()
   nr   = a%get_nrows()
@@ -6565,7 +6569,13 @@ 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
+      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 
         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 89950c4c..74c3f2cb 100644
--- a/base/serial/impl/psb_z_coo_impl.F90
+++ b/base/serial/impl/psb_z_coo_impl.F90
@@ -601,9 +601,12 @@ subroutine  psb_z_coo_clean_zeros(a, info)
   nzin = a%get_nzeros()
   j = 0
   do i=1, nzin
-    cpy = (a%val(i) /= zzero)
-    ! Always keep the diagonal, even if numerically zero
-    if (.not.cpy) cpy = (a%ia(i) == a%ja(i))
+    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
       j = j + 1
       a%val(j) = a%val(i)
@@ -5930,12 +5933,19 @@ 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
+    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 
       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 a37e04e5..7ceff47f 100644
--- a/base/serial/impl/psb_z_csc_impl.F90
+++ b/base/serial/impl/psb_z_csc_impl.F90
@@ -2422,9 +2422,12 @@ subroutine  psb_z_csc_clean_zeros(a, info)
   j        = a%icp(1)
   do i=1, nc
     do k = ilcp(i), ilcp(i+1) -1
-      cpy = (a%val(k) /= zzero)
-      ! Always keep the diagonal, even if numerically zero
-      if (.not.cpy) cpy = (i == a%ia(k))
+      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
         a%val(j) = a%val(k)
         a%ia(j)  = a%ia(k)
@@ -4317,7 +4320,8 @@ 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()
   nc   = a%get_ncols()
@@ -4326,7 +4330,13 @@ 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
+      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
         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 57d31f80..54659def 100644
--- a/base/serial/impl/psb_z_csr_impl.F90
+++ b/base/serial/impl/psb_z_csr_impl.F90
@@ -3643,9 +3643,12 @@ subroutine  psb_z_csr_clean_zeros(a, info)
   j        = a%irp(1)
   do i=1, nr
     do k = ilrp(i), ilrp(i+1) -1
-      cpy = (a%val(k) /= zzero)
-      ! Always keep the diagonal, even if numerically zero
-      if (.not.cpy) cpy = (i == a%ja(k))
+      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
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
@@ -6556,7 +6559,8 @@ 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()
   nr   = a%get_nrows()
@@ -6565,7 +6569,13 @@ 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
+      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 
         a%val(j) = a%val(k)
         a%ja(j)  = a%ja(k)
         j = j + 1