我有一个二维实数数组，我想找到 n 个最大值并将这些最大值分配给 1，将所有其他值分配给 0。
以下代码通过在循环内使用 MAXLOC 找到最大值，将其更改为 -9999，从而将其排除在循环的下一次迭代之外，从而正确地完成了此操作。最后，所有 -9999 值都分配给 1。问题是这种方法非常慢。 这个样本数据集很好，它有 8604 个单元格，其中选择了最高的 50 个，但真实数据有 108183756 个单元格，我需要找到最高的 1672137 个单元格。我会很感激任何帮助。
PROGRAM mapalgebra
IMPLICIT NONE
REAL, DIMENSION(64,126) :: n
REAL, DIMENSION(64,126) :: a
REAL, DIMENSION(64,126) :: r
REAL, DIMENSION(64,126) :: ine
REAL, DIMENSION(64,126) :: s
REAL, DIMENSION(64,126) :: p
REAL, DIMENSION(64,126) :: TTP
INTEGER, DIMENSION(64,126) :: newlu
INTEGER :: row,col,max_rows,max_cols,i, j,demand, si
INTEGER, ALLOCATABLE :: AR1(:)
newlu = 0
max_rows=64
max_cols=126
OPEN(UNIT=1, FILE="urb.txt")
OPEN(UNIT=2, FILE="acc.txt")
OPEN(UNIT=4, FILE="ran.txt")
OPEN(UNIT=5, FILE="lu1.txt")
OPEN(UNIT=7, FILE="suit.txt")
OPEN(UNIT=8, FILE="pob.txt")
DO row = 1,max_rows
READ(1,*) (n(row,col),col=1,max_cols)
END DO
DO row = 1,max_rows
READ(2,*) (a(row,col),col=1,max_cols)
END DO
DO row = 1,max_rows
READ(4,*) (r(row,col),col=1,max_cols)
END DO
DO row = 1,max_rows
READ(5,*) (ine(row,col),col=1,max_cols)
END DO
DO row = 1,max_rows
READ(7,*) (s(row,col),col=1,max_cols)
END DO
DO row = 1,max_rows
READ(8,*) (p(row,col),col=1,max_cols)
END DO
demand = 50
print *, "Land use demand is : ", demand
TTP = (ine+n+r+a+s)*p !population weighted model, (inertia+nhood+randomness+accessibility+suitability)*population
si = SIZE(SHAPE(TTP)) ! Get number of dimensions
! in array
print *, "TTP has : ", si , " dimensions"
ALLOCATE (AR1(si)) ! Allocate AR1 to number
! of dimensions in array
DO i=1,demand
AR1=MAXLOC (TTP)
!print *, "MAXLOC (TTP) : ", AR1
TTP(AR1(1),AR1(2)) = -9999
newlu(AR1(1),AR1(2)) = 1
END DO
OPEN(UNIT=3, FILE="output.txt", ACTION="write", STATUS="replace")
WRITE(3,11) 'ncols 126'
WRITE(3,11) 'nrows 64'
WRITE(3,11) 'xllcorner 461229.21505206'
WRITE(3,11) 'yllcorner 4478822.89048722'
WRITE(3,11) 'cellsize 99.38142966 '
WRITE(3,11) 'NODATA_value 0 '
11 format(A,I3)
DO i=1,max_rows
WRITE(3,*) (newlu(i,j), j=1,max_cols)
END DO
CLOSE (1)
CLOSE (2)
CLOSE (3)
CLOSE (4)
CLOSE (5)
CLOSE (7)
CLOSE (8)
END PROGRAM mapalgebra
这是一个解决方案。思路是统计
>= v0v0kO(N*log(N))Program test
implicit none
real, allocatable :: a(:)
integer, parameter :: n=108183756, k=1672137
integer :: k0, i, kk
real :: v0, v1, v2, vmin, vmax
allocate(a(n))
call random_number(a)
vmax = maxval(a,dim=1)
vmin = minval(a,dim=1)
write(*,*) vmin,vmax
if (count(a == vmax) > k) then
! special case
kk = k
do i = 1, n
if (a(i) == vmax .and. kk > 0) then
a(i) = 1.0
kk = kk - 1
else
a(i) = 0.0
end if
end do
else
! ...
v1 = vmax
v2 = vmin ! the v2 initial value could be refined
! in the case where k << n
! Reduce the [v1,v2] interval by dichotomy
v0 = (v1+v2)/2.0
k0 = count(a >= v0)
do
write(*,*) "-----",v1,v2,v0,k0
if (k0 == k) exit ! perfect value
if (v0 == v2 .or. v0 == v1) exit ! see below
if (k0 > k) then
v2 = v0
else
v1 = v0
end if
v0 = (v1+v2)/2.0
k0 = count(a >= v0)
end do
! If the previous loop exited because v0 was equal to
! either v1 or v2, then it means that there were no
! more possible real value between v1 and v2.
! In turns it means that there are several elements of
! a(:) that are equal to the minimum value (v2) among
! the k largest. This has to be taken into account
if (k0 == k) then
kk = 0
else
v0 = v1
kk = k-count(a >= v0)
end if
do i = 1, n
if (a(i) >= v0) then
a(i) = 1.0
else if (a(i) >= v2 .and. kk > 0) then
a(i) = 1.0
kk = kk - 1
else
a(i) = 0.0
end if
end do
end if
write(*,*) count(a == 1.0)
End