Хранение переменной с многомерным индексом в Fortran

Вопрос

Рассмотрим следующий код:

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_indices порядка 100 (хотя 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], вызовите

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

Однако, помимо того, что это очень расточительное использование памяти, это решение не позволит установить для n_indices число больше 7, учитывая ограничение в 7 измерений для массива.

OBS: этот вопрос является побочным продуктом этого. Я попытался задать вопрос более точно, потерпев неудачу с первой попытки, и мне очень помог ответ @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

storage fortran fortran90
person mzp    schedule 25.05.2018    source источник
comment
Не важно: но теперь ограничение на массивы Fortran составляет 15 рангов.   -  person francescalus    schedule 25.05.2018
comment
Какую реальную проблему вы пытаетесь смоделировать с помощью этого?   -  person Rodrigo Rodrigues    schedule 26.05.2018
comment
@RodrigoRodrigues Еще раз спасибо за вашу помощь. Я отредактировал вопрос, чтобы включить фактический код.   -  person mzp    schedule 26.05.2018

Ответы (1)


arrow_upward
2
arrow_downward

Для любого 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 хотите/нужно.

Изменить (после фактического кода из OP):

Если я правильно расшифровал ваш код (и учитывая то, что вы сказали в своем предыдущем вопросе), вы заявляете:

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

(где L_max = 9, поэтому size(coeff_grid) = 21233664 =~160 МБ), а затем заполнить его:

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

(где count имеет порядок 1000, т. е. 0,005% его элементов), поэтому таким образом вы можете получить значения по его 4 индексам с помощью записи массива.

Пожалуйста, не делай этого. В этом случае вам также не нужна разреженная матрица. Предложенный вами новый подход намного лучше: сохранение индексов в каждой строке меньшего массива и выборка массива коэффициентов по соответствующему расположению этих индексов в его собственном массиве. Это намного быстрее (избегая большого распределения) и намного эффективнее с точки зрения использования памяти.

PS: Обязательно ли придерживаться Fortran 90? Это очень старая версия стандарта, и есть вероятность, что используемый вами компилятор реализует более новую версию. Вы можете значительно улучшить качество своего кода с помощью встроенного move_alloc (для меньшего количества копий массива), типовых констант из встроенного модуля iso_fortran_env (для переносимости), нотации [], >, <, <=,... (для удобочитаемости)...

person Rodrigo Rodrigues    schedule 26.05.2018