Fortran中等效于unique的函数

Question

我发现很多问题围绕这个问题转，但没有一个直接回答这个问题的:

-在fortran中，消除整数列表中重复的方式有哪些(a)最快(挂钟)和(b)最优雅(简明清晰)

肯定有比我微弱尝试更好的方法:

Program unique
 implicit none
 !   find "indices", the list of unique numbers in "list"
 integer( kind = 4 ) :: kx, list(10)
 integer( kind = 4 ),allocatable :: indices(:)
 logical :: mask(10)
 !!$    list=(/3,2,5,7,3,1,4,7,3,3/)
 list=(/1,(kx,kx=1,9)/)
 mask(1)=.true.
 do kx=10,2,-1
   mask(kx)= .not.(any(list(:kx-1)==list(kx)))
 end do
 indices=pack([(kx,kx=1,10)],mask)
 print *,indices
End Program unique

我原本期望列表是有序的，但如果放宽这个要求会更好。

- Clinton Winant

---

你看过Rosetta Code吗？ - Arya McCarthy

请参考以下两个网址：https://stackoverflow.com/questions/14137610/remove-repeated-elements-on-an-2d-array-in-fortran-90 和 https://stackoverflow.com/questions/12147864/finding-duplicate-records-in-fortran。 - Vladimir F Героям слава

@AryaMcCarthy Rosetta Code可以找到重复项，但我不能称其为快速或优雅的解决方案：它涉及嵌套的do循环，并且对于相同大小的数组需要比我的示例多一倍的时间... - Clinton Winant

请注意，您自己的示例进行了比较，这是不必要的。但我认为您无法更简洁地编写它（当然，我看不出您可能考虑的向量化方法）。如果您有疑问，没有内在的if语句。您可以使用任何方法编写自己的函数，并且对该函数的调用将很短。 - Vladimir F Героям слава

@HighPerformanceMark 你关于outarr=...的评论让我想起了那个声称有一行CFD代码的人：调用CFD。说真的，虽然对你来说这个剔除任务（去重？）可能不具有挑战性，但它是数百万Fortran代码中的一部分。 - Clinton Winant

Answer 1

我忍不住要写出一个你可能会喜欢的答案。下面的代码将为未排序的整数输入数组按升序返回唯一值数组。请注意，输出结果是实际值，而不仅仅是索引。

program unique_sort
    implicit none
    integer :: i = 0, min_val, max_val
    integer, dimension(10) :: val, unique
    integer, dimension(:), allocatable :: final

    val = [ 3,2,5,7,3,1,4,7,3,3 ]
    min_val = minval(val)-1
    max_val = maxval(val)
    do while (min_val<max_val)
        i = i+1
        min_val = minval(val, mask=val>min_val)
        unique(i) = min_val
    enddo
    allocate(final(i), source=unique(1:i))   !<-- Or, just use unique(1:i) 
    print "(10i5:)", final
end program unique_sort
! output:    1    2    3    4    5    7

请参阅此代码片段，对比使用（unique_sort）的速度、您提供的示例（unique_indices）和Rosetta Code的示例（remove_dups）以及其他几种变体。我想测试@High Performance Mark的代码，但尚未进行。

Run program 1,000,000 times, 100 integers 0<=N<=50
- unique_sort     t~2.1 sec  input: unsorted, w/duplicates  output: sorted unique values
- remove_dup      t~1.4      input: unsorted, w/duplicates  output: unsorted unique values
- unique_indices  t~1.0      input: sorted, w/duplicates    output: unsorted indices for unique values
- BONUS!(Python)  t~4.1      input: unsorted, w/duplicates  output: sorted unique values

底线：在我的机器上（i7 8GB笔记本电脑）unique_indices比remove_dups略快。然而，remove_dups不需要输入数组进行预排序，实际上返回的是值而不是索引（请参见gist中修改后返回值的unique_indices版本，这似乎并没有减慢它的速度）。另一方面，unique_sort需要大约两倍的时间，但是设计用于处理未排序的输入，并以排序顺序返回值，在8个LOC（减去var声明）中。所以这似乎是一个公平的权衡。总之，我相信可以使用某种掩码语句来优化unique_sort以获得更高的速度，但那是另一天的事情。

更新 上面显示的时间是从一个测试程序中获得的，其中每个子程序都放置在一个模块中，并通过过程调用执行。但是，当unique_sort直接放置在主程序中时，我发现性能有了惊人的提升，仅需约0.08秒即可完成100万次运行。通过不使用程序，速度提高了约25倍，这对我来说似乎很奇怪-通常，我认为编译器会将过程调用的成本优化掉。例如，我发现使用过程或将其直接放置在主程序中执行对remove_dup或unique_indices的性能没有影响。

- Matt P

Answer 2

在 @VladimirF 指出我进行了过度比较后，我发现可以对原始代码进行向量化处理（删除 do 循环 do kx....）。我将“unique”函数与基于维基百科的归并排序算法相结合。核心内容包含在 SortUnique 模块中。

Module SortUnique
contains
  Recursive Subroutine MergeSort(temp, Begin, Finish, list)
    ! 1st 3 arguments are input, 4th is output sorted list
    implicit none
    integer(kind=4),intent(inout) :: Begin,list(:),temp(:)
    integer(kind=4),intent(in) :: Finish
    integer(kind=4) :: Middle
    if (Finish-Begin<2) then    !if run size =1
       return                   !it is sorted
    else
       ! split longer runs into halves
       Middle = (Finish+Begin)/2
       ! recursively sort both halves from list into temp
       call MergeSort(list, Begin, Middle, temp)
       call MergeSort(list, Middle, Finish, temp)
       ! merge sorted runs from temp into list
       call Merge(temp, Begin, Middle, Finish, list)
     endif
  End Subroutine MergeSort

  Subroutine Merge(list, Begin, Middle, Finish, temp)
    implicit none
    integer(kind=4),intent(inout) :: list(:),temp(:)
    integer(kind=4),intent(in) ::Begin,Middle,Finish
    integer(kind=4)    :: kx,ky,kz
    ky=Begin
    kz=Middle
    !! While there are elements in the left or right runs...
    do kx=Begin,Finish-1
       !! If left run head exists and is <= existing right run head.
       if (ky.lt.Middle.and.(kz.ge.Finish.or.list(ky).le.list(kz))) then
          temp(kx)=list(ky)
          ky=ky+1
       else
          temp(kx)=list(kz)
          kz = kz + 1
       end if
    end do

  End Subroutine Merge

  Function Unique(list)
    !! usage sortedlist=Unique(list)
    implicit none
    integer(kind=4) :: strt,fin,N
    integer(kind=4), intent(inout) :: list(:)
    integer(kind=4), allocatable  :: unique(:),work(:)
    logical,allocatable :: mask(:)
    ! sort
    work=list;strt=1;N=size(list);fin=N+1
    call MergeSort(work,strt,fin,list) 
    ! cull duplicate indices
    allocate(mask(N));
    mask=.false.
    mask(1:N-1)=list(1:N-1)==list(2:N)
    unique=pack(list,.not.mask)

  End Function Unique

End Module SortUnique

Program TestUnique
  use SortUnique
  implicit none
  !   find "indices", the list of unique numbers in "list"
  integer (kind=4),allocatable :: list(:),newlist(:)
  integer (kind=4)  :: kx,N=100000 !N  even
  real (kind=4) :: start,finish,myrandom

  allocate(list(N))
  do kx=1,N
     call random_number(myrandom)
     list(kx)=ifix(float(N)/2.*myrandom)
  end do

  call cpu_time(start)

  newlist=unique(list)

  call cpu_time(finish)
  print *,"cull duplicates: ",finish-start
  print *,"size(newlist) ",size(newlist)

End Program TestUnique

在@HighPerformanceMark的建议下，该函数被简单地调用为newlist=unique(list)。以上内容肯定不够简洁，但似乎很清晰，并且比我的原始解决方案或其他提出的解决方案快大约200倍。

- Clinton Winant

---

当我尝试运行使用ifort 17编译的这个程序时，我遇到了forrtl: severe (157): Program Exception - access violation的错误。 - Matt P

@MattP 我从gfortran那里没有收到任何投诉，编译后的程序也可以正确运行。我想我会遵循弗拉基米尔·F所说的规则，即如果在你的编译器上可以工作，就不要担心它。说真的，我没有任何ifort版本。 - Clinton Winant

可以的。反正我有机会时可能会试着用ifort编译它。为了达到你所观察到的性能，这可能值得一试。 - Matt P

@MattP 非常感谢，我已经编辑了代码以修复你发现的错误。 - Clinton Winant

@ClintonWinant 当我尝试使用gfortran运行您的代码时，我收到以下错误信息：在sortunique_mod.f90文件的第39行 Fortran运行时错误：数组'list'的第1维的索引“100001”超出了100000的上限。 - geom