将未知大小的数组(子例程输出)传递给另一个子例程

Question

我是英特尔 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
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是不可能的。

Answer 1

您可以从

手册页或 MKL 安装目录中的包含文件mkl_spblas.fi（例如/path/to/mkl/include/）中找到mkl_ddnscsr的界面。

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或此（或此）页面中找到PARDISO的界面。