将未知大小的数组（子程序输出）传递给另一个子程序

Posted 2023-03-06

【中文标题】将未知大小的数组（子程序输出）传递给另一个子程序

【英文标题】：Passing array of unknown size (subroutine output) to another subroutine

【发布时间】：2016-07-15 23:28:59

【问题描述】：

我是英特尔 MKL 的新手。这是我遇到的一个问题——显然这个问题与 MKL 本身无关，而是与如何声明一个迄今未知大小的数组并将其作为子程序的输出传递给另一个子程序的问题。

我正在尝试使用 mkl_ddnscsr 将矩阵转换为适合 Pardiso 调用的 CSR 格式：

CALL mkl_ddnscsr(job,Nt,Nt,Adns,Nt,Acsr,ja,ia,info)

CALL PARDISO(pt,1,1,11,13,Nt,Acsr,ia,ja,perm,1,iparm,0,b,x,errr)

问题是，在调用 mkl_ddnscsr 子例程之前，我不知道 CSR 形式 Acsr 和索引向量 ja 的长度是多少。 Acsr 和 ja 在主程序中，或者这两行所在的子程序中应该如何声明？

我尝试了类似的东西

INTERFACE
SUBROUTINE mkl_ddnscsr(job, m, n, Adns, lda, Acsr, ja, ia, info)
IMPLICIT NONE
INTEGER :: job(8)
INTEGER :: m, n, lda, info
INTEGER, ALLOCATABLE :: ja(:)
INTEGER :: ia(m+1)
REAL(KIND=8), ALLOCATABLE :: Acsr(:)
REAL(KIND=8) :: Adns(:)
END SUBROUTINE
END INTERFACE

紧随其后

INTEGER, ALLOCATABLE :: ja(:)
REAL(KIND=8), ALLOCATABLE :: Acsr(:)

在INTERFACE之外，在主程序中。但是这个配置在运行时给了我分段错误。

另一方面，如果我尝试类似的东西

INTERFACE
SUBROUTINE mkl_ddnscsr(job, m, n, Adns, lda, Acsr, ja, ia, info)
IMPLICIT NONE
INTEGER :: job(8)
INTEGER :: m, n, lda, info
INTEGER :: ja(:), ia(m+1)
REAL(KIND=8) :: Acsr(:), Adns(:)
END SUBROUTINE
END INTERFACE

然后

INTEGER, DIMENSION(:) :: ja
REAL(KIND=8), DIMENSION(:) :: Acsr

然后ifort会给我以下信息：

error #6596: If a deferred-shape array is intended, then the ALLOCATABLE or POINTER attribute is missing; if an assumed-shape array is intended, the array must be a dummy argument.

有人知道如何解决这个问题吗？在主程序（或主子程序）中声明 ja 和 Acsr 并传递它们的正确方法是什么？

请注意，子例程是英特尔 MKL 包的一部分，不是我自己编写的，所以module 似乎是不可能的。

【问题讨论】：

使用标签 fortran。您不必在标题中重复。 kind=8 是难看的代码味道

【参考方案1】：

您可以从 manual page 或 MKL 安装目录中的包含文件 mkl_spblas.fi 中找到 mkl_ddnscsr 的接口（例如，/path/to/mkl/include/）。

INTERFACE
    subroutine mkl_ddnscsr ( job, m, n, Adns, lda, Acsr, AJ, AI, info )
    integer            job(8)
    integer            m, n, lda, info
    integer            AJ(*), AI(m+1)
    double precision   Adns(*), Acsr(*)
    end
END INTERFACE

因为这个例程只有 Fortran77 样式的虚拟参数（即显式形状数组 AI(m+1) 或假定大小的数组，如 Adns(*)），所以您可以传递任何本地或可分配数组（在调用方分配后）作为实际的论据。此外，显式编写接口块不是强制性的，但它应该对include（在调用方）检测潜在的接口不匹配很有用。

根据手册，mkl_ddnscsr（将密集矩阵转换为稀疏矩阵的例程）看起来像这样：

program main
    implicit none
    ! include 'mkl_spblas.fi'   !! or mkl.fi (not mandatory but recommended)
    integer :: nzmax, nnz, job( 8 ), m, n, lda, info, irow, k
    double precision :: A( 10, 20 )
    double precision, allocatable :: Asparse(:)
    integer, allocatable :: ia(:), ja(:)

    A( :, : ) =  0.0d0
    A( 2, 3 ) = 23.0d0
    A( 2, 7 ) = 27.0d0
    A( 5, 4 ) = 54.0d0
    A( 9, 9 ) = 99.0d0

    !! Give an estimate of the number of non-zeros.
    nzmax = 10

    !! Or assume that non-zeros occupy at most 2% of A(:,:), for example.
    ! nzmax = size( A ) / 50

    !! Or count the number of non-zeros directly.
    ! nzmax = count( abs( A ) > 0.0d0 )

    print *, "nzmax = ", nzmax

    m   = size( A, 1 )   !! number of rows
    n   = size( A, 2 )   !! number of columns
    lda = m              !! leading dimension of A

    allocate( Asparse( nzmax ) )
    allocate( ja( nzmax ) )   !! <-> columns(:)
    allocate( ia( m + 1 ) )   !! <-> rowIndex(:)

    job( 1 )   = 0       !! convert dense to sparse A
    job( 2:3 ) = 1       !! use 1-based indices
    job( 4 )   = 2       !! use the whole A as input
    job( 5 )   = nzmax   !! maximum allowed number of non-zeros
    job( 6 )   = 1       !! generate Asparse, ia, and ja as output

    call mkl_ddnscsr( job, m, n, A, lda, Asparse, ja, ia, info )

    if ( info /= 0 ) then
        print *, "insufficient nzmax (stopped at ", info, "row)"; stop
    endif

    nnz = ia( m+1 ) - 1
    print *, "number of non-zero elements = ", nnz

    do irow = 1, m
        !! This loop runs only for rows having nonzero elements.
        do k = ia( irow ), ia( irow + 1 ) - 1
            print "(2i5, f15.8)", irow, ja( k ), Asparse( k )
        enddo
    enddo

end program

使用ifort -mkl test.f90（使用 ifort14.0）编译会得到预期的结果

 nzmax =           10
 number of non-zero elements =            4
    2    3    23.00000000
    2    7    27.00000000
    5    4    54.00000000
    9    9    99.00000000

至于nzmax的判断，我认为至少有三种方式：（1）只用一个猜测值（如上）； (2) 假设整个数组中非零元素的比例；或 (3) 直接计算密集数组中非零的数量。在任何情况下，因为我们有准确数量的非零作为输出 (nnz)，我们可以重新分配 Asparse 和 ja 以获得准确的大小（如果需要）。

同样，您可以从包含文件mkl_pardiso.fi 或this（或this）页面中找到PARDISO 的接口。

【讨论】：

此页面也可能有助于了解有关 rowIndex(:) hpc.ut.ee/dokumendid/ips_xe_2015/composerxe/Documentation/en_US/…987654324@

的更多信息

这似乎现在可以工作——没有编译错误并且代码运行。但是在计算过程中，我遇到了以下错误消息：内存不足，无法分配 Fortran RTL 消息缓冲区，消息 #174 = hex 000000ae。我的代码反复使用 PARDISO 求解器来求解一组微分方程。关于出了什么问题，我已经没有想法了......

嗯...我不能确定原因，但是这个页面可能有帮助吗？ ***.com/questions/9751996/… 提出一个单独的问题（包括英特尔论坛）可能会很有用。

谷歌搜索该错误消息似乎表明可能的原因是一些内存泄漏（例如，通过使用数组指针而不是可分配数组），尽管完全不确定。在限制迭代次数等后，valgrind 可能会有所帮助。