Smearing estimate with Fortran 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
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