From 241b078e72be0510748fa878792abc2a01991d3d Mon Sep 17 00:00:00 2001 From: Stack-1 Date: Wed, 30 Sep 2026 18:28:08 +0200 Subject: [PATCH] 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 | 3 ++- ext/impl/psb_c_cp_hdia_from_coo.f90 | 6 ++---- ext/impl/psb_c_hdia_csmv.f90 | 2 +- ext/impl/psb_d_cp_dia_to_coo.f90 | 3 ++- ext/impl/psb_d_cp_hdia_from_coo.f90 | 6 ++---- ext/impl/psb_d_hdia_csmv.f90 | 2 +- ext/impl/psb_s_cp_dia_to_coo.f90 | 3 ++- ext/impl/psb_s_cp_hdia_from_coo.f90 | 6 ++---- ext/impl/psb_s_hdia_csmv.f90 | 2 +- ext/impl/psb_z_cp_dia_to_coo.f90 | 3 ++- ext/impl/psb_z_cp_hdia_from_coo.f90 | 6 ++---- ext/impl/psb_z_hdia_csmv.f90 | 2 +- ext/impl/psi_c_xtr_coo_from_dia.f90 | 10 ++++++---- ext/impl/psi_d_xtr_coo_from_dia.f90 | 10 ++++++---- ext/impl/psi_s_xtr_coo_from_dia.f90 | 10 ++++++---- ext/impl/psi_z_xtr_coo_from_dia.f90 | 10 ++++++---- 16 files changed, 44 insertions(+), 40 deletions(-) diff --git a/ext/impl/psb_c_cp_dia_to_coo.f90 b/ext/impl/psb_c_cp_dia_to_coo.f90 index 8d2093221..d668700a9 100644 --- a/ext/impl/psb_c_cp_dia_to_coo.f90 +++ b/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) diff --git a/ext/impl/psb_c_cp_hdia_from_coo.f90 b/ext/impl/psb_c_cp_hdia_from_coo.f90 index e787e70a5..8244aec34 100644 --- a/ext/impl/psb_c_cp_hdia_from_coo.f90 +++ b/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 diff --git a/ext/impl/psb_c_hdia_csmv.f90 b/ext/impl/psb_c_hdia_csmv.f90 index c486cbebf..116bd59e3 100644 --- a/ext/impl/psb_c_hdia_csmv.f90 +++ b/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 diff --git a/ext/impl/psb_d_cp_dia_to_coo.f90 b/ext/impl/psb_d_cp_dia_to_coo.f90 index bace6af91..9c98a9638 100644 --- a/ext/impl/psb_d_cp_dia_to_coo.f90 +++ b/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) diff --git a/ext/impl/psb_d_cp_hdia_from_coo.f90 b/ext/impl/psb_d_cp_hdia_from_coo.f90 index 4d9440089..2140d4d0f 100644 --- a/ext/impl/psb_d_cp_hdia_from_coo.f90 +++ b/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 diff --git a/ext/impl/psb_d_hdia_csmv.f90 b/ext/impl/psb_d_hdia_csmv.f90 index 4944dff82..8715a9a08 100644 --- a/ext/impl/psb_d_hdia_csmv.f90 +++ b/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 diff --git a/ext/impl/psb_s_cp_dia_to_coo.f90 b/ext/impl/psb_s_cp_dia_to_coo.f90 index 4c103287d..030bdb7e4 100644 --- a/ext/impl/psb_s_cp_dia_to_coo.f90 +++ b/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) diff --git a/ext/impl/psb_s_cp_hdia_from_coo.f90 b/ext/impl/psb_s_cp_hdia_from_coo.f90 index 5cd7353c8..bfb143581 100644 --- a/ext/impl/psb_s_cp_hdia_from_coo.f90 +++ b/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 diff --git a/ext/impl/psb_s_hdia_csmv.f90 b/ext/impl/psb_s_hdia_csmv.f90 index 8f45a82b1..4638b1952 100644 --- a/ext/impl/psb_s_hdia_csmv.f90 +++ b/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 diff --git a/ext/impl/psb_z_cp_dia_to_coo.f90 b/ext/impl/psb_z_cp_dia_to_coo.f90 index 570aeca49..0cf2a06ab 100644 --- a/ext/impl/psb_z_cp_dia_to_coo.f90 +++ b/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) diff --git a/ext/impl/psb_z_cp_hdia_from_coo.f90 b/ext/impl/psb_z_cp_hdia_from_coo.f90 index 4f46bafb9..274633175 100644 --- a/ext/impl/psb_z_cp_hdia_from_coo.f90 +++ b/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 diff --git a/ext/impl/psb_z_hdia_csmv.f90 b/ext/impl/psb_z_hdia_csmv.f90 index 34e736f38..3feb41de3 100644 --- a/ext/impl/psb_z_hdia_csmv.f90 +++ b/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 diff --git a/ext/impl/psi_c_xtr_coo_from_dia.f90 b/ext/impl/psi_c_xtr_coo_from_dia.f90 index 59bc97e5a..7a40cb529 100644 --- a/ext/impl/psi_c_xtr_coo_from_dia.f90 +++ b/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 diff --git a/ext/impl/psi_d_xtr_coo_from_dia.f90 b/ext/impl/psi_d_xtr_coo_from_dia.f90 index 4e19ba783..7502906b5 100644 --- a/ext/impl/psi_d_xtr_coo_from_dia.f90 +++ b/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 diff --git a/ext/impl/psi_s_xtr_coo_from_dia.f90 b/ext/impl/psi_s_xtr_coo_from_dia.f90 index db58cd4c1..db0dab1e4 100644 --- a/ext/impl/psi_s_xtr_coo_from_dia.f90 +++ b/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 diff --git a/ext/impl/psi_z_xtr_coo_from_dia.f90 b/ext/impl/psi_z_xtr_coo_from_dia.f90 index 61cc65751..2af01a181 100644 --- a/ext/impl/psi_z_xtr_coo_from_dia.f90 +++ b/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