Получение шага по массиву Фортрана

Question

В настоящее время я модернизирую некоторые программы на Фортране, и я хочу написать обертку вокруг подпрограммы старого стиля Фортран 77 с подписью

SUBROUTINE INITMATRIX( M , N , A , LDA) 
 INTEGER M, N, LDA 
 DOUBLE PRECISION A(LDA,*) 

Для этой подпрограммы я написал обертку, принимающую матрицу стиля Фортран 90 какinput

SUBROUTINE INITMATRIX_F90( A ) 
  DOUBLE PRECISION, INTENT(INOUT) :: A(:,:) 
  INTEGER :: M, N, LDA 
  M = SIZE(A,1) 
  N = SIZE(A,2) 
  LDA = SIZE(A,1) 
  CALL INITMATRIX(M, N, A, LDA)  
END SUBROUTINE

Это работает хорошо, если я не передам часть в процедуру.Например, у меня есть матрица 20 x 20, и я хочу инициализировать только первые 10 строк.Тогда я бы назвал

DOUBLE PRECISION A(20,20) 
CALL INITMATRIX_F90(A(1:10,1:20))

, что приведет к ошибке, так как моя оболочка получает неправильный начальный размер массива.В примере у меня есть LDA=10 вместо LDA = 20.Есть ли способ получить доступ к шагу / расширению массива, чтобы восстановить ведущее измерение?Что касается заголовочного файла ISO_Fortran_binding.h для совместимости с C, информация хранится внутри дескриптора массива.

Чтобы визуализировать проблему, вот MWE, демонстрирующий проблему.

PROGRAM MAIN
    INTERFACE
        SUBROUTINE INITMATRIX_F90(A)
            DOUBLE PRECISION, INTENT(INOUT) :: A(:,:)
        END SUBROUTINE

    END INTERFACE

    DOUBLE PRECISION :: A(20,20)
    A = 1.0D0

    CALL INITMATRIX_F90(A(1:10,1:20))

END PROGRAM MAIN

SUBROUTINE INITMATRIX_F90( A )
    USE ISO_C_BINDING
    IMPLICIT NONE


    DOUBLE PRECISION, INTENT(INOUT), POINTER :: A(:,:)
    INTEGER :: M, N, LDA
    TYPE(C_PTR) :: LOC1, LOC2
    INTEGER*16 :: LOCX1, LOCX2
    CHARACTER*32 :: TMP
    M = SIZE(A,1)
    N = SIZE(A,2)

    WRITE(TMP, *) C_LOC(A(1,1))
    READ(TMP, *) LOCX1
    WRITE(TMP, *) C_LOC(A(1,2))
    READ(TMP, *) LOCX2

    LDA = (LOCX2-LOCX1) / C_SIZEOF(A(1,1))

    WRITE(*,*) "M = ", M
    WRITE(*,*) "N = ", N
    WRITE(*,*) "LOC = ", LOCX1, LOCX2
    WRITE(*,*) "LDA(COMPUTED) = ", LDA

END SUBROUTINE

(я знаю, что указатель отсутствует в интерфейсе, он только для того, чтобы заставить работать C_LOC.)

Выходные данные тогда

 M =           10
 N =           20
 LOC =  140721770410864 140721770411024
 LDA(COMPUTED) =           20

где очевидно, что ведущее измерение правильно вычисляется с помощью грязного хака.Внутренние структуры, используемые компилятором GNU Fortran или привязкой ISO C <-> Fortran (отличной от используемой в GNU), содержат информацию о том, как получить к ним доступ из Fortran без подвохов.

Другой MWE - это следующая обертка вокруг DLASET LAPACK:

PROGRAM MAIN
    INTERFACE
        SUBROUTINE INITMATRIX2_F90(A)
            DOUBLE PRECISION, INTENT(INOUT) :: A(:,:)
        END SUBROUTINE

    END INTERFACE

    DOUBLE PRECISION :: A(20,20)
    A = 0.0D0

    CALL INITMATRIX2_F90(A(1:10,1:20))

    WRITE(*,*) A(1:20,1)

END PROGRAM MAIN

SUBROUTINE INITMATRIX2_F90( A )
    USE ISO_C_BINDING
    IMPLICIT NONE
    DOUBLE PRECISION, INTENT(INOUT) :: A(:,:)
    INTEGER :: M, N, LDA
    EXTERNAL DLASET
    M = SIZE(A,1)
    N = SIZE(A,2)
    CALL DLASET( "All", M, N, 1.0D0, 1.0D0, A(1,1) , M)
END SUBROUTINE

Это дает 20 единиц в первом столбце A вместо десяти единиц и десяти нулей.

Answer 1

Боюсь, что ваши усилия по модернизации приведут к полной запутанности.Если вам нужно позвонить своему INITMATRIX с чем-то, что является ломтиком, не используйте свой F90wrapper, вы ничего не получите.Конечно, не по тому, что вы планируете, это было бы незаконно.

Что произойдет в

DOUBLE PRECISION, INTENT(INOUT) :: A(:,:) 

CALL INITMATRIX(M, N, A, LDA)  

, когда A не является смежным, так это то, что компилятор сделает копию A и передасткопия.Поэтому попытка использовать дескриптор исходного массива будет бесполезной.Даже если он сработает, код, который вы получите в итоге, будет хуже, чем оригинал.

Я рекомендую либо модернизировать сам INITMATRIX, либо просто вызывать его так, как вы его называете.

---

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

CALL INITMATRIX(M, N, A(1,1), LDA) 

Если вы на самом деле делаете это в INITMATRIX_F90, вы должны поместить его в первый INITMATRIX_F90 пример, чтобы прояснить его).

Что вы делаетев вашем новом примере с получением разницы адресов первых элементов каждого столбца, который фактически иногда запускается.Вы можете сделать это, это должно работать.Проще, если вы либо 1. используете общие расширения LOC (и, необязательно, SIZEOF), либо 2. используете transfer(), чтобы получить целочисленное значение, вместо подпрограмм ввода-вывода.Обратите внимание, что 8-байтового целого числа достаточно, и лучше всего использовать INTEGER(C_INTPTR_T) (или ptrdiff).

---

Рассмотрим ваш MWE после того, как я исправил ошибочный POINTER и удалил ненужные биты:

PROGRAM MAIN
    INTERFACE
        SUBROUTINE INITMATRIX_F90(A)
            DOUBLE PRECISION, INTENT(INOUT) :: A(:,:)
        END SUBROUTINE
    END INTERFACE

    DOUBLE PRECISION :: A(20,20)
    A = 1.0D0
    print *,"loc in main",loc(A(1,1))
    CALL INITMATRIX_F90(A(1:10,1:10))
END PROGRAM MAIN

SUBROUTINE INITMATRIX_F90( A )
    IMPLICIT NONE
    DOUBLE PRECISION, INTENT(INOUT) :: A(:,:)

    print *,"LOC in INITMATRIX_F90:",loc(A(1,1))
    call external(A)
END SUBROUTINE

subroutine external(A)
  double precision :: A(*)
  print *,"LOC in external:", loc(A(1))
end subroutine

вывод:

> ./a.out 
 loc in main      140721998532864
 LOC in INITMATRIX_F90:      140721998532864
 LOC in external:             37291664

Как видите, компилятор сделал копию при передаче A во внешнюю процедуру.