修改prolong
This commit is contained in:
jaunatisblue
@@ -1943,14 +1943,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((/&
|
||||
@@ -1962,6 +1954,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"
|
||||
@@ -2076,8 +2072,6 @@
|
||||
call symmetry_bd(3,extc,func,funcc,SoA)
|
||||
|
||||
!~~~~~~> prolongation start...
|
||||
! ~~~> 循环体开始 <~~~
|
||||
#if 0
|
||||
do k = kmino, kmaxo
|
||||
pz = piz(k)
|
||||
kc = ciz(k)
|
||||
@@ -2107,38 +2101,8 @@
|
||||
WC(5,py)*tmp_yz(ii, 5) + WC(6,py)*tmp_yz(ii, 6)
|
||||
end do
|
||||
|
||||
i = imino
|
||||
do while (i < imaxo)
|
||||
! 针对 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
|
||||
! 3. 【降维：X 向】最后在最内层只处理 X 方向的 6 点加权
|
||||
! 此时每个点的计算量从原来的 200+ 次乘法降到了仅 6 次
|
||||
do i = imino, imaxo
|
||||
px = pix(i)
|
||||
ic = cix(i)
|
||||
@@ -2152,70 +2116,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
|
||||
|
||||
Reference in New Issue
Block a user