Passer un tableau de taille inconnue (sortie de sous-programme) vers un autre sous-programme

Question

Je suis nouveau sur Intel MKL. Voici un problème que j'ai rencontré - apparemment un problème non lié à MKL lui-même, mais au problème de comment déclarer et passer un tableau de taille jusqu'ici inconnue comme sortie d'un sous-programme vers un autre sous-programme.

J'essaie d'utiliser mkl_ddnscsr pour convertir une matrice à son format RSE adapté pour les appels par 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) 

Le problème est, je ne sais pas ce que la longueur du formulaire CSR Acsr et l'indice vector ja avant d'appeler le sous-programme mkl_ddnscsr. Comment devrait-on déclarer Acsr et ja dans le programme principal, ou le sous-programme où se trouvent ces deux lignes?

J'ai essayé quelque chose comme

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 

suivie

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

dehors de l'interface, dans le programme principal. Mais cette configuration me donne la faute de segmentation en cours d'exécution.

D'autre part, si je tente quelque chose comme

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 

puis

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

Puis ifort me donner le message suivant:

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. 

Quelqu'un at-il une idée comment travailler autour de cela? Quelle est la bonne façon de déclarer ja et Acsr dans le programme principal (ou le sous-programme principal) et de les transmettre?

Notez que les sous-programmes font partie du package Intel MKL, pas quelque chose que j'écris seul, donc il semble que module serait hors de question.

Answer 1

Vous pouvez trouver l'interface pour mkl_ddnscsr à partir du manual page ou du fichier d'inclusion mkl_spblas.fi dans votre répertoire d'installation MKL (par exemple,/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 

Parce que cette routine n'a que des arguments factices Fortran77 de style (c.-à-tableau-forme explicite AI(m+1) ou un tableau de taille supposée comme Adns(*)), vous pouvez passer des réseaux locaux ou affectables (après alloués dans le côté de l'appelant) en tant qu'arguments réels. En outre, il n'est pas obligatoire d'écrire un bloc d'interface de manière explicite, mais il devrait être utile de le détecter (du côté de l'appelant) pour détecter une non-concordance d'interface potentielle.

Selon le manuel, il ressemble à mkl_ddnscsr (une routine pour la conversion d'une dense matrice creuse) fonctionne quelque chose comme ceci:

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 

avec ifort -mkl test.f90 Compiler (avec ifort14.0) donne le résultat attendu

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

En ce qui concerne la détermination de nzmax, je pense qu'il ya au moins trois façons pour cela: (1) il suffit d'utiliser une valeur de deviner (comme ci-dessus); (2) supposer la fraction d'éléments non-zéro dans le tableau entier; ou (3) compter directement le nombre de nonzeros dans le tableau dense. Dans tous les cas, parce que nous avons le nombre exact de nonzeros en sortie (nnz), nous pourrions réallouer Asparse et ja pour avoir la taille exacte (si nécessaire).

De même, vous pouvez trouver l'interface pour PARDISO à partir du fichier comprennent mkl_pardiso.fi ou de this (ou this) page.