1
\$\begingroup\$

I wrote this Fortran module and I need its subroutines to run as fast as possible. My hope is to get some advice to optimize this code, and after having the serial version running as fast as it can run I will use MPI to make it parallel.

Module smearingmodule
IMPLICIT NONE
contains
subroutine smearing(i, nMCd, nPeriodT, nPeriods, DF, theta, wm) bind(C, name="smearing_")
    use, intrinsic                         :: iso_c_binding, only : c_double, c_int
    integer(c_int), intent(in)                                 :: i, nMCd, nPeriods
    integer(c_int), intent(in), dimension(nPeriods)            :: nPeriodT
    real(c_double), intent(in), dimension(i,5)        :: DF
    real(c_double), intent(in), dimension(i,nMCd)     :: theta
    real(c_double), dimension(nMCd, nPeriods), intent(out)  :: wm
    integer                                             :: d

    do d=1,nMCd
        CALL subsmearing(d, i, nMCd, nPeriodT, nPeriods, DF, theta, wm(d,:))
    end do

end subroutine smearing


subroutine subsmearing(d, i, nMCd, nPeriodT, nPeriods, DF, theta, wm)
    integer, intent(in)                                 :: d, i, nMCd, nPeriods
    integer, intent(in), dimension(nPeriods)            :: nPeriodT
    double precision, intent(in), dimension(i,5)        :: DF
    double precision, intent(in), dimension(i,nMCd)     :: theta
    double precision, dimension(nPeriods), intent(out)  :: wm
    double precision, allocatable, dimension(:,:)       :: out
    double precision, dimension(i,8)                    :: bigDF
    double precision, dimension(i)                      :: epredC, epredT, diff, A1, A2
    double precision, dimension(i,i)                    :: B1, B2
    integer                                             :: j, Period, jj
    double precision                                    :: sumWeights

    diff = 0.0d0
    wm = 0.0d0
    epredC = exp(DF(:,4) + (theta(:,d) * DF(:,1)))
    epredT = exp(DF(:,4) - (theta(:,d) * (1-DF(:,1))))

    do j=1,i
        B1(:,j)=epredT(j)
        B2(:,j)=epredC(j)
    end do

    A1 = matmul(DF(:,5), B1)
    A2 = matmul(DF(:,5), B2)

    diff = (A1-A2)/DBLE(i)

    bigDF(:,1:5) = DF
    bigDF(:,6) = epredC
    bigDF(:,7) = epredT
    bigDF(:,8) = diff

    do jj =1, nPeriods
        ALLOCATE(out(nPeriodT(jj),7))
        CALL subsetPeriod(bigDF, i, 8, nPeriodT(jj), jj, out) !Filters data
        sumWeights = sum(out(:,2))
        do j=1, nPeriodT(jj)
            wm(jj) = wm(jj) + (out(j,7)*out(j,2)/sumWeights)
        end do
        DEALLOCATE(out)
    end do

end subroutine subsmearing

subroutine subsetPeriod(A, rowA, colA, rowB, Period, B)
    integer, intent(in)                                    :: rowA, colA, rowB, Period
    double precision, dimension(rowA, colA), intent(in)    :: A
    double precision, dimension(rowB,colA-1), intent(out)  :: B
    integer                                                :: i, pos

    pos = 1
    do i = 1, size(A,1)
        if(A(i,2)==Period)then
            B(pos,1) = A(i,1)
            B(pos,2:) = A(i,3:)
            pos = pos+1
        end if
    end do
end subroutine subsetPeriod

end module smearingmodule

DF is going to be millions by 5, theta millions by 4000. That is why I want to optimize the code and make it parallel.

Ultimately, I will be calling this code from R.


Link in case you are curious about smearing

\$\endgroup\$
2
\$\begingroup\$

I don't know, I might be totally mistaken, but from parsing your code and what it supposedly does, I think it could be simplified quite a lot. When thinking about performance, one of the main resources to treat carefully is the memory, you should minimize its access where ever possible.

One point, which totally confused me here, was the matrix multiplication you set up with a matrix of just duplicated entries. I also would change some commenting and don't see the point in splitting this routine up. Especially your last routine should be covered by a forall loop. Thus, I'd suggest something like this:

!> Using a comment to describe the purpose of this module.
!! Deploying an underscore to separate words
module smearing_module
  use, intrinsic :: iso_c_binding, only: c_double, c_int

  implicit none

  private ! Make everything private

  ! Use a kind parameter to use for your reals (instead of the
  ! old declaration with "double precision").
  integer, parameter :: rk = c_double

  ! Only export things, that should be used from outside.
  public :: smearing, rk


contains ! I intend everything in modules, except the contains.


  !> Employ some description comment here.
  !! Do you really need to mimic the name-mangling, or would it be okay
  !! to leave the underscore out in the exported name?
  subroutine smearing(i, nMCd, nPeriodT, nPeriods, DF, theta, wm) bind(C, name="smearing_")
      !> It helps the reader to provide a short explanation on the parameters.
      !! Thus, I think, parameters should one line each with some comment.
      integer(kind=c_int), intent(in) :: i

      !> Though, this is a little more verbose, it can help a lot,
      !! and you can make use of doxygen to automatically create
      !! an API documentation for your code.
      integer(kind=c_int), intent(in) :: nMCd

      !> I also like to add the kind keyword.
      integer(kind=c_int), intent(in) :: nPeriods

      !> I prefer to have the array declarations after the variable,
      !! especially if it applies just to a single variable:
      integer(kind=c_int), intent(in) :: nPeriodT(nPeriods)

      !> With the array definition on the right side, single declarations
      !! get shorter, as the dimension keyword can be omitted.
      real(kind=rk), intent(in) :: DF(i,5)

      !> It also puts an emphasis on the other attributes.
      real(kind=rk), intent(in) :: theta(i,nMCd)

      !> Also the reader just needs to look at the variable name, to
      !! directly see that this is supposed to be an array. He does
      !! not need to parse the list of attributes for this. 
      real(kind=rk), intent(out) :: wm(nMCd, nPeriods)

      ! local variables:
      ! I like to separate local variables a little from the parameters.
      integer       :: d, j, Period, jj
      real(kind=rk) :: sumWeights
      real(kind=rk) :: diff(i)
      real(kind=rk) :: dsum

      ! Summation independent of d. Thus, can be pulled out here.
      dsum = sum(DF(:,5))/real(i,kind=rk)

      wm = 0.0_rk

      do d=1,nMCd
          ! Only the difference between epredT and epredC is needed later on.
          ! Therefore, it is more efficient to just compute this, saves memory
          ! and bandwidth.
          do j=1,i
              diff(j) = exp(DF(j,4) - (theta(j,d) * (1.0_rk-DF(j,1)))) &
                  & - exp(DF(j,4) + (theta(j,d) * DF(j,1)))
          end do

          ! Essentially, your matrix multiplication boiled down to this scalar
          ! multiplication here. It might even be better, to pull this into
          ! the loop above to avoid another iteration over j.
          diff = diff * dsum

          do jj=1,nPeriods
              sumWeights = 0.0_rk
              ! The forall construct lets you do the masking, you put into
              ! a separate routine:
              forall(j=1:i, DF(j,2) == jj)
                  sumWeights = sumWeights + DF(j,3)
                  wm(d,jj) = wm(d,jj) + diff(j)*DF(j,3)
              end forall
              wm(d,jj) = wm(d,jj) / sumWeights
          end do
      end do

  end subroutine smearing

end module smearing_module

Depending on nMCd it might also be better to change the ordering of the loops or storing diff for all d beforehand and things like that. Note, how much less memory needs to be touched in my suggestion. Especially you can avoid copying all the data from DF to bigDF.

\$\endgroup\$
2
\$\begingroup\$

I added some comments to your code. I haven't tested anything, but that's the first things I would do to optimize your code. Under Linux, if you compile with g, you can use perf top to see live where your code takes time. Good luck!

Module smearingmodule
IMPLICIT NONE
contains
subroutine smearing(i, nMCd, nPeriodT, nPeriods, DF, theta, wm) bind(C, name="smearing_")
    use, intrinsic                         :: iso_c_binding, only : c_double, c_int
    integer(c_int), intent(in)                                 :: i, nMCd, nPeriods
    integer(c_int), intent(in), dimension(nPeriods)            :: nPeriodT
    real(c_double), intent(in), dimension(i,5)        :: DF
    real(c_double), intent(in), dimension(i,nMCd)     :: theta
    real(c_double), dimension(nMCd, nPeriods), intent(out)  :: wm
    integer                                             :: d

!   To use more efficiently vectorization in the subsmearing sburoutine,
!   DF(:,4), theta(:,d) and DF(:,1) need to be 256-bit aligned.
!   This can be done by :
!   1) Setting the size of the 1st dimension of DF and theta to a multiple of 4
!      double-precision elements (padding at the end of each column)
!   2) Forcing the arrays to start at the proper alignment. With the
!      ifort compiler, this can be done after the allocation using
!      !DIR$ ATTRIBUTES ALIGN : 32 :: DF, theta
!      
    do d=1,nMCd
    ! * Now you have to pass the physical sizes of the arrays DF and 
    !   theta because they don''t correspond to the logical size
    !   because of the padding
    ! * If you can, use wf dimensioned as (nPeriods,nMCd) to avoid
    !   the implicit matrix transposition in (d,:)
        CALL subsmearing(d, i, nMCd, nPeriodT, nPeriods, DF, theta, wm(d,:))
    end do

end subroutine smearing


subroutine subsmearing(d, i, nMCd, nPeriodT, nPeriods, DF, theta, wm)
    integer, intent(in)                                 :: d, i, nMCd, nPeriods
    integer, intent(in), dimension(nPeriods)            :: nPeriodT
    double precision, intent(in), dimension(i,5)        :: DF
    double precision, intent(in), dimension(i,nMCd)     :: theta
    double precision, dimension(nPeriods), intent(out)  :: wm
    double precision, allocatable, dimension(:,:)       :: out
    double precision, dimension(i,8)                    :: bigDF
    double precision, dimension(i)                      :: epredC, epredT, diff, A1, A2
    double precision, dimension(i,i)                    :: B1, B2
    integer                                             :: j, Period, jj
    double precision                                    :: sumWeights

!   In general, put explicitly the boundaries in your array syntax. This
!   makes the code more explicit and in case the physical dimensions
!   don''t correspond to the logical dimensions (because of padding for
!   instance), your code will be more robust.

    diff = 0.0d0
    wm = 0.0d0

!   1) Use 1.d0 insetad of 1 to avoid conversion
!   2) Remove array syntax here so that DF streams only once through the CPU
!     epredC = 0.d0
!     epredT = 0.d0
!     do k=1, i ! and not size(DF,1) because of padding
!      epredC = epredC + exp(DF(k,4) + (theta(k,d) * DF(k,1)))
!      epredT = epredT + exp(DF(k,4) - (theta(k,d) * (1.d0-DF(k,1))))
!     end do
    epredC = exp(DF(:,4) + (theta(:,d) * DF(:,1)))
    epredT = exp(DF(:,4) - (theta(:,d) * (1-DF(:,1))))


    do j=1,i
        B1(:,j)=epredT(j)
        B2(:,j)=epredC(j)
    end do

    ! Here, replace matmul by Lapack matrix multiply (dgemm).
    ! Matmul is good only for small matrices.
    ! Otherwise, put the boundaries of DF as DF(1:i,5)
    call dgemm('N','N',
    A1 = matmul(DF(:,5), B1)
    A2 = matmul(DF(:,5), B2)

    diff = (A1-A2)/DBLE(i)

    ! It would be better to dimension bigDF as (8,i) here for 
    ! the subsetPeriod subroutine
    ! In that case, you would need:
    ! do k=1,i
    !   bigDF(1:5,k) = DF(k)
    !   bigDF(  6,k) = epredC
    !   bigDF(  7,k) = epredT
    !   bigDF(  8,k) = diff
    ! end do
    ! If bigDF is properly aligned, you will have a very efficient store
    ! here

    bigDF(:,1:5) = DF
    bigDF(:,6) = epredC
    bigDF(:,7) = epredT
    bigDF(:,8) = diff

    ! 1) Put the allocate and deallocate statements outside of the
    !     loop.
    ! 2) Dimension the array as ( 8, maxval(nPeriodT) ). 8 is better
    !    than 7 for vectorization, and the transposed version is
    !    preferable for the called loop
    ! 3) Compute sumWeights_inv = 1.d0/sumWeights outside of the loop
    !    and multiply by sumWeights_inv in the loop : a division is a
    !    evry expensive operation
    ! 4) sumWeights could be computed on the fly for free in
    !    subsetPeriod. You would need to pass it as an argument
    do jj =1, nPeriods
        ALLOCATE(out(nPeriodT(jj),7))
        CALL subsetPeriod(bigDF, i, 8, nPeriodT(jj), jj, out) !Filters data
        sumWeights = sum(out(:,2))
        do j=1, nPeriodT(jj)
            ! Multiply by inverse
            wm(jj) = wm(jj) + (out(j,7)*out(j,2)/sumWeights)
        end do
        DEALLOCATE(out)
    end do

end subroutine subsmearing

subroutine subsetPeriod(A, rowA, colA, rowB, Period, B)
    integer, intent(in)                                    :: rowA, colA, rowB, Period
    double precision, dimension(rowA, colA), intent(in)    :: A
    double precision, dimension(rowB,colA-1), intent(out)  :: B
    integer                                                :: i, pos

    pos = 1
    do i = 1, size(A,1)
        if(A(i,2)==Period)then
            B(pos,1) = A(i,1)
    ! Transposing bigDF and out is important for better
    ! data access in this line
            B(pos,2:) = A(i,3:)
            pos = pos+1
        end if
    end do
end subroutine subsetPeriod

end module smearingmodule
\$\endgroup\$
  • \$\begingroup\$ Thanks @Anthony Scemama! I'm trying to follow your suggestions one by one and already broke my code when trying to implement the first one, transposing BigDF. I edited my question to reflect this. What am I doing wrong? Thanks again \$\endgroup\$ – Ignacio Jul 29 '15 at 13:02
  • \$\begingroup\$ What you may and may not do after receiving answers. I've rolled back Rev 4 → 2. \$\endgroup\$ – 200_success Jul 29 '15 at 15:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.