From f7db723fdb8fd6e294df12007c551e8454c173b5 Mon Sep 17 00:00:00 2001
From: Stack-1 <simone.staccone@virgilio.it>
Date: Sun, 12 Jul 2026 17:07:54 +0200
Subject: [PATCH] [ADD] Added loop unrolling using 4 accumulators in csr SpMV

---
 base/serial/impl/psb_d_csr_impl.F90 | 94 +++++++++++++++++++++++------
 1 file changed, 76 insertions(+), 18 deletions(-)

diff --git a/base/serial/impl/psb_d_csr_impl.F90 b/base/serial/impl/psb_d_csr_impl.F90
index 763bd4f22..d74bd6b99 100644
--- a/base/serial/impl/psb_d_csr_impl.F90
+++ b/base/serial/impl/psb_d_csr_impl.F90
@@ -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