Convolution in time domain

conv_lib.F90

! gfortran -ffree-line-length-none -fmax-identifier-length=63 -cpp -C -Wall -c conv_lib.F90

#ifdef __GFORTRAN__
#  define quote(x) "x"
#else
#  ifdef __INTEL_COMPILER
#    define quote(x) #x
#  else
#    define quote(x) #x
#  endif
#endif
#define WHERE_AM_I __FILE__, " ", __LINE__
#define RAISE(message) write(ERROR_UNIT, *) "RAISE: ", WHERE_AM_I, (message); error stop
#define RAISE_IF(isBad) \
   if(isBad)then; \
     RAISE(quote(isBad)); \
   end if
#define ASSERT(isOk) RAISE_IF(.not.(isOk))

module conv_lib
   use, intrinsic:: iso_fortran_env, only: INT8, INT16, INT32, INT64, REAL32, REAL64, REAL128, ERROR_UNIT

   implicit none

   private
   public:: conv

   interface conv
      module procedure convRealDim1KindREAL32
      module procedure convRealDim1KindREAL64
      module procedure convRealDim1KindREAL128
      module procedure convIntegerDim1KindINT8
      module procedure convIntegerDim1KindINT16
      module procedure convIntegerDim1KindINT32
      module procedure convIntegerDim1KindINT64
   end interface conv

   interface conv_short_with_long
      module procedure conv_short_with_longRealDim1KindREAL32
      module procedure conv_short_with_longRealDim1KindREAL64
      module procedure conv_short_with_longRealDim1KindREAL128
      module procedure conv_short_with_longIntegerDim1KindINT8
      module procedure conv_short_with_longIntegerDim1KindINT16
      module procedure conv_short_with_longIntegerDim1KindINT32
      module procedure conv_short_with_longIntegerDim1KindINT64
   end interface conv_short_with_long

contains

   ! zs = xs*ys (* is a conv operator)
   ! ASSUMPTIONS:
   ! - dX = dY = dZ
   ! CAUTIONS:
   ! - dX (dY) is not multipled.
   subroutine convRealDim1KindREAL32(xs, ys, zs)
      Real(kind=REAL32), dimension(:), intent(in):: xs, ys
      Real(kind=REAL32), dimension(:), intent(out):: zs

      if(size(xs, 1) >= size(ys, 1))then
         call conv_short_with_long(ys, xs, zs)
      else
         call conv_short_with_long(xs, ys, zs)
      end if
   end subroutine convRealDim1KindREAL32

   subroutine conv_short_with_longRealDim1KindREAL32(xs, ys, zs)
      Real(kind=REAL32), dimension(:), intent(in):: xs, ys
      Real(kind=REAL32), dimension(:), intent(out):: zs
      Integer:: nXs, nYs, nZs
      Integer:: i
      Integer:: iMax
      Logical:: isReachedEndOfZs

      nXs = size(xs, 1)
      nYs = size(ys, 1)
      nZs = size(zs, 1)
      ASSERT(nXs <= nYs)

      if(nZs <= 0) return
      zs = 0
      if(nXs <= 0) return ! this condition includes nYs <= 0 since nXs <= nYs
      if(nXs == 1)then
         zs = xs(1)*ys(1:nZs)
         return
      end if

      iMax = min(nXs - 1, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent (i = 1:iMax)
         zs(i) = dot_product(ys(1:i), xs(i:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nYs, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent(i = nXs:iMax)
         zs(i) = dot_product(ys(i-nXs+1:i), xs(nXs:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nXs + nYs, nZs)
      do concurrent(i = nYs + 1:iMax)
         zs(i) = dot_product(ys(i-nXs+1:nYs), xs(nXs:i-nYs+1:-1))
      end do
   end subroutine conv_short_with_longRealDim1KindREAL32

   subroutine convRealDim1KindREAL64(xs, ys, zs)
      Real(kind=REAL64), dimension(:), intent(in):: xs, ys
      Real(kind=REAL64), dimension(:), intent(out):: zs

      if(size(xs, 1) >= size(ys, 1))then
         call conv_short_with_long(ys, xs, zs)
      else
         call conv_short_with_long(xs, ys, zs)
      end if
   end subroutine convRealDim1KindREAL64

   subroutine conv_short_with_longRealDim1KindREAL64(xs, ys, zs)
      Real(kind=REAL64), dimension(:), intent(in):: xs, ys
      Real(kind=REAL64), dimension(:), intent(out):: zs
      Integer:: nXs, nYs, nZs
      Integer:: i
      Integer:: iMax
      Logical:: isReachedEndOfZs

      nXs = size(xs, 1)
      nYs = size(ys, 1)
      nZs = size(zs, 1)
      ASSERT(nXs <= nYs)

      if(nZs <= 0) return
      zs = 0
      if(nXs <= 0) return ! this condition includes nYs <= 0 since nXs <= nYs
      if(nXs == 1)then
         zs = xs(1)*ys(1:nZs)
         return
      end if

      iMax = min(nXs - 1, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent (i = 1:iMax)
         zs(i) = dot_product(ys(1:i), xs(i:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nYs, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent(i = nXs:iMax)
         zs(i) = dot_product(ys(i-nXs+1:i), xs(nXs:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nXs + nYs, nZs)
      do concurrent(i = nYs + 1:iMax)
         zs(i) = dot_product(ys(i-nXs+1:nYs), xs(nXs:i-nYs+1:-1))
      end do
   end subroutine conv_short_with_longRealDim1KindREAL64

   subroutine convRealDim1KindREAL128(xs, ys, zs)
      Real(kind=REAL128), dimension(:), intent(in):: xs, ys
      Real(kind=REAL128), dimension(:), intent(out):: zs

      if(size(xs, 1) >= size(ys, 1))then
         call conv_short_with_long(ys, xs, zs)
      else
         call conv_short_with_long(xs, ys, zs)
      end if
   end subroutine convRealDim1KindREAL128

   subroutine conv_short_with_longRealDim1KindREAL128(xs, ys, zs)
      Real(kind=REAL128), dimension(:), intent(in):: xs, ys
      Real(kind=REAL128), dimension(:), intent(out):: zs
      Integer:: nXs, nYs, nZs
      Integer:: i
      Integer:: iMax
      Logical:: isReachedEndOfZs

      nXs = size(xs, 1)
      nYs = size(ys, 1)
      nZs = size(zs, 1)
      ASSERT(nXs <= nYs)

      if(nZs <= 0) return
      zs = 0
      if(nXs <= 0) return ! this condition includes nYs <= 0 since nXs <= nYs
      if(nXs == 1)then
         zs = xs(1)*ys(1:nZs)
         return
      end if

      iMax = min(nXs - 1, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent (i = 1:iMax)
         zs(i) = dot_product(ys(1:i), xs(i:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nYs, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent(i = nXs:iMax)
         zs(i) = dot_product(ys(i-nXs+1:i), xs(nXs:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nXs + nYs, nZs)
      do concurrent(i = nYs + 1:iMax)
         zs(i) = dot_product(ys(i-nXs+1:nYs), xs(nXs:i-nYs+1:-1))
      end do
   end subroutine conv_short_with_longRealDim1KindREAL128

   subroutine convIntegerDim1KindINT8(xs, ys, zs)
      Integer(kind=INT8), dimension(:), intent(in):: xs, ys
      Integer(kind=INT8), dimension(:), intent(out):: zs

      if(size(xs, 1) >= size(ys, 1))then
         call conv_short_with_long(ys, xs, zs)
      else
         call conv_short_with_long(xs, ys, zs)
      end if
   end subroutine convIntegerDim1KindINT8

   subroutine conv_short_with_longIntegerDim1KindINT8(xs, ys, zs)
      Integer(kind=INT8), dimension(:), intent(in):: xs, ys
      Integer(kind=INT8), dimension(:), intent(out):: zs
      Integer:: nXs, nYs, nZs
      Integer:: i
      Integer:: iMax
      Logical:: isReachedEndOfZs

      nXs = size(xs, 1)
      nYs = size(ys, 1)
      nZs = size(zs, 1)
      ASSERT(nXs <= nYs)

      if(nZs <= 0) return
      zs = 0
      if(nXs <= 0) return ! this condition includes nYs <= 0 since nXs <= nYs
      if(nXs == 1)then
         zs = xs(1)*ys(1:nZs)
         return
      end if

      iMax = min(nXs - 1, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent (i = 1:iMax)
         zs(i) = dot_product(ys(1:i), xs(i:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nYs, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent(i = nXs:iMax)
         zs(i) = dot_product(ys(i-nXs+1:i), xs(nXs:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nXs + nYs, nZs)
      do concurrent(i = nYs + 1:iMax)
         zs(i) = dot_product(ys(i-nXs+1:nYs), xs(nXs:i-nYs+1:-1))
      end do
   end subroutine conv_short_with_longIntegerDim1KindINT8

   subroutine convIntegerDim1KindINT16(xs, ys, zs)
      Integer(kind=INT16), dimension(:), intent(in):: xs, ys
      Integer(kind=INT16), dimension(:), intent(out):: zs

      if(size(xs, 1) >= size(ys, 1))then
         call conv_short_with_long(ys, xs, zs)
      else
         call conv_short_with_long(xs, ys, zs)
      end if
   end subroutine convIntegerDim1KindINT16

   subroutine conv_short_with_longIntegerDim1KindINT16(xs, ys, zs)
      Integer(kind=INT16), dimension(:), intent(in):: xs, ys
      Integer(kind=INT16), dimension(:), intent(out):: zs
      Integer:: nXs, nYs, nZs
      Integer:: i
      Integer:: iMax
      Logical:: isReachedEndOfZs

      nXs = size(xs, 1)
      nYs = size(ys, 1)
      nZs = size(zs, 1)
      ASSERT(nXs <= nYs)

      if(nZs <= 0) return
      zs = 0
      if(nXs <= 0) return ! this condition includes nYs <= 0 since nXs <= nYs
      if(nXs == 1)then
         zs = xs(1)*ys(1:nZs)
         return
      end if

      iMax = min(nXs - 1, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent (i = 1:iMax)
         zs(i) = dot_product(ys(1:i), xs(i:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nYs, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent(i = nXs:iMax)
         zs(i) = dot_product(ys(i-nXs+1:i), xs(nXs:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nXs + nYs, nZs)
      do concurrent(i = nYs + 1:iMax)
         zs(i) = dot_product(ys(i-nXs+1:nYs), xs(nXs:i-nYs+1:-1))
      end do
   end subroutine conv_short_with_longIntegerDim1KindINT16

   subroutine convIntegerDim1KindINT32(xs, ys, zs)
      Integer(kind=INT32), dimension(:), intent(in):: xs, ys
      Integer(kind=INT32), dimension(:), intent(out):: zs

      if(size(xs, 1) >= size(ys, 1))then
         call conv_short_with_long(ys, xs, zs)
      else
         call conv_short_with_long(xs, ys, zs)
      end if
   end subroutine convIntegerDim1KindINT32

   subroutine conv_short_with_longIntegerDim1KindINT32(xs, ys, zs)
      Integer(kind=INT32), dimension(:), intent(in):: xs, ys
      Integer(kind=INT32), dimension(:), intent(out):: zs
      Integer:: nXs, nYs, nZs
      Integer:: i
      Integer:: iMax
      Logical:: isReachedEndOfZs

      nXs = size(xs, 1)
      nYs = size(ys, 1)
      nZs = size(zs, 1)
      ASSERT(nXs <= nYs)

      if(nZs <= 0) return
      zs = 0
      if(nXs <= 0) return ! this condition includes nYs <= 0 since nXs <= nYs
      if(nXs == 1)then
         zs = xs(1)*ys(1:nZs)
         return
      end if

      iMax = min(nXs - 1, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent (i = 1:iMax)
         zs(i) = dot_product(ys(1:i), xs(i:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nYs, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent(i = nXs:iMax)
         zs(i) = dot_product(ys(i-nXs+1:i), xs(nXs:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nXs + nYs, nZs)
      do concurrent(i = nYs + 1:iMax)
         zs(i) = dot_product(ys(i-nXs+1:nYs), xs(nXs:i-nYs+1:-1))
      end do
   end subroutine conv_short_with_longIntegerDim1KindINT32

   subroutine convIntegerDim1KindINT64(xs, ys, zs)
      Integer(kind=INT64), dimension(:), intent(in):: xs, ys
      Integer(kind=INT64), dimension(:), intent(out):: zs

      if(size(xs, 1) >= size(ys, 1))then
         call conv_short_with_long(ys, xs, zs)
      else
         call conv_short_with_long(xs, ys, zs)
      end if
   end subroutine convIntegerDim1KindINT64

   subroutine conv_short_with_longIntegerDim1KindINT64(xs, ys, zs)
      Integer(kind=INT64), dimension(:), intent(in):: xs, ys
      Integer(kind=INT64), dimension(:), intent(out):: zs
      Integer:: nXs, nYs, nZs
      Integer:: i
      Integer:: iMax
      Logical:: isReachedEndOfZs

      nXs = size(xs, 1)
      nYs = size(ys, 1)
      nZs = size(zs, 1)
      ASSERT(nXs <= nYs)

      if(nZs <= 0) return
      zs = 0
      if(nXs <= 0) return ! this condition includes nYs <= 0 since nXs <= nYs
      if(nXs == 1)then
         zs = xs(1)*ys(1:nZs)
         return
      end if

      iMax = min(nXs - 1, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent (i = 1:iMax)
         zs(i) = dot_product(ys(1:i), xs(i:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nYs, nZs)
      isReachedEndOfZs = (iMax == nZs)
      do concurrent(i = nXs:iMax)
         zs(i) = dot_product(ys(i-nXs+1:i), xs(nXs:1:-1))
      end do
      if(isReachedEndOfZs) return

      iMax = min(nXs + nYs, nZs)
      do concurrent(i = nYs + 1:iMax)
         zs(i) = dot_product(ys(i-nXs+1:nYs), xs(nXs:i-nYs+1:-1))
      end do
   end subroutine conv_short_with_longIntegerDim1KindINT64
end module conv_lib