Compare commits
4 Commits
yx-prolong
...
hxh-omp
| Author | SHA1 | Date | |
|---|---|---|---|
| e11363e06e | |||
| f70e90f694 | |||
|
75dd5353b0 | ||
|
|
23a82d063b |
@@ -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) {
|
||||
|
||||
@@ -1958,9 +1958,8 @@
|
||||
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,jc_min,jc_max
|
||||
integer :: ic, jc, kc, ix_offset,ix,iy,iz
|
||||
real*8 :: res_line
|
||||
real*8 :: tmp_z_slab(extc(1), extc(2)) ! 分配在 k 循环外
|
||||
if(wei.ne.3)then
|
||||
write(*,*)"prolongrestrict.f90::prolong3: this routine only surport 3 dimension"
|
||||
write(*,*)"dim = ",wei
|
||||
@@ -2072,34 +2071,8 @@
|
||||
endif
|
||||
|
||||
call symmetry_bd(3,extc,func,funcc,SoA)
|
||||
! 对每个 k（pz, kc 固定）预计算 Z 向插值的 2D 切片
|
||||
jc_min = minval(ciy(jmino:jmaxo))
|
||||
jc_max = maxval(ciy(jmino:jmaxo))
|
||||
|
||||
do k = kmino, kmaxo
|
||||
pz = piz(k); kc = ciz(k)
|
||||
! --- Pass 1: Z 方向，只算一次 ---
|
||||
do iy = jc_min-3, jc_max+3 ! 仅需的 iy 范围
|
||||
do ii = imini-3, imaxi+3 ! 仅需的 ii 范围
|
||||
tmp_z_slab(ii, iy) = sum(WC(:,pz) * funcc(ii, iy, kc-2:kc+3))
|
||||
end do
|
||||
end do
|
||||
|
||||
do j = jmino, jmaxo
|
||||
py = piy(j); jc = ciy(j)
|
||||
! --- Pass 2: Y 方向 ---
|
||||
do ii = imini-3, imaxi+3
|
||||
tmp_xyz_line(ii) = sum(WC(:,py) * tmp_z_slab(ii, jc-2:jc+3))
|
||||
end do
|
||||
! --- Pass 3: X 方向 ---
|
||||
do i = imino, imaxo
|
||||
funf(i,j,k) = sum(WC(:,pix(i)) * tmp_xyz_line(cix(i)-2:cix(i)+3))
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
|
||||
!~~~~~~> prolongation start...
|
||||
#if 0
|
||||
do k = kmino, kmaxo
|
||||
pz = piz(k)
|
||||
kc = ciz(k)
|
||||
@@ -2133,7 +2106,28 @@ end do
|
||||
|
||||
|
||||
|
||||
#if 0
|
||||
! 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
|
||||
|
||||
! 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
|
||||
@@ -2151,7 +2145,7 @@ end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
#endif
|
||||
|
||||
return
|
||||
|
||||
end subroutine prolong3
|
||||
@@ -2350,11 +2344,7 @@ end do
|
||||
integer::imino,imaxo,jmino,jmaxo,kmino,kmaxo
|
||||
|
||||
real*8,dimension(3) :: CD,FD
|
||||
|
||||
real*8 :: tmp_xz_plane(extf(1), 6)
|
||||
real*8 :: tmp_x_line(extf(1))
|
||||
integer :: fi, fj, fk, ii, jj, kk
|
||||
|
||||
|
||||
if(wei.ne.3)then
|
||||
write(*,*)"prolongrestrict.f90::restrict3: this routine only surport 3 dimension"
|
||||
write(*,*)"dim = ",wei
|
||||
@@ -2436,56 +2426,6 @@ end do
|
||||
call symmetry_bd(2,extf,funf,funff,SoA)
|
||||
|
||||
!~~~~~~> restriction start...
|
||||
do k = kmino, kmaxo
|
||||
fk = 2*(k + lbc(3) - 1) - 1 - lbf(3) + 1
|
||||
|
||||
do j = jmino, jmaxo
|
||||
fj = 2*(j + lbc(2) - 1) - 1 - lbf(2) + 1
|
||||
|
||||
! 优化点 1: 显式展开 Z 方向计算，减少循环开销
|
||||
! 确保 ii 循环是最内层且连续访问
|
||||
!DIR$ VECTOR ALWAYS
|
||||
do ii = 1, extf(1)
|
||||
! 预计算当前 j 对应的 6 行在 Z 方向的压缩结果
|
||||
! 这里直接硬编码 jj 的偏移，彻底消除一层循环
|
||||
tmp_xz_plane(ii, 1) = C1*(funff(ii,fj-2,fk-2)+funff(ii,fj-2,fk+3)) + &
|
||||
C2*(funff(ii,fj-2,fk-1)+funff(ii,fj-2,fk+2)) + &
|
||||
C3*(funff(ii,fj-2,fk )+funff(ii,fj-2,fk+1))
|
||||
tmp_xz_plane(ii, 2) = C1*(funff(ii,fj-1,fk-2)+funff(ii,fj-1,fk+3)) + &
|
||||
C2*(funff(ii,fj-1,fk-1)+funff(ii,fj-1,fk+2)) + &
|
||||
C3*(funff(ii,fj-1,fk )+funff(ii,fj-1,fk+1))
|
||||
tmp_xz_plane(ii, 3) = C1*(funff(ii,fj ,fk-2)+funff(ii,fj ,fk+3)) + &
|
||||
C2*(funff(ii,fj ,fk-1)+funff(ii,fj ,fk+2)) + &
|
||||
C3*(funff(ii,fj ,fk )+funff(ii,fj ,fk+1))
|
||||
tmp_xz_plane(ii, 4) = C1*(funff(ii,fj+1,fk-2)+funff(ii,fj+1,fk+3)) + &
|
||||
C2*(funff(ii,fj+1,fk-1)+funff(ii,fj+1,fk+2)) + &
|
||||
C3*(funff(ii,fj+1,fk )+funff(ii,fj+1,fk+1))
|
||||
tmp_xz_plane(ii, 5) = C1*(funff(ii,fj+2,fk-2)+funff(ii,fj+2,fk+3)) + &
|
||||
C2*(funff(ii,fj+2,fk-1)+funff(ii,fj+2,fk+2)) + &
|
||||
C3*(funff(ii,fj+2,fk )+funff(ii,fj+2,fk+1))
|
||||
tmp_xz_plane(ii, 6) = C1*(funff(ii,fj+3,fk-2)+funff(ii,fj+3,fk+3)) + &
|
||||
C2*(funff(ii,fj+3,fk-1)+funff(ii,fj+3,fk+2)) + &
|
||||
C3*(funff(ii,fj+3,fk )+funff(ii,fj+3,fk+1))
|
||||
end do
|
||||
|
||||
! 优化点 2: 同样向量化 Y 方向压缩
|
||||
!DIR$ VECTOR ALWAYS
|
||||
do ii = 1, extf(1)
|
||||
tmp_x_line(ii) = C1*(tmp_xz_plane(ii, 1) + tmp_xz_plane(ii, 6)) + &
|
||||
C2*(tmp_xz_plane(ii, 2) + tmp_xz_plane(ii, 5)) + &
|
||||
C3*(tmp_xz_plane(ii, 3) + tmp_xz_plane(ii, 4))
|
||||
end do
|
||||
|
||||
! 优化点 3: 最终写入，利用已经缓存在 tmp_x_line 的数据
|
||||
do i = imino, imaxo
|
||||
fi = 2*(i + lbc(1) - 1) - 1 - lbf(1) + 1
|
||||
func(i, j, k) = C1*(tmp_x_line(fi-2) + tmp_x_line(fi+3)) + &
|
||||
C2*(tmp_x_line(fi-1) + tmp_x_line(fi+2)) + &
|
||||
C3*(tmp_x_line(fi ) + tmp_x_line(fi+1))
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
#if 0
|
||||
do k = kmino,kmaxo
|
||||
do j = jmino,jmaxo
|
||||
do i = imino,imaxo
|
||||
@@ -2509,7 +2449,7 @@ end do
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
#endif
|
||||
|
||||
return
|
||||
|
||||
end subroutine restrict3
|
||||
|
||||
Reference in New Issue
Block a user