Fortran 中带掩码的移动平均线

Question

我必须在 Fortran 中计算维度为 (7320,8520) 的屏蔽数据集的移动平均值。我写了一个接收数据（TS）并输出平均数据（TS_NEW）的子程序。问题是代码运行花费的时间太长（它实际上从未完成，尽管没有运行进入内存问题）。我想知道是否有一种方法可以使代码更高效。下面是我写的代码：

SUBROUTINE avgwin(ts,winsize,size1,size2,sizelat,sizelon,ts_new)
implicit none
double precision, dimension(size1,size2),INTENT(IN) :: ts
double precision, dimension(winsize,winsize) :: store
double precision, dimension(sizelon,sizelat),INTENT(OUT) :: ts_new
integer :: j,k
integer :: A, B
integer,INTENT(IN) :: winsize,size1,size2,sizelat,sizelon
logical, dimension(size1,size2) :: mask,mask2
double precision :: SUMVAR, COUNTVAR
A=1
B=1

mask = ts > 0 !Mask to highlight all the OK values
mask2 = ts < 0 !Mask to highlight all the values to be discarded
do j=1,sizelat !Looping through latitude
    do k=1,sizelon !Looping through longitude
        if (ALL(mask2(k:k+winsize-1,j:j+winsize-1)) .eqv. .true.) then
            ts_new(B,A) = -100 !Adds a fill value if all the elements are to be discarded
            B=B+1
        else
            SUMVAR = sum(ts(k:k+winsize-1,j:j+winsize-1),     MASK=mask(k:k+winsize-1,j:j+winsize-1))
            COUNTVAR = count(mask(k:k+winsize-1,j:j+winsize-1))
            ts_new(B,A) = SUMVAR/COUNTVAR
            B=B+1
        end if
    end do
    B=1
    A=A+1
end do
END SUBROUTINE

program test
implicit none
double precision, dimension(7320,8520) :: DATA
double precision, dimension(:,:),allocatable :: DATA_NEW
integer :: sizelat, sizelon, i, j, len1, len2, winsize
integer, dimension(3) :: sizes

len1 = 7320
len2 = 8520
do i=1,8520
    do j=1,7320
        DATA(j,i)= i !Just for testing purposes
    end do
end do

sizes(1:3) = (/300,301,302/)
    do w=1,3
       winsize = sizes(w)
       sizelon = len1-winsize+1
       sizelat = len2-winsize+1
       allocate(DATA_NEW(sizelon,sizelat))
       CALL avgwin(DATA,winsize,len1,len2,sizelat,sizelon,DATANEW)
    end do

end program test

Answer 1

虽然不确定这是否符合OP的目的，但如何先沿一个维度收集数据，然后再沿另一维度收集处理后的数据（即部分求和）？例如，如果我们考虑一个更简单的求和 data( 1:L, 1:L ) 而不是移动大小 w 的 window 的问题，可能有三种不同的方法可以实现这一点：

program main
    implicit none
    real, allocatable, dimension(:,:) :: data, direct, part1, part2
    integer :: i1, i2, L, S, w
    real :: t1, t2

    L = 2000
    w = 50
    S = L - w + 1

    allocate( data( L, L ), direct( S, S ), &
              part1( L, S ), part2( S, S ) )

!> test data

    do i2 = 1, L
    do i1 = 1, L
        data( i1, i2 ) = mod( i1 + i2, 2 )
    enddo
    enddo

!> method 1: direct sum (cost = O( S^2 * w^2 ))

    call cpu_time( t1 )
    do i2 = 1, S
    do i1 = 1, S
        direct( i1, i2 ) = sum( data( i1:(i1 + w - 1), i2:(i2 + w - 1) ) )
    enddo
    enddo
    call cpu_time( t2 )
    print *, "time (s) = ", t2 - t1

!> method 2: partial sum (cost = O( S^2 * w * 2 ))

    call cpu_time( t1 )
    do i2 = 1, S
    do i1 = 1, L
        part1( i1, i2 ) = sum( data( i1, i2:(i2 + w - 1) ) )
    enddo
    enddo

    do i2 = 1, S
    do i1 = 1, S
        part2( i1, i2 ) = sum( part1( i1:(i1 + w - 1), i2 ) )
    enddo
    enddo
    call cpu_time( t2 )
    print *, "time (s) = ", t2 - t1

    print *, "error = ", maxval( abs( part2 - direct ) )

!> method 3: an improved version of method 2 (cost = O( S^2 ))

    call cpu_time( t1 )
    do i1 = 1, L
        part1( i1, 1 ) = sum( data( i1, 1:w ) )
        do i2 = 2, S
            part1( i1, i2 ) = part1( i1, i2-1 ) &
                    - data( i1, i2-1 ) + data( i1, i2+w-1 )
        enddo
    enddo

    do i2 = 1, S
        part2( 1, i2 ) = sum( part1( 1:w, i2 ) )
        do i1 = 2, S
            part2( i1, i2 ) = part2( i1-1, i2 ) &
                    - part1( i1-1, i2 ) + part1( i1+w-1, i2 )
        enddo
    enddo
    call cpu_time( t2 )
    print *, "time (s) = ", t2 - t1

    print *, "error = ", maxval( abs( part2 - direct ) )

end program

然后，gfortran-7.2 -O3 test.f90 似乎提供了一些不错的加速：

 time (s) =    9.64789867
 time (s) =   0.345023155
 error =    0.00000000
 time (s) =    8.60958099E-02
 error =    0.00000000

要使用掩码计算移动平均线，类似的方法可能会以某种方式起作用。如果我们搜索网络，可能还有其他（更好的）approaches/libraries这种移动平均线，因为它是很常见的计算...

Fortran 中带掩码的移动平均线

Moving average with mask in Fortran

fortran

moving-average