Optimize fdderivs: skip redundant 2nd-order work in 4th-order overlap

AMSS_NCKU_source/fdderivs_c.C

@@ -141,12 +141,26 @@ void fdderivs(const int ex[3],
     const int j4_hi = ex2 - 3;
     const int k4_hi = ex3 - 3;
 
+    /*
+     * Strategy A:
+     * Avoid redundant work in overlap of 2nd/4th-order regions.
+     * Only compute 2nd-order on shell points that are NOT overwritten by
+     * the 4th-order pass.
+     */
+    const int has4 = (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi);
+
     if (i2_lo <= i2_hi && j2_lo <= j2_hi && k2_lo <= k2_hi) {
         for (int k0 = k2_lo; k0 <= k2_hi; ++k0) {
             const int kF = k0 + 1;
             for (int j0 = j2_lo; j0 <= j2_hi; ++j0) {
                 const int jF = j0 + 1;
                 for (int i0 = i2_lo; i0 <= i2_hi; ++i0) {
+                    if (has4 &&
+                        i0 >= i4_lo && i0 <= i4_hi &&
+                        j0 >= j4_lo && j0 <= j4_hi &&
+                        k0 >= k4_lo && k0 <= k4_hi) {
+                        continue;
+                    }
                     const int iF = i0 + 1;
                     const size_t p = idx_ex(i0, j0, k0, ex);
 
@@ -193,7 +207,7 @@ void fdderivs(const int ex[3],
         }
     }
 
-    if (i4_lo <= i4_hi && j4_lo <= j4_hi && k4_lo <= k4_hi) {
+    if (has4) {
         for (int k0 = k4_lo; k0 <= k4_hi; ++k0) {
             const int kF = k0 + 1;
             for (int j0 = j4_lo; j0 <= j4_hi; ++j0) {

AMSS_NCKU_source/prolongrestrict_cell.f90

@@ -1955,7 +1955,11 @@
   integer::maxcx,maxcy,maxcz
 
   real*8,dimension(3) :: CD,FD
-  
+  real*8 :: tmp_yz(extc(1), 6)      ! 存储整条 X 线上 6 个 Y 轴偏置的 Z 向插值结果
+  real*8 :: tmp_xyz_line(extc(1))   ! 存储整条 X 线上完成 Y 向融合后的结果
+  real*8 :: v1, v2, v3, v4, v5, v6
+  integer :: ic, jc, kc, ix_offset,ix,iy,iz
+  real*8 :: res_line
   if(wei.ne.3)then
      write(*,*)"prolongrestrict.f90::prolong3: this routine only surport 3 dimension"
      write(*,*)"dim = ",wei
@@ -2070,114 +2074,74 @@
      
 !~~~~~~> prolongation start...
  do k = kmino, kmaxo
-   do j = jmino,jmaxo
-    do i = imino,imaxo
-       cxI(1) = cix(i)
-       cxI(2) = ciy(j)
-       cxI(3) = ciz(k)
-       px = pix(i)
-       py = piy(j)
      pz = piz(k)
+     kc = ciz(k)
+
+     do j = jmino, jmaxo
+        py = piy(j)
+        jc = ciy(j)
+
+! --- 步骤 1 & 2 融合：分段处理 X 轴，提升 Cache 命中率 ---
+        ! 我们将 ii 循环逻辑重组，减少对 funcc 的跨行重复访问
+        do ii = 1, extc(1)
+           ! 1. 先做 Z 方向的 6 条线插值（针对当前的 ii 和当前的 6 个 iy）
+           ! 我们直接在这里把 Y 方向的加权也做了，省去 tmp_yz 数组
+           ! 这样 funcc 的数据读进来后立即完成所有维度的贡献，不再写回内存
+           
+           res_line = 0.0d0
+           do jj = 1, 6
+              iy = jc - 3 + jj
+              ! 这一行代码是核心：一次性完成 Z 插值并加上 Y 的权重
+              ! 编译器会把 WC(jj, py) 存在寄存器里
+              res_line = res_line + WC(jj, py) * ( &
+                         WC(1, pz) * funcc(ii, iy, kc-2) + &
+                         WC(2, pz) * funcc(ii, iy, kc-1) + &
+                         WC(3, pz) * funcc(ii, iy, kc  ) + &
+                         WC(4, pz) * funcc(ii, iy, kc+1) + &
+                         WC(5, pz) * funcc(ii, iy, kc+2) + &
+                         WC(6, pz) * funcc(ii, iy, kc+3) )
+           end do
+           tmp_xyz_line(ii) = res_line
+        end do
+
+
+
 #if 0
-       if(ii/2*2==ii)then
-         if(jj/2*2==jj)then
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)=  C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           endif
-         else
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)= C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)=  C1*tmp1(1)+C2*tmp1(2)+C3*tmp1(3)+C4*tmp1(4)+C5*tmp1(5)+C6*tmp1(6)
-           endif
-         endif
-       else
-         if(jj/2*2==jj)then
-           if(kk/2*2==kk)then               
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C1*tmp2(:,1)+C2*tmp2(:,2)+C3*tmp2(:,3)+C4*tmp2(:,4)+C5*tmp2(:,5)+C6*tmp2(:,6)
-             funf(i,j,k)=  C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           endif
-         else
-           if(kk/2*2==kk)then
-             tmp2= C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)= C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           else
-             tmp2= C6*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-                   C5*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-                   C4*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-                   C3*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-                   C2*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-                   C1*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
-             tmp1= C6*tmp2(:,1)+C5*tmp2(:,2)+C4*tmp2(:,3)+C3*tmp2(:,4)+C2*tmp2(:,5)+C1*tmp2(:,6)
-             funf(i,j,k)=  C6*tmp1(1)+C5*tmp1(2)+C4*tmp1(3)+C3*tmp1(4)+C2*tmp1(5)+C1*tmp1(6)
-           endif
-         endif
-       endif
-#else 
-       tmp2= WC(1,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-2)+&
-             WC(2,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)-1)+&
-             WC(3,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)  )+&
-             WC(4,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+1)+&
-             WC(5,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+2)+&
-             WC(6,pz)*funcc(cxI(1)-2:cxI(1)+3,cxI(2)-2:cxI(2)+3,cxI(3)+3)
+        ! 1. 【降维：Z 向】对当前 (j,k) 相关的 6 条 Y 偏置线进行 Z 向插值
+        ! 结果存入 tmp_yz(x_index, y_offset)
+        do jj = 1, 6
+           iy = jc - 3 + jj
+           do ii = 1, extc(1)
+              tmp_yz(ii, jj) = WC(1,pz)*funcc(ii, iy, kc-2) + &
+                               WC(2,pz)*funcc(ii, iy, kc-1) + &
+                               WC(3,pz)*funcc(ii, iy, kc  ) + &
+                               WC(4,pz)*funcc(ii, iy, kc+1) + &
+                               WC(5,pz)*funcc(ii, iy, kc+2) + &
+                               WC(6,pz)*funcc(ii, iy, kc+3)
+           end do
+        end do
 
-       tmp1= WC(1,py)*tmp2(:,1)+WC(2,py)*tmp2(:,2)+WC(3,py)*tmp2(:,3)+&
-             WC(4,py)*tmp2(:,4)+WC(5,py)*tmp2(:,5)+WC(6,py)*tmp2(:,6)
-
-       funf(i,j,k)= WC(1,px)*tmp1(1)+WC(2,px)*tmp1(2)+WC(3,px)*tmp1(3)+&
-                    WC(4,px)*tmp1(4)+WC(5,px)*tmp1(5)+WC(6,px)*tmp1(6)
+        ! 2. 【降维：Y 向】将 Z 向结果合并，得到整条 X 轴线上的 Y-Z 融合值
+        do ii = 1, extc(1)
+           tmp_xyz_line(ii) = WC(1,py)*tmp_yz(ii, 1) + WC(2,py)*tmp_yz(ii, 2) + &
+                              WC(3,py)*tmp_yz(ii, 3) + WC(4,py)*tmp_yz(ii, 4) + &
+                              WC(5,py)*tmp_yz(ii, 5) + WC(6,py)*tmp_yz(ii, 6)
+        end do
 #endif
+        ! 3. 【降维：X 向】最后在最内层只处理 X 方向的 6 点加权
+        ! 此时每个点的计算量从原来的 200+ 次乘法降到了仅 6 次
+        do i = imino, imaxo
+           px = pix(i)
+           ic = cix(i)
+           
+           ! 直接从预计算好的 line 中读取连续的 6 个点
+           ! ic-2 到 ic+3 对应原始 6 点算子
+           funf(i,j,k) = WC(1,px)*tmp_xyz_line(ic-2) + &
+                         WC(2,px)*tmp_xyz_line(ic-1) + &
+                         WC(3,px)*tmp_xyz_line(ic  ) + &
+                         WC(4,px)*tmp_xyz_line(ic+1) + &
+                         WC(5,px)*tmp_xyz_line(ic+2) + &
+                         WC(6,px)*tmp_xyz_line(ic+3)
         end do
      end do
   end do