fdurastante/psblas3
1
0
Fork 0
mirror of https://github.com/sfilippone/psblas3.git synced 2026-10-16 10:52:08 +00:00
Code Issues Packages Projects Releases Wiki Activity

[ADD] Added loop unrolling using 4 accumulators in csr SpMV

Browse Source
This commit is contained in:
Stack-1
2026-07-12 17:07:54 +02:00
parent d7f3bf536d
commit f7db723fdb
1 changed files with 76 additions and 18 deletions
Download Patch File Download Diff File Expand all files Collapse all files
base/serial/impl/psb_d_csr_impl.F90
+76 -18
View File
@@ -127,6 +127,14 @@ contains
integer(psb_ipk_) :: i,j,ir
real(psb_dpk_) :: acc
integer(psb_ipk_) :: ir0, ir1, navg
real(psb_dpk_) :: acc0, acc1, acc2, acc3
! Whole-matrix dispatch on the mean row length: long-row matrices use
! a 4-accumulator unrolled loop that breaks the FMA latency chain;
! short-row matrices keep the original loop untouched. The branch sits
! outside the row loop because a per-row branch mispredicts on
! irregular matrices whose rows straddle the threshold.
integer(psb_ipk_), parameter :: ithresh = 16
if (alpha == dzero) then
if (beta == dzero) then
@@ -146,18 +154,46 @@ contains
if ((.not.tra).and.(.not.ctra)) then
if (m > 0) then
navg = (irp(m+1)-irp(1))/m
else
navg = 0
end if
if (beta == dzero) then
if (alpha == done) then
!$omp parallel do private(i,j, acc) schedule(static)
do i=1,m
acc = dzero
!$omp simd reduction(+:acc)
do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j))
enddo
y(i) = acc
end do
if (navg >= ithresh) then
!$omp parallel do private(i,j,ir0,ir1,acc0,acc1,acc2,acc3) schedule(static)
do i=1,m
ir0 = irp(i)
ir1 = irp(i+1)-1
acc0 = dzero
acc1 = dzero
acc2 = dzero
acc3 = dzero
do j=ir0, ir1-3, 4
acc0 = acc0 + val(j) * x(ja(j))
acc1 = acc1 + val(j+1) * x(ja(j+1))
acc2 = acc2 + val(j+2) * x(ja(j+2))
acc3 = acc3 + val(j+3) * x(ja(j+3))
enddo
do j=ir0 + 4*((ir1-ir0+1)/4), ir1
acc0 = acc0 + val(j) * x(ja(j))
enddo
y(i) = (acc0+acc1) + (acc2+acc3)
end do
else
!$omp parallel do private(i,j, acc) schedule(static)
do i=1,m
acc = dzero
!$omp simd reduction(+:acc)
do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j))
enddo
y(i) = acc
end do
end if
else if (alpha == -done) then
@@ -189,15 +225,37 @@ contains
else if (beta == done) then
if (alpha == done) then
!$omp parallel do private(i,j,acc)
do i=1,m
acc = dzero
!$omp simd reduction(+:acc)
do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j))
enddo
y(i) = y(i) + acc
end do
if (navg >= ithresh) then
!$omp parallel do private(i,j,ir0,ir1,acc0,acc1,acc2,acc3) schedule(static)
do i=1,m
ir0 = irp(i)
ir1 = irp(i+1)-1
acc0 = dzero
acc1 = dzero
acc2 = dzero
acc3 = dzero
do j=ir0, ir1-3, 4
acc0 = acc0 + val(j) * x(ja(j))
acc1 = acc1 + val(j+1) * x(ja(j+1))
acc2 = acc2 + val(j+2) * x(ja(j+2))
acc3 = acc3 + val(j+3) * x(ja(j+3))
enddo
do j=ir0 + 4*((ir1-ir0+1)/4), ir1
acc0 = acc0 + val(j) * x(ja(j))
enddo
y(i) = y(i) + (acc0+acc1) + (acc2+acc3)
end do
else
!$omp parallel do private(i,j,acc)
do i=1,m
acc = dzero
!$omp simd reduction(+:acc)
do j=irp(i), irp(i+1)-1
acc = acc + val(j) * x(ja(j))
enddo
y(i) = y(i) + acc
end do
end if
else if (alpha == -done) then