获得Fortran数组的步幅

Question

我目前正在对一些Fortran程序进行现代化改造，我想用一个带有签名的旧式Fortran 77例程编写一个包装器

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

对于这个例程，我写了一个以Fortran 90样式矩阵作为输入的包装器

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。有没有办法访问数组的步幅/扩展以恢复领先的维度？关于与C互操作的ISO_Fortran_binding.h头文件，信息存储在数组描述符中。

为了使问题可视化，这里有一个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是围绕LAPACK的DLASET的以下包装器：

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

这在A的第一列中提供了20个而不是10个和10个零。

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示例中以使其清楚）。

您在新示例中执行的操作是获取每列第一个元素的地址差异，实际上有时会运行。你可以做到，它应该工作。如果您使用公共扩展LOC（以及可选的SIZEOF），或者使用transfer()获取整数值而不是I / O例程，则更容易。注意，一个8字节的整数就足够了，最好是使用INTEGER(C_INTPTR_T)（或ptrdiff）。

---

在我修复了有问题的POINTER并删除了不必要的位之后，请考虑您的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
    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传递给外部过程时制作了副本。