在 Fortran 中存储具有多维索引的变量

Question

考虑以下代码：

program example

  implicit none

  integer, parameter :: n_coeffs = 1000
  integer, parameter :: n_indices = 5
  integer :: i
  real(8), dimension(n_coeffs) :: coeff
  integer, dimension(n_coeffs,n_indices) :: index

  do i = 1, n_coeffs
    coeff(i)   = real(i*3,8)
    index(i,:) = [2,4,8,16,32]*i
  end do

end

对于任何 5 维索引我需要获得相关系数，在不知道或计算的情况下 i。例如，给定[2,4,8,16,32]我需要获得3.0无需计算i.

是否有一个合理的解决方案，也许使用稀疏矩阵，适用于n_indices100 左右（尽管n_coeffs仍然是 1000 的数量级）？

一个糟糕的解决方案

一种解决方案是定义一个 5 维数组，如下所示

real(8), dimension(2000,4000,8000,16000,32000) :: coeff2

do i = 1, ncoeffs
    coeff2(index(i,1),index(i,2),index(i,3),index(i,4),index(i,5)) = coeff(i)
end do

然后，得到与相关的系数[2,4,8,16,32], call

coeff2(2,4,8,16,32)

然而，除了非常浪费内存之外，这个解决方案不允许n_indices考虑到数组的 7 维限制，要设置为大于 7 的数字。

OBS：这个问题是this one。我试图更准确地提出这个问题，但在第一次尝试中失败了，这一努力极大地受益于@Rodrigo_Rodrigues 的回答。

实际代码

如果它有帮助，这里是我试图解决的实际问题的代码。它是一种用于逼近函数的自适应稀疏网格方法。主要目标是使插值在 和 处尽可能快：

MODULE MOD_PARAMETERS

    IMPLICIT NONE
    SAVE

    INTEGER, PARAMETER :: d     = 2     ! number of dimensions
    INTEGER, PARAMETER :: L_0   = 4     ! after this adaptive grid kicks in, for L <= L_0 usual sparse grid
    INTEGER, PARAMETER :: L_max = 9    ! maximum level
    INTEGER, PARAMETER :: bound = 0     ! 0 -> for f = 0 at boundary
                                        ! 1 -> adding grid points at boundary
                                        ! 2 -> extrapolating close to boundary
    INTEGER, PARAMETER :: max_error = 1
    INTEGER, PARAMETER :: L2_error  = 1
    INTEGER, PARAMETER :: testing_sample = 1000000

    REAL(8), PARAMETER :: eps = 0.01D0   ! epsilon for adaptive grid

END MODULE MOD_PARAMETERS

PROGRAM MAIN

USE MOD_PARAMETERS
IMPLICIT NONE

INTEGER, DIMENSION(d,d)                :: ident
REAL(8), DIMENSION(d)                  :: xd
INTEGER, DIMENSION(2*d)                :: temp
INTEGER, DIMENSION(:,:),   ALLOCATABLE :: grid_index, temp_grid_index, grid_index_new, J_index
REAL(8), DIMENSION(:),     ALLOCATABLE :: coeff, temp_coeff, J_coeff

REAL(8) :: temp_min, temp_max, V, T, B, F, x1
INTEGER :: k, k_1, k_2, h, i, j, L, n, dd, L1, L2, dsize, count, first, repeated, add, ind

INTEGER :: time1, time2, clock_rate, clock_max

REAL(8), DIMENSION(L_max,L_max,2**(L_max),2**(L_max)) :: coeff_grid
INTEGER, DIMENSION(d) :: level, LL, ii

REAL(8), DIMENSION(testing_sample,d) :: x_rand
REAL(8), DIMENSION(testing_sample)   :: interp1, interp2

! ============================================================================
! EXECUTABLE
! ============================================================================

ident = 0
DO i = 1,d
    ident(i,i) = 1
ENDDO

! Initial grid point
dsize = 1
ALLOCATE(grid_index(dsize,2*d),grid_index_new(dsize,2*d))
grid_index(1,:) = 1
grid_index_new  = grid_index

ALLOCATE(coeff(dsize))
xd = (/ 0.5D0, 0.5D0 /)
CALL FF(xd,coeff(1))
CALL FF(xd,coeff_grid(1,1,1,1))

L = 1
n = SIZE(grid_index_new,1)
ALLOCATE(J_index(n*2*d,2*d))
ALLOCATE(J_coeff(n*2*d))

CALL SYSTEM_CLOCK (time1,clock_rate,clock_max)

DO WHILE (L .LT. L_max)

    L       = L+1
    n       = SIZE(grid_index_new,1)
    count   = 0
    first   = 1
    DEALLOCATE(J_index,J_coeff)
    ALLOCATE(J_index(n*2*d,2*d))
    ALLOCATE(J_coeff(n*2*d))
    J_index = 0
    J_coeff = 0.0D0
    DO k = 1,n
        DO i = 1,d
            DO j = 1,2
                IF ((bound .EQ. 0) .OR. (bound .EQ. 2)) THEN
                    temp = grid_index_new(k,:)+(/ident(i,:),ident(i,:)*(grid_index_new(k,d+i)-(-1)**j)/)
                ELSEIF (bound .EQ. 1) THEN
                    IF (grid_index_new(k,i) .EQ. 1) THEN
                        temp = grid_index_new(k,:)+(/ident(i,:),ident(i,:)*(-(-1)**j)/)
                    ELSE
                        temp = grid_index_new(k,:)+(/ident(i,:),ident(i,:)*(grid_index_new(k,d+i)-(-1)**j)/)
                    ENDIF
                ENDIF
                CALL XX(d,temp(1:d),temp(d+1:2*d),xd)
                temp_min = MINVAL(xd)
                temp_max = MAXVAL(xd)
                IF ((temp_min .GE. 0.0D0) .AND. (temp_max .LE. 1.0D0)) THEN
                    IF (first .EQ. 1) THEN
                        first = 0
                        count = count+1
                        J_index(count,:) = temp
                        V = 0.0D0
                        DO k_1 = 1,SIZE(grid_index,1)
                            T = 1.0D0
                            DO k_2 = 1,d
                                CALL XX(1,temp(k_2),temp(d+k_2),x1)
                                CALL BASE(x1,grid_index(k_1,k_2),grid_index(k_1,k_2+d),B)
                                T = T*B
                            ENDDO
                            V = V+coeff(k_1)*T
                        ENDDO
                        CALL FF(xd,F)
                        J_coeff(count) = F-V
                    ELSE
                        repeated = 0
                        DO h = 1,count
                            IF (SUM(ABS(J_index(h,:)-temp)) .EQ. 0) THEN
                                repeated = 1
                            ENDIF
                        ENDDO
                        IF (repeated .EQ. 0) THEN
                            count = count+1
                            J_index(count,:) = temp
                            V = 0.0D0
                            DO k_1 = 1,SIZE(grid_index,1)
                                T = 1.0D0
                                DO k_2 = 1,d
                                    CALL XX(1,temp(k_2),temp(d+k_2),x1)
                                    CALL BASE(x1,grid_index(k_1,k_2),grid_index(k_1,k_2+d),B)
                                    T = T*B
                                ENDDO
                                V = V+coeff(k_1)*T
                            ENDDO
                            CALL FF(xd,F)
                            J_coeff(count) = F-V
                        ENDIF
                    ENDIF
                ENDIF
            ENDDO
        ENDDO
    ENDDO

    ALLOCATE(temp_grid_index(dsize,2*d))
    ALLOCATE(temp_coeff(dsize))
    temp_grid_index = grid_index
    temp_coeff      = coeff
    DEALLOCATE(grid_index,coeff)
    ALLOCATE(grid_index(dsize+count,2*d))
    ALLOCATE(coeff(dsize+count))
    grid_index(1:dsize,:) = temp_grid_index
    coeff(1:dsize)        = temp_coeff
    DEALLOCATE(temp_grid_index,temp_coeff)
    grid_index(dsize+1:dsize+count,:) = J_index(1:count,:)
    coeff(dsize+1:dsize+count)        = J_coeff(1:count)
    dsize = dsize + count    

    DO i = 1,count
        coeff_grid(J_index(i,1),J_index(i,2),J_index(i,3),J_index(i,4)) = J_coeff(i)
    ENDDO

    IF (L .LE. L_0) THEN
        DEALLOCATE(grid_index_new)
        ALLOCATE(grid_index_new(count,2*d))
        grid_index_new = J_index(1:count,:)
    ELSE
        add = 0
        DO h = 1,count
            IF (ABS(J_coeff(h)) .GT. eps) THEN
                add = add + 1
                J_index(add,:) = J_index(h,:)
            ENDIF
        ENDDO
        DEALLOCATE(grid_index_new)
        ALLOCATE(grid_index_new(add,2*d))
        grid_index_new = J_index(1:add,:)
    ENDIF

ENDDO

CALL SYSTEM_CLOCK (time2,clock_rate,clock_max)
PRINT *, 'Elapsed real time1 = ', DBLE(time2-time1)/DBLE(clock_rate)
PRINT *, 'Grid Points        = ', SIZE(grid_index,1)

! ============================================================================
! Compute interpolated values:
! ============================================================================

CALL RANDOM_NUMBER(x_rand)
CALL SYSTEM_CLOCK (time1,clock_rate,clock_max)
DO i = 1,testing_sample
    V = 0.0D0
    DO L1=1,L_max
        DO L2=1,L_max
            IF (L1+L2 .LE. L_max+1) THEN
                level = (/L1,L2/)
                T = 1.0D0
                DO dd = 1,d
                    T = T*(1.0D0-ABS(x_rand(i,dd)/2.0D0**(-DBLE(level(dd)))-DBLE(2*FLOOR(x_rand(i,dd)*2.0D0**DBLE(level(dd)-1))+1)))
                ENDDO
                V = V + coeff_grid(L1,L2,2*FLOOR(x_rand(i,1)*2.0D0**DBLE(L1-1))+1,2*FLOOR(x_rand(i,2)*2.0D0**DBLE(L2-1))+1)*T
            ENDIF
        ENDDO
    ENDDO
    interp2(i) = V
ENDDO
CALL SYSTEM_CLOCK (time2,clock_rate,clock_max)
PRINT *, 'Elapsed real time2 = ', DBLE(time2-time1)/DBLE(clock_rate)

END PROGRAM

对于任何 5 维索引，我需要获取关联的 系数，在不知道或计算的情况下 i。例如，给定[2,4,8,16,32]我需要获得3.0无需计算i.

  function findloc_vector(matrix, vector) result(out)
    integer, intent(in) :: matrix(:, :)
    integer, intent(in) :: vector(size(matrix, dim=2))
    integer :: out, i

    do i = 1, size(matrix, dim=1)
      if (all(matrix(i, :) == vector)) then
        out = i
        return
      end if
    end do
    stop "No match for this vector"
  end

这就是你如何使用它：

print*, coeff(findloc_vector(index, [2,4,8,16,32])) ! outputs 3.0

我必须承认我不愿意发布这段代码，因为即使这样回答你的问题，老实说我认为这是not你真正想要/需要什么，但你没有提供足够的信息让我知道你真正想要什么do想要/需要。

Edit（在OP的实际代码之后）：

如果我正确解密了您的代码（并考虑到您在上一个问题中所说的内容），您将声明：

REAL(8), DIMENSION(L_max,L_max,2**(L_max),2**(L_max)) :: coeff_grid

(where L_max = 9, so size(coeff_grid) = 21233664=~160MB），然后用以下内容填充：

DO i = 1,count
    coeff_grid(J_index(i,1),J_index(i,2),J_index(i,3),J_index(i,4)) = J_coeff(i)
ENDDO

(where count的数量级为 1000，即其元素的 0.005%），因此这样您可以使用数组表示法通过其 4 个索引来获取值。

请不要这样做。在这种情况下，您也不需要稀疏矩阵。您提出的新方法要好得多：将索引存储在较小数组的每一行中，并通过这些索引在其自己的数组中的相应位置获取系数数组。这更快（避免大量分配）并且内存效率更高。

PS：你必须坚持使用 Fortran 90 吗？它是标准的一个非常旧的版本，您使用的编译器很可能实现了更新的版本。您可以通过内在的方法大大提高代码质量move_alloc（对于较少的数组副本），内部模块中的类型常量iso_fortran_env（为了便携性），[], >, <, <=,...符号（为了可读性）...