refactor: Merge all transpose functions.

Author: semi-h

src/backend.f90

@@ -29,11 +29,7 @@ module m_base_backend
       procedure(transeq_ders), deferred :: transeq_y
       procedure(transeq_ders), deferred :: transeq_z
       procedure(tds_solve), deferred :: tds_solve
-      procedure(trans_d2d), deferred :: trans_x2y
-      procedure(trans_d2d), deferred :: trans_x2z
-      procedure(trans_d2d), deferred :: trans_y2z
-      procedure(trans_d2d), deferred :: trans_z2y
-      procedure(trans_d2d), deferred :: trans_y2x
+      procedure(trans_d2d), deferred :: trans_d2d
       procedure(sum9into3), deferred :: sum_yzintox
       procedure(vecadd), deferred :: vecadd
       procedure(get_fields), deferred :: get_fields
@@ -82,7 +78,7 @@ end subroutine tds_solve
    end interface
 
    abstract interface
-      subroutine trans_d2d(self, u_, u)
+      subroutine trans_d2d(self, u_, u, direction)
          !! transposer subroutines are straightforward, they rearrange
          !! data into our specialist data structure so that regardless
          !! of the direction tridiagonal systems are solved efficiently
@@ -94,6 +90,7 @@ subroutine trans_d2d(self, u_, u)
          class(base_backend_t) :: self
          class(field_t), intent(inout) :: u_
          class(field_t), intent(in) :: u
+         integer, intent(in) :: direction
       end subroutine trans_d2d
    end interface
 

src/cuda/backend.f90

@@ -34,11 +34,7 @@ module m_cuda_backend
       procedure :: transeq_y => transeq_y_cuda
       procedure :: transeq_z => transeq_z_cuda
       procedure :: tds_solve => tds_solve_cuda
-      procedure :: trans_x2y => trans_x2y_cuda
-      procedure :: trans_x2z => trans_x2z_cuda
-      procedure :: trans_y2z => trans_y2z_cuda
-      procedure :: trans_z2y => trans_z2y_cuda
-      procedure :: trans_y2x => trans_y2x_cuda
+      procedure :: trans_d2d => trans_d2d_cuda
       procedure :: sum_yzintox => sum_yzintox_cuda
       procedure :: vecadd => vecadd_cuda
       procedure :: set_fields => set_fields_cuda
@@ -421,107 +417,50 @@ subroutine tds_solve_dist(self, du, u, dirps, tdsops, blocks, threads)
 
    end subroutine tds_solve_dist
 
-   subroutine trans_x2y_cuda(self, u_y, u)
+   subroutine trans_d2d_cuda(self, u_o, u_i, direction)
       implicit none
 
       class(cuda_backend_t) :: self
-      class(field_t), intent(inout) :: u_y
-      class(field_t), intent(in) :: u
-
-      real(dp), device, pointer, dimension(:, :, :) :: u_d, u_y_d
-      type(dim3) :: blocks, threads
-
-      select type(u); type is (cuda_field_t); u_d => u%data_d; end select
-      select type(u_y); type is (cuda_field_t); u_y_d => u_y%data_d; end select
-
-      blocks = dim3(self%nx_loc/SZ, self%nz_loc, self%ny_loc/SZ)
-      threads = dim3(SZ, SZ, 1)
-
-      call trans_x2y_k<<<blocks, threads>>>(u_y_d, u_d, self%nz_loc)
-
-   end subroutine trans_x2y_cuda
-
-   subroutine trans_x2z_cuda(self, u_z, u)
-      implicit none
-
-      class(cuda_backend_t) :: self
-      class(field_t), intent(inout) :: u_z
-      class(field_t), intent(in) :: u
+      class(field_t), intent(inout) :: u_o
+      class(field_t), intent(in) :: u_i
+      integer, intent(in) :: direction
 
-      real(dp), device, pointer, dimension(:, :, :) :: u_d, u_z_d
+      real(dp), device, pointer, dimension(:, :, :) :: u_o_d, u_i_d
       type(dim3) :: blocks, threads
 
-      select type(u); type is (cuda_field_t); u_d => u%data_d; end select
-      select type(u_z); type is (cuda_field_t); u_z_d => u_z%data_d; end select
-
-      blocks = dim3(self%nx_loc, self%ny_loc/SZ, 1)
-      threads = dim3(SZ, 1, 1)
-
-      call trans_x2z_k<<<blocks, threads>>>(u_z_d, u_d, self%nz_loc)
-
-   end subroutine trans_x2z_cuda
-
-   subroutine trans_y2z_cuda(self, u_z, u_y)
-      implicit none
-
-      class(cuda_backend_t) :: self
-      class(field_t), intent(inout) :: u_z
-      class(field_t), intent(in) :: u_y
-
-      real(dp), device, pointer, dimension(:, :, :) :: u_z_d, u_y_d
-      type(dim3) :: blocks, threads
-
-      select type(u_z); type is (cuda_field_t); u_z_d => u_z%data_d; end select
-      select type(u_y); type is (cuda_field_t); u_y_d => u_y%data_d; end select
-
-      blocks = dim3(self%nx_loc/SZ, self%ny_loc/SZ, self%nz_loc)
-      threads = dim3(SZ, SZ, 1)
-
-      call trans_y2z_k<<<blocks, threads>>>(u_z_d, u_y_d, &
-                                            self%nx_loc, self%nz_loc)
-
-   end subroutine trans_y2z_cuda
-
-   subroutine trans_z2y_cuda(self, u_y, u_z)
-      implicit none
-
-      class(cuda_backend_t) :: self
-      class(field_t), intent(inout) :: u_y
-      class(field_t), intent(in) :: u_z
-
-      real(dp), device, pointer, dimension(:, :, :) :: u_y_d, u_z_d
-      type(dim3) :: blocks, threads
-
-      select type(u_y); type is (cuda_field_t); u_y_d => u_y%data_d; end select
-      select type(u_z); type is (cuda_field_t); u_z_d => u_z%data_d; end select
-
-      blocks = dim3(self%nx_loc/SZ, self%ny_loc/SZ, self%nz_loc)
-      threads = dim3(SZ, SZ, 1)
-
-      call trans_z2y_k<<<blocks, threads>>>(u_y_d, u_z_d, &
-                                            self%nx_loc, self%nz_loc)
-
-   end subroutine trans_z2y_cuda
-
-   subroutine trans_y2x_cuda(self, u_x, u_y)
-      implicit none
-
-      class(cuda_backend_t) :: self
-      class(field_t), intent(inout) :: u_x
-      class(field_t), intent(in) :: u_y
-
-      real(dp), device, pointer, dimension(:, :, :) :: u_x_d, u_y_d
-      type(dim3) :: blocks, threads
-
-      select type(u_x); type is (cuda_field_t); u_x_d => u_x%data_d; end select
-      select type(u_y); type is (cuda_field_t); u_y_d => u_y%data_d; end select
-
-      blocks = dim3(self%nx_loc/SZ, self%ny_loc/SZ, self%nz_loc)
-      threads = dim3(SZ, SZ, 1)
-
-      call trans_y2x_k<<<blocks, threads>>>(u_x_d, u_y_d, self%nz_loc)
+      select type(u_o); type is (cuda_field_t); u_o_d => u_o%data_d; end select
+      select type(u_i); type is (cuda_field_t); u_i_d => u_i%data_d; end select
+
+      select case (direction)
+      case (12) ! x2y
+         blocks = dim3(self%nx_loc/SZ, self%nz_loc, self%ny_loc/SZ)
+         threads = dim3(SZ, SZ, 1)
+         call trans_x2y_k<<<blocks, threads>>>(u_o_d, u_i_d, self%nz_loc)
+      case (13) ! x2z
+         blocks = dim3(self%nx_loc, self%ny_loc/SZ, 1)
+         threads = dim3(SZ, 1, 1)
+         call trans_x2z_k<<<blocks, threads>>>(u_o_d, u_i_d, self%nz_loc)
+      case (21) ! y2x
+         blocks = dim3(self%nx_loc/SZ, self%ny_loc/SZ, self%nz_loc)
+         threads = dim3(SZ, SZ, 1)
+         call trans_y2x_k<<<blocks, threads>>>(u_o_d, u_i_d, self%nz_loc)
+      case (23) ! y2z
+         blocks = dim3(self%nx_loc/SZ, self%ny_loc/SZ, self%nz_loc)
+         threads = dim3(SZ, SZ, 1)
+         call trans_y2z_k<<<blocks, threads>>>(u_o_d, u_i_d, &
+                                               self%nx_loc, self%nz_loc)
+      case (32) ! z2y
+         blocks = dim3(self%nx_loc/SZ, self%ny_loc/SZ, self%nz_loc)
+         threads = dim3(SZ, SZ, 1)
+
+         call trans_z2y_k<<<blocks, threads>>>(u_o_d, u_i_d, &
+                                               self%nx_loc, self%nz_loc)
+      case default
+         print *, 'Transpose direction is undefined.'
+         stop
+      end select
 
-   end subroutine trans_y2x_cuda
+   end subroutine trans_d2d_cuda
 
    subroutine sum_yzintox_cuda(self, du, dv, dw, &
                                du_y, dv_y, dw_y, du_z, dv_z, dw_z)

src/omp/backend.f90

@@ -24,11 +24,7 @@ module m_omp_backend
       procedure :: transeq_y => transeq_y_omp
       procedure :: transeq_z => transeq_z_omp
       procedure :: tds_solve => tds_solve_omp
-      procedure :: trans_x2y => trans_x2y_omp
-      procedure :: trans_x2z => trans_x2z_omp
-      procedure :: trans_y2z => trans_y2z_omp
-      procedure :: trans_z2y => trans_z2y_omp
-      procedure :: trans_y2x => trans_y2x_omp
+      procedure :: trans_d2d => trans_d2d_omp
       procedure :: sum_yzintox => sum_yzintox_omp
       procedure :: vecadd => vecadd_omp
       procedure :: set_fields => set_fields_omp
@@ -164,50 +160,15 @@ subroutine tds_solve_omp(self, du, u, dirps, tdsops)
 
    end subroutine tds_solve_omp
 
-   subroutine trans_x2y_omp(self, u_, u)
+   subroutine trans_d2d_omp(self, u_, u, direction)
       implicit none
 
       class(omp_backend_t) :: self
       class(field_t), intent(inout) :: u_
       class(field_t), intent(in) :: u
+      integer, intent(in) :: direction
 
-   end subroutine trans_x2y_omp
-
-   subroutine trans_x2z_omp(self, u_, u)
-      implicit none
-
-      class(omp_backend_t) :: self
-      class(field_t), intent(inout) :: u_
-      class(field_t), intent(in) :: u
-
-   end subroutine trans_x2z_omp
-
-   subroutine trans_y2z_omp(self, u_, u)
-      implicit none
-
-      class(omp_backend_t) :: self
-      class(field_t), intent(inout) :: u_
-      class(field_t), intent(in) :: u
-
-   end subroutine trans_y2z_omp
-
-   subroutine trans_z2y_omp(self, u_, u)
-      implicit none
-
-      class(omp_backend_t) :: self
-      class(field_t), intent(inout) :: u_
-      class(field_t), intent(in) :: u
-
-   end subroutine trans_z2y_omp
-
-   subroutine trans_y2x_omp(self, u_, u)
-      implicit none
-
-      class(omp_backend_t) :: self
-      class(field_t), intent(inout) :: u_
-      class(field_t), intent(in) :: u
-
-   end subroutine trans_y2x_omp
+   end subroutine trans_d2d_omp
 
    subroutine sum_yzintox_omp(self, du, dv, dw, &
                                du_y, dv_y, dw_y, du_z, dv_z, dw_z)

src/solver.f90

@@ -185,9 +185,9 @@ subroutine transeq(self, du, dv, dw, u, v, w)
       dw_y => self%backend%allocator%get_block()
 
       ! reorder data from x orientation to y orientation
-      call self%backend%trans_x2y(u_y, u)
-      call self%backend%trans_x2y(v_y, v)
-      call self%backend%trans_x2y(w_y, w)
+      call self%backend%trans_d2d(u_y, u, 12)
+      call self%backend%trans_d2d(v_y, v, 12)
+      call self%backend%trans_d2d(w_y, w, 12)
 
       ! similar to the x direction, obtain derivatives in y.
       call self%backend%transeq_y(du_y, dv_y, dw_y, u_y, v_y, w_y, self%ydirps)
@@ -209,9 +209,9 @@ subroutine transeq(self, du, dv, dw, u, v, w)
       dw_z => self%backend%allocator%get_block()
 
       ! reorder from x to z
-      call self%backend%trans_x2z(u_z, u)
-      call self%backend%trans_x2z(v_z, v)
-      call self%backend%trans_x2z(w_z, w)
+      call self%backend%trans_d2d(u_z, u, 13)
+      call self%backend%trans_d2d(v_z, v, 13)
+      call self%backend%trans_d2d(w_z, w, 13)
 
       ! get the derivatives in z
       call self%backend%transeq_z(du_z, dv_z, dw_z, u_z, v_z, w_z, self%zdirps)
@@ -267,9 +267,9 @@ subroutine divergence(self, div_u, u, v, w)
       w_y => self%backend%allocator%get_block()
 
       ! reorder data from x orientation to y orientation
-      call self%backend%trans_x2y(u_y, du_x)
-      call self%backend%trans_x2y(v_y, dv_x)
-      call self%backend%trans_x2y(w_y, dw_x)
+      call self%backend%trans_d2d(u_y, du_x, 12)
+      call self%backend%trans_d2d(v_y, dv_x, 12)
+      call self%backend%trans_d2d(w_y, dw_x, 12)
 
       call self%backend%allocator%release_block(du_x)
       call self%backend%allocator%release_block(dv_x)
@@ -303,8 +303,8 @@ subroutine divergence(self, div_u, u, v, w)
       call self%backend%vecadd(1._dp, dw_y, 1._dp, dv_y)
 
       ! reorder from y to z
-      call self%backend%trans_y2z(u_z, du_y)
-      call self%backend%trans_y2z(w_z, dw_y)
+      call self%backend%trans_d2d(u_z, du_y, 23)
+      call self%backend%trans_d2d(w_z, dw_y, 23)
 
       ! release all the unnecessary blocks.
       call self%backend%allocator%release_block(du_y)
@@ -358,8 +358,8 @@ subroutine gradient(self, dpdx, dpdy, dpdz, pressure)
       dpdz_sxy_y => self%backend%allocator%get_block()
 
       ! reorder data from z orientation to y orientation
-      call self%backend%trans_z2y(p_sxy_y, p_sxy_z)
-      call self%backend%trans_z2y(dpdz_sxy_y, dpdz_sxy_z)
+      call self%backend%trans_d2d(p_sxy_y, p_sxy_z, 32)
+      call self%backend%trans_d2d(dpdz_sxy_y, dpdz_sxy_z, 32)
 
       call self%backend%allocator%release_block(p_sxy_z)
       call self%backend%allocator%release_block(dpdz_sxy_z)
@@ -386,9 +386,9 @@ subroutine gradient(self, dpdx, dpdy, dpdz, pressure)
       dpdz_sx_x => self%backend%allocator%get_block()
 
       ! reorder from y to x
-      call self%backend%trans_y2x(p_sx_x, p_sx_y)
-      call self%backend%trans_y2x(dpdy_sx_x, dpdy_sx_y)
-      call self%backend%trans_y2x(dpdz_sx_x, dpdz_sx_y)
+      call self%backend%trans_d2d(p_sx_x, p_sx_y, 21)
+      call self%backend%trans_d2d(dpdy_sx_x, dpdy_sx_y, 21)
+      call self%backend%trans_d2d(dpdz_sx_x, dpdz_sx_y, 21)
 
       ! release all the y directional fields.
       call self%backend%allocator%release_block(p_sx_y)