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

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_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

Answer 1

Для любого 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 (для переносимости), обозначения [], >, <, <=,... (для удобочитаемости) ...