修改prolong

AMSS_NCKU_source/prolongrestrict_cell.f90

@@ -1912,7 +1912,7 @@
 !   63/8192*f_6 - 495/8192*f_5 + 1155/4096*f_4
 !--------------------------------------------------------------------------
 #if 1
-   subroutine prolong3(wei,llbc,uubc,extc,func,&
+  subroutine prolong3(wei,llbc,uubc,extc,func,&
                       llbf,uubf,extf,funf,&
                       llbp,uubp,SoA,Symmetry)
   implicit none
@@ -1944,14 +1944,6 @@
   integer, dimension(extf(2)) :: piy
   integer, dimension(extf(3)) :: piz
 
-! ~~~> 在子程序开头新增声明 <~~~
-  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, iy, ix_offset
-  real*8 :: v_m2, v_m1, v_0, v_p1, v_p2, v_p3
-  real*8 :: res_line
-! ~~~> 声明end <~~~
   real*8, parameter :: C1=7.7d1/8.192d3,C2=-6.93d2/8.192d3,C3=3.465d3/4.096d3
   real*8, parameter :: C6=6.3d1/8.192d3,C5=-4.95d2/8.192d3,C4=1.155d3/4.096d3
   real*8, dimension(6,2), parameter :: WC = reshape((/&
@@ -1963,6 +1955,10 @@
   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
 
   if(wei.ne.3)then
      write(*,*)"prolongrestrict.f90::prolong3: this routine only surport 3 dimension"
@@ -2077,9 +2073,7 @@
   call symmetry_bd(3,extc,func,funcc,SoA)
      
 !~~~~~~> prolongation start...
-  ! ~~~> 循环体开始 <~~~
+ do k = kmino, kmaxo
-#if 0
-  do k = kmino, kmaxo
      pz = piz(k)
      kc = ciz(k)
 
@@ -2108,38 +2102,8 @@
                               WC(5,py)*tmp_yz(ii, 5) + WC(6,py)*tmp_yz(ii, 6)
         end do
 
-        i = imino
+        ! 3. 【降维：X 向】最后在最内层只处理 X 方向的 6 点加权
-        do while (i < imaxo)
+        ! 此时每个点的计算量从原来的 200+ 次乘法降到了仅 6 次
-           ! 针对 i 
-           ic = cix(i)
-           v_m2 = tmp_xyz_line(ic-2); v_m1 = tmp_xyz_line(ic-1)
-           v_0  = tmp_xyz_line(ic  ); v_p1 = tmp_xyz_line(ic+1)
-           v_p2 = tmp_xyz_line(ic+2); v_p3 = tmp_xyz_line(ic+3)
-
-           ! 这里的 WC(1:6, pix(i)) 被直接硬编码
-           ! 注意：根据你的 pix 定义，1 和 2 分别对应 C1..C6 和 C6..C1
-           if (pix(i) == 2) then
-              funf(i, j, k) = C6*v_m2 + C5*v_m1 + C4*v_0 + C3*v_p1 + C2*v_p2 + C1*v_p3
-              funf(i+1, j, k) = C1*v_m2 + C2*v_m1 + C3*v_0 + C4*v_p1 + C5*v_p2 + C6*v_p3
-           else
-              funf(i, j, k) = C1*v_m2 + C2*v_m1 + C3*v_0 + C4*v_p1 + C5*v_p2 + C6*v_p3
-              funf(i+1, j, k) = C6*v_m2 + C5*v_m1 + C4*v_0 + C3*v_p1 + C2*v_p2 + C1*v_p3
-           end if
-           
-           i = i + 2
-        end do
-        
-        ! 处理可能存在的最后一个孤立点
-        if (i == imaxo) then
-           ic = cix(i); px = pix(i)
-           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 if
-
-
-
-#if 0
         do i = imino, imaxo
            px = pix(i)
            ic = cix(i)
@@ -2153,70 +2117,9 @@
                          WC(5,px)*tmp_xyz_line(ic+2) + &
                          WC(6,px)*tmp_xyz_line(ic+3)
         end do
-#endif
      end do
   end do
-#endif
 
-
-! ~~~> 修改后的核心循环块 <~~~
-  do k = kmino, kmaxo
-     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
-
-        ! --- 步骤 3：X 向插值（保持你现在的双点展开，这部分已经很快了） ---
-        i = imino
-        do while (i < imaxo)
-           ic = cix(i)
-           v_m2 = tmp_xyz_line(ic-2); v_m1 = tmp_xyz_line(ic-1)
-           v_0  = tmp_xyz_line(ic  ); v_p1 = tmp_xyz_line(ic+1)
-           v_p2 = tmp_xyz_line(ic+2); v_p3 = tmp_xyz_line(ic+3)
-
-           if (pix(i) == 2) then
-              funf(i,   j, k) = C6*v_m2 + C5*v_m1 + C4*v_0 + C3*v_p1 + C2*v_p2 + C1*v_p3
-              funf(i+1, j, k) = C1*v_m2 + C2*v_m1 + C3*v_0 + C4*v_p1 + C5*v_p2 + C6*v_p3
-           else
-              funf(i,   j, k) = C1*v_m2 + C2*v_m1 + C3*v_0 + C4*v_p1 + C5*v_p2 + C6*v_p3
-              funf(i+1, j, k) = C6*v_m2 + C5*v_m1 + C4*v_0 + C3*v_p1 + C2*v_p2 + C1*v_p3
-           end if
-           i = i + 2
-        end do
-        
-        if (i == imaxo) then
-           ic = cix(i); px = pix(i)
-           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 if
-     end do
-  end do
   return
 
   end subroutine prolong3