Passing array of unknown size (subroutine output) to another subroutine

Question

I'm new to Intel MKL. Here's a problem I've come across -- apparently a problem not related to MKL itself, but to the problem of how to declare and pass an array of hitherto unknown size as an output of a subroutine to another subroutine.

I'm trying to use mkl_ddnscsr to convert a matrix to its CSR format suitable for calling by Pardiso:

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)

Problem is, I have no idea what the length of the CSR form Acsr and the index vector ja before calling the mkl_ddnscsr subroutine. How should one declare Acsr and ja in the main program, or the subroutine where these two lines are located?

I tried something like

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

followed by

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

outside the INTERFACE, in the main program. But this configuration gives me the segmentation fault upon running.

On the other hand, if I try something like

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

and then

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

Then ifort would give me the following message:

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.

Anyone got any idea how to work around this? What's the right way to declare ja and Acsr in the main program (or main subroutine) and pass them around?

Note that the subroutines are part of the Intel MKL package, not something I write on my own, so it appears that module would be out of the question.

Answer 1

You can find the interface for mkl_ddnscsr from the manual page, or from the include file mkl_spblas.fi in your MKL install directory (e.g., /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

Because this routine has only Fortran77-style dummy arguments (i.e., explicit-shape array AI(m+1) or assumed-size arrays like Adns(*)), you can pass any local or allocatable arrays (after allocated in the caller side) as actual arguments. Also, it is not mandatory to write an interface block explicitly, but it should be useful to include it (in the caller side) to detect potential interface mismatch.

According to the manual, it looks like mkl_ddnscsr (a routine for converting a dense to sparse matrix) works something like this:

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

Compiling with ifort -mkl test.f90 (with ifort14.0) gives the expected result

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

As for the determination of nzmax, I think there are at least three ways for this: (1) just use a guess value (as above); (2) assume the fraction of nonzero elements in the whole array; or (3) directly count the number of nonzeros in the dense array. In any case, because we have the exact number of nonzeros as output (nnz), we could re-allocate Asparse and ja to have the exact size (if necessary).

Similarly, you can find the interface for PARDISO from the include file mkl_pardiso.fi or from this (or this) page.