ext: fix DIA/HDIA correctness, conversion and storage bugs

1. psb_?_hdia_csmv.f90: with beta /= 0 the kernel scaled y(i) instead of
   y(rdisp+i), so every hack after the first used the wrong element and
   the result was silently wrong. Verified against CSR: the error was
   1.08e+02 on a 1000x1000 tridiagonal test, now 3.6e-15.

2. psb_?_cp_dia_to_coo.f90: the file did not use psi_ext_util_mod, so the
   call to psi_?_xtr_coo_from_dia had no explicit interface and the
   optional rdisp argument was seen as present with a garbage value.
   Any conversion out of DIA crashed. Added the module, as the HDIA
   version already does.

3. psi_?_xtr_coo_from_dia.f90: the extraction wrote one COO entry per
   stored position, padding zeros included, while the caller sizes the
   buffer with get_nzeros(). Skip the zeros, which also keeps the
   round trip DIA -> COO -> DIA exact in the number of nonzeros.
   psb_?_cp_dia_to_coo.f90 now calls set_nzeros with the count that the
   extraction returns.

4. psb_?_cp_hdia_from_coo.f90: diaOffsets holds one integer per diagonal
   per hack, but was allocated with hacksize*iszd, so 31/32 of the array
   was never used. On a 40^3 pdegen matrix the format occupied 5,314,884
   bytes instead of 3,601,204 (-32%); at 500^3 it would waste about
   3.4 GB per process.

5. psb_?_cp_hdia_from_coo.f90: a%nzeros counted the padded positions
   rather than the nonzeros, so get_nzeros() returned 442,222 instead of
   438,400 on the same matrix, inflating the MFLOPS of every benchmark
   and breaking the comparison with the other formats.

ext/impl/psb_c_cp_dia_to_coo.f90

@@ -34,6 +34,7 @@ subroutine psb_c_cp_dia_to_coo(a,b,info)
   
   use psb_base_mod
   use psb_c_dia_mat_mod, psb_protect_name => psb_c_cp_dia_to_coo
+  use psi_ext_util_mod
   implicit none 
 
   class(psb_c_dia_sparse_mat), intent(in)    :: a
@@ -58,7 +59,7 @@ subroutine psb_c_cp_dia_to_coo(a,b,info)
        & size(a%data,1),size(a%data,2),&
        & a%data,a%offset,info)
   
-  call b%set_nzeros(nza)
+  call b%set_nzeros(nzd)
   call b%set_host()
   call b%fix(info)
 

ext/impl/psb_c_cp_hdia_from_coo.f90

@@ -136,7 +136,7 @@ contains
       write(*,*) 'Hackcount ',nhacks,' Allocation height ',iszd
       write(*,*) 'Hackoffsets ',a%hackOffsets(:)
     end if
-    if (info == psb_success_) call psb_realloc(hacksize*iszd,a%diaOffsets,info)
+    if (info == psb_success_) call psb_realloc(iszd,a%diaOffsets,info)
     if (info == psb_success_) call psb_realloc(hacksize*iszd,a%val,info)
     if (info /= psb_success_) return
     klast1 = 1 
@@ -163,9 +163,7 @@ contains
            & a%val((hacksize*hackfirst)+1:hacksize*hacknext),info,&
            & initdata=.true.,rdisp=(i-1))
           
-      call countnz(nr,nc,(i-1),hacksize,(hacknext-hackfirst),&
-           & a%diaOffsets(hackfirst+1:hacknext),nzout)
-      a%nzeros = a%nzeros + nzout
+      a%nzeros = a%nzeros + (klast1-kfirst)
       call cleand(nr,(hacknext-hackfirst),d,a%diaOffsets(hackfirst+1:hacknext))
       
     end do

ext/impl/psb_c_hdia_csmv.f90

@@ -138,7 +138,7 @@ contains
       enddo
     else
       do  i = 1, min(nrd,nr-rdisp)
-        y(rdisp+i) = beta*y(i)
+        y(rdisp+i) = beta*y(rdisp+i)
       end do
     endif
     do j=1, ncd

ext/impl/psb_d_cp_dia_to_coo.f90

@@ -34,6 +34,7 @@ subroutine psb_d_cp_dia_to_coo(a,b,info)
   
   use psb_base_mod
   use psb_d_dia_mat_mod, psb_protect_name => psb_d_cp_dia_to_coo
+  use psi_ext_util_mod
   implicit none 
 
   class(psb_d_dia_sparse_mat), intent(in)    :: a
@@ -58,7 +59,7 @@ subroutine psb_d_cp_dia_to_coo(a,b,info)
        & size(a%data,1),size(a%data,2),&
        & a%data,a%offset,info)
   
-  call b%set_nzeros(nza)
+  call b%set_nzeros(nzd)
   call b%set_host()
   call b%fix(info)
 

ext/impl/psb_d_cp_hdia_from_coo.f90

@@ -136,7 +136,7 @@ contains
       write(*,*) 'Hackcount ',nhacks,' Allocation height ',iszd
       write(*,*) 'Hackoffsets ',a%hackOffsets(:)
     end if
-    if (info == psb_success_) call psb_realloc(hacksize*iszd,a%diaOffsets,info)
+    if (info == psb_success_) call psb_realloc(iszd,a%diaOffsets,info)
     if (info == psb_success_) call psb_realloc(hacksize*iszd,a%val,info)
     if (info /= psb_success_) return
     klast1 = 1 
@@ -163,9 +163,7 @@ contains
            & a%val((hacksize*hackfirst)+1:hacksize*hacknext),info,&
            & initdata=.true.,rdisp=(i-1))
           
-      call countnz(nr,nc,(i-1),hacksize,(hacknext-hackfirst),&
-           & a%diaOffsets(hackfirst+1:hacknext),nzout)
-      a%nzeros = a%nzeros + nzout
+      a%nzeros = a%nzeros + (klast1-kfirst)
       call cleand(nr,(hacknext-hackfirst),d,a%diaOffsets(hackfirst+1:hacknext))
       
     end do

ext/impl/psb_d_hdia_csmv.f90

@@ -138,7 +138,7 @@ contains
       enddo
     else
       do  i = 1, min(nrd,nr-rdisp)
-        y(rdisp+i) = beta*y(i)
+        y(rdisp+i) = beta*y(rdisp+i)
       end do
     endif
     do j=1, ncd

ext/impl/psb_s_cp_dia_to_coo.f90

@@ -34,6 +34,7 @@ subroutine psb_s_cp_dia_to_coo(a,b,info)
   
   use psb_base_mod
   use psb_s_dia_mat_mod, psb_protect_name => psb_s_cp_dia_to_coo
+  use psi_ext_util_mod
   implicit none 
 
   class(psb_s_dia_sparse_mat), intent(in)    :: a
@@ -58,7 +59,7 @@ subroutine psb_s_cp_dia_to_coo(a,b,info)
        & size(a%data,1),size(a%data,2),&
        & a%data,a%offset,info)
   
-  call b%set_nzeros(nza)
+  call b%set_nzeros(nzd)
   call b%set_host()
   call b%fix(info)
 

ext/impl/psb_s_cp_hdia_from_coo.f90

@@ -136,7 +136,7 @@ contains
       write(*,*) 'Hackcount ',nhacks,' Allocation height ',iszd
       write(*,*) 'Hackoffsets ',a%hackOffsets(:)
     end if
-    if (info == psb_success_) call psb_realloc(hacksize*iszd,a%diaOffsets,info)
+    if (info == psb_success_) call psb_realloc(iszd,a%diaOffsets,info)
     if (info == psb_success_) call psb_realloc(hacksize*iszd,a%val,info)
     if (info /= psb_success_) return
     klast1 = 1 
@@ -163,9 +163,7 @@ contains
            & a%val((hacksize*hackfirst)+1:hacksize*hacknext),info,&
            & initdata=.true.,rdisp=(i-1))
           
-      call countnz(nr,nc,(i-1),hacksize,(hacknext-hackfirst),&
-           & a%diaOffsets(hackfirst+1:hacknext),nzout)
-      a%nzeros = a%nzeros + nzout
+      a%nzeros = a%nzeros + (klast1-kfirst)
       call cleand(nr,(hacknext-hackfirst),d,a%diaOffsets(hackfirst+1:hacknext))
       
     end do

ext/impl/psb_s_hdia_csmv.f90

@@ -138,7 +138,7 @@ contains
       enddo
     else
       do  i = 1, min(nrd,nr-rdisp)
-        y(rdisp+i) = beta*y(i)
+        y(rdisp+i) = beta*y(rdisp+i)
       end do
     endif
     do j=1, ncd

ext/impl/psb_z_cp_dia_to_coo.f90

@@ -34,6 +34,7 @@ subroutine psb_z_cp_dia_to_coo(a,b,info)
   
   use psb_base_mod
   use psb_z_dia_mat_mod, psb_protect_name => psb_z_cp_dia_to_coo
+  use psi_ext_util_mod
   implicit none 
 
   class(psb_z_dia_sparse_mat), intent(in)    :: a
@@ -58,7 +59,7 @@ subroutine psb_z_cp_dia_to_coo(a,b,info)
        & size(a%data,1),size(a%data,2),&
        & a%data,a%offset,info)
   
-  call b%set_nzeros(nza)
+  call b%set_nzeros(nzd)
   call b%set_host()
   call b%fix(info)
 

ext/impl/psb_z_cp_hdia_from_coo.f90

@@ -136,7 +136,7 @@ contains
       write(*,*) 'Hackcount ',nhacks,' Allocation height ',iszd
       write(*,*) 'Hackoffsets ',a%hackOffsets(:)
     end if
-    if (info == psb_success_) call psb_realloc(hacksize*iszd,a%diaOffsets,info)
+    if (info == psb_success_) call psb_realloc(iszd,a%diaOffsets,info)
     if (info == psb_success_) call psb_realloc(hacksize*iszd,a%val,info)
     if (info /= psb_success_) return
     klast1 = 1 
@@ -163,9 +163,7 @@ contains
            & a%val((hacksize*hackfirst)+1:hacksize*hacknext),info,&
            & initdata=.true.,rdisp=(i-1))
           
-      call countnz(nr,nc,(i-1),hacksize,(hacknext-hackfirst),&
-           & a%diaOffsets(hackfirst+1:hacknext),nzout)
-      a%nzeros = a%nzeros + nzout
+      a%nzeros = a%nzeros + (klast1-kfirst)
       call cleand(nr,(hacknext-hackfirst),d,a%diaOffsets(hackfirst+1:hacknext))
       
     end do

ext/impl/psb_z_hdia_csmv.f90

@@ -138,7 +138,7 @@ contains
       enddo
     else
       do  i = 1, min(nrd,nr-rdisp)
-        y(rdisp+i) = beta*y(i)
+        y(rdisp+i) = beta*y(rdisp+i)
       end do
     endif
     do j=1, ncd

ext/impl/psi_c_xtr_coo_from_dia.f90

@@ -69,10 +69,12 @@ subroutine psi_c_xtr_coo_from_dia(nr,nc,ia,ja,val,nz,nrd,ncd,data,offsets,info,r
       ir = i + rdisp_
       ic = i + rdisp_ + offsets(j)
       if (debug) write(0,*) ' Loop  I',i,ir,ic
-      nz = nz + 1
-      ia(nz) = ir
-      ja(nz) = ic 
-      val(nz) = data(i,j)
+      if (data(i,j) /= czero) then 
+        nz = nz + 1
+        ia(nz) = ir
+        ja(nz) = ic 
+        val(nz) = data(i,j)
+      end if
     enddo
   end do
 

ext/impl/psi_d_xtr_coo_from_dia.f90

@@ -69,10 +69,12 @@ subroutine psi_d_xtr_coo_from_dia(nr,nc,ia,ja,val,nz,nrd,ncd,data,offsets,info,r
       ir = i + rdisp_
       ic = i + rdisp_ + offsets(j)
       if (debug) write(0,*) ' Loop  I',i,ir,ic
-      nz = nz + 1
-      ia(nz) = ir
-      ja(nz) = ic 
-      val(nz) = data(i,j)
+      if (data(i,j) /= dzero) then 
+        nz = nz + 1
+        ia(nz) = ir
+        ja(nz) = ic 
+        val(nz) = data(i,j)
+      end if
     enddo
   end do
 

ext/impl/psi_s_xtr_coo_from_dia.f90

@@ -69,10 +69,12 @@ subroutine psi_s_xtr_coo_from_dia(nr,nc,ia,ja,val,nz,nrd,ncd,data,offsets,info,r
       ir = i + rdisp_
       ic = i + rdisp_ + offsets(j)
       if (debug) write(0,*) ' Loop  I',i,ir,ic
-      nz = nz + 1
-      ia(nz) = ir
-      ja(nz) = ic 
-      val(nz) = data(i,j)
+      if (data(i,j) /= szero) then 
+        nz = nz + 1
+        ia(nz) = ir
+        ja(nz) = ic 
+        val(nz) = data(i,j)
+      end if
     enddo
   end do
 

ext/impl/psi_z_xtr_coo_from_dia.f90

@@ -69,10 +69,12 @@ subroutine psi_z_xtr_coo_from_dia(nr,nc,ia,ja,val,nz,nrd,ncd,data,offsets,info,r
       ir = i + rdisp_
       ic = i + rdisp_ + offsets(j)
       if (debug) write(0,*) ' Loop  I',i,ir,ic
-      nz = nz + 1
-      ia(nz) = ir
-      ja(nz) = ic 
-      val(nz) = data(i,j)
+      if (data(i,j) /= zzero) then 
+        nz = nz + 1
+        ia(nz) = ir
+        ja(nz) = ic 
+        val(nz) = data(i,j)
+      end if
     enddo
   end do