From cd0fa59647dfa07823ddfc8531c5fa060da3d927 Mon Sep 17 00:00:00 2001
From: Stack-1 <simone.staccone@virgilio.it>
Date: Thu, 6 Aug 2026 16:12:07 +0200
Subject: [PATCH] [UPDATE] Extended the 4-accumulator SpMV unrolling to s/c/z

Generated from template-psblas, where the unrolled csmv_inner kernel was
added to the type-generic X_csr_impl.P90. The d version produced by the
template is byte-identical to the hand-written one in f7db723f, so this
commit only carries s/c/z.

Same dispatch as d: mean row length >= 16 nonzeros/row uses four
accumulators to break the FMA latency chain, shorter rows keep the
original loop. Applies to the beta==0/alpha==1 and beta==1/alpha==1
branches of psb_?_csr_csmv_inner.

Checked with mpif90 -O3 -fopenmp on all three types; the unrolled and
reference loops agree to rounding on rows of every length mod 4,
including empty rows.
---
 base/serial/impl/psb_c_csr_impl.F90 | 94 +++++++++++++++++++++++------
 base/serial/impl/psb_s_csr_impl.F90 | 94 +++++++++++++++++++++++------
 base/serial/impl/psb_z_csr_impl.F90 | 94 +++++++++++++++++++++++------
 3 files changed, 228 insertions(+), 54 deletions(-)

diff --git a/base/serial/impl/psb_c_csr_impl.F90 b/base/serial/impl/psb_c_csr_impl.F90
index d0fe429d8..d53c402c8 100644
--- a/base/serial/impl/psb_c_csr_impl.F90
+++ b/base/serial/impl/psb_c_csr_impl.F90
@@ -127,6 +127,14 @@ contains
 
     integer(psb_ipk_) :: i,j,ir
     complex(psb_spk_) :: acc
+    integer(psb_ipk_) :: ir0, ir1, navg
+    complex(psb_spk_) :: 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 == czero) then
       if (beta == czero) 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 == czero) then
 
         if (alpha == cone) then
-          !$omp parallel do private(i,j, acc) schedule(static)
-          do i=1,m
-            acc  = czero
-            !$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 = czero
+              acc1 = czero
+              acc2 = czero
+              acc3 = czero
+              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  = czero
+              !$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 == -cone) then
 
@@ -189,15 +225,37 @@ contains
       else if (beta == cone) then
 
         if (alpha == cone) then
-          !$omp parallel do private(i,j,acc)
-          do i=1,m
-            acc  = czero
-            !$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 = czero
+              acc1 = czero
+              acc2 = czero
+              acc3 = czero
+              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  = czero
+              !$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 == -cone) then
 
diff --git a/base/serial/impl/psb_s_csr_impl.F90 b/base/serial/impl/psb_s_csr_impl.F90
index 92bdc0f02..4b295c7e1 100644
--- a/base/serial/impl/psb_s_csr_impl.F90
+++ b/base/serial/impl/psb_s_csr_impl.F90
@@ -127,6 +127,14 @@ contains
 
     integer(psb_ipk_) :: i,j,ir
     real(psb_spk_) :: acc
+    integer(psb_ipk_) :: ir0, ir1, navg
+    real(psb_spk_) :: 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 == szero) then
       if (beta == szero) 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 == szero) then
 
         if (alpha == sone) then
-          !$omp parallel do private(i,j, acc) schedule(static)
-          do i=1,m
-            acc  = szero
-            !$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 = szero
+              acc1 = szero
+              acc2 = szero
+              acc3 = szero
+              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  = szero
+              !$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 == -sone) then
 
@@ -189,15 +225,37 @@ contains
       else if (beta == sone) then
 
         if (alpha == sone) then
-          !$omp parallel do private(i,j,acc)
-          do i=1,m
-            acc  = szero
-            !$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 = szero
+              acc1 = szero
+              acc2 = szero
+              acc3 = szero
+              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  = szero
+              !$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 == -sone) then
 
diff --git a/base/serial/impl/psb_z_csr_impl.F90 b/base/serial/impl/psb_z_csr_impl.F90
index 17dbd4abc..1a6d8be04 100644
--- a/base/serial/impl/psb_z_csr_impl.F90
+++ b/base/serial/impl/psb_z_csr_impl.F90
@@ -127,6 +127,14 @@ contains
 
     integer(psb_ipk_) :: i,j,ir
     complex(psb_dpk_) :: acc
+    integer(psb_ipk_) :: ir0, ir1, navg
+    complex(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 == zzero) then
       if (beta == zzero) 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 == zzero) then
 
         if (alpha == zone) then
-          !$omp parallel do private(i,j, acc) schedule(static)
-          do i=1,m
-            acc  = zzero
-            !$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 = zzero
+              acc1 = zzero
+              acc2 = zzero
+              acc3 = zzero
+              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  = zzero
+              !$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 == -zone) then
 
@@ -189,15 +225,37 @@ contains
       else if (beta == zone) then
 
         if (alpha == zone) then
-          !$omp parallel do private(i,j,acc)
-          do i=1,m
-            acc  = zzero
-            !$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 = zzero
+              acc1 = zzero
+              acc2 = zzero
+              acc3 = zzero
+              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  = zzero
+              !$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 == -zone) then