mumps en subroutine

Bonjour,
Je fais appel au package mumps dans un programme principal. Le code est bon, il affiche les bons résultats.
Le code est:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
 
PROGRAM MUMPS_TEST
implicit none
INCLUDE 'mpif.h'
INCLUDE 'dmumps_struc.h'
TYPE (DMUMPS_STRUC) mumps_par
 
     	INTEGER                IERR, II,k,l,i,j,nblignes,jj,MZ,Nb,NbZ
     	double precision       Q(1:5,1:5) 
   	double precision       BB(1:5), param(1:43)
	double precision       :: valeur
 
  	open (unit = 60,file = 'entrees.txt', status='unknown') ! fichier des input qui doit contenir :N,NZ, les éléments de la matrice A et du vecteur B
 	open (unit = 61,file = 'param.txt', status='unknown') 
 
 
 
 
 
!remplissage du fichier entrees.txt
 
Nb=5
NbZ=12
 
do k=1,Nb 
do l=1,Nb
Q(k,l)=0
enddo
enddo
!2) remplissage des éléments diagonaux
Q(1,2)= 3.0
Q(2,3)=-3.0
Q(4,3)= 2.0
Q(5,5)= 1.0
Q(2,1)= 3.0
Q(1,1)= 2.0
Q(5,2)= 4.0
Q(3,4)= 2.0
Q(2,5)= 6.0
Q(3,2)= -1.0
Q(1,3)= 4.0
Q(3,3)= 1.0 
 
!vecteur B (2 ème membre du système linéaire) tel que Ax=B
BB(1)=20.0
BB(2)=24.0
BB(3)=9.0
BB(4)=6.0
BB(5)=13.0
 
 
!AFFICHAGE 
write(60,'(I1)') Nb !affichage du rang de Q
 
write(60,'(I2)') NbZ  !affichage du nbre des nonzéro de A sur 4 longueurs pour éviter les espaces en début de ligne du fichier input_1_iter (même si cet estpace n'est pas gênant)
 
! affichage de Q
do k=1,5 !je parcous A jusqu'à la ligne MN
	do l=1,5 !je parcours toutes les colonnes
 
		if (Q(k,l)/=0) then
		write(60,100) k,l,Q(k,l)
		100 FORMAT(I1,/,I1,/,f4.1) !format réel réparti sur 12 longueurs avec 10 chiffres après la virgule 
		end if
 
	end do
end do
 
!affichage du vecteur B
do i=1,5
	write(60,105)BB(i)
	105 FORMAT(f4.1)
end do
 
 close(60)
 
 
! lecture du fichier entrees.txt et attribution de ses valeurs au tableau param
open (60,file='entrees.txt',action='read') ! on ouvre le fichier (,status='old')
nblignes=43 !nblignes=2+(NZ*3)+(dimension du vecteur B)=2+(12*3)+5=43
 
   do jj=1,nblignes
   ! on lit une ligne du fichier
   read(60,*) valeur
  ! on stocke la donnee lue dans la case correspondante du tableau   
  param(jj)=valeur
write (61, *) param(jj)
  end do
 close(60)
 
!!!Début mumps    
 
      CALL MPI_INIT(IERR)
!do k=1,5 !boucle dans le temps
 
 
!les composantes: COMM, SYM et PAR doivent être définis par l'utilisateur avant a call with JOB=-1
! Define a communicator for the package.
      mumps_par%COMM = MPI_COMM_WORLD !mis avant l'initialisation (JOB=-1)
!  Initialize an instance of the package
!  for L U factorization (sym = 0, with working host)
 
      mumps_par%JOB = -1 !initialise l'exemple par des valeurs par défaults (comme ICNTL) avant tout call mumps
      mumps_par%SYM = 0! matrice asymétrique
      mumps_par%PAR = 1 !PAR=1 obligatoire dans la version séquentielle
      CALL DMUMPS(mumps_par)
      !open (unit = 17,file = 'input_simpletest_real')
      IF (mumps_par%INFOG(1).LT.0) THEN !informations de tous les processeurs
       WRITE(6,'(A,A,I6,A,I9)') " ERROR RETURN: ","  mumps_par%INFOG(1)= ", mumps_par%INFOG(1), &
            "  mumps_par%INFOG(2)= ", mumps_par%INFOG(2) 
       GOTO 500
      END IF
!  Define problem on the host (processor 0)
	IF ( mumps_par%MYID .eq. 0 ) THEN
!définir une matrice centralized assembled avant la factorisation (Job=2) et l'analyse (job=1)
        mumps_par%N=Nb
        mumps_par%NZ=NbZ
        ALLOCATE( mumps_par%IRN ( NbZ ) )
        ALLOCATE( mumps_par%JCN ( NbZ ) )
        ALLOCATE( mumps_par%A( NbZ ) )
        ALLOCATE( mumps_par%RHS ( Nb ) )
 
		j=1
		DO II = 1, NbZ
		mumps_par%IRN(II)=param(j+2)
		mumps_par%JCN(II)=param(j+3)
		mumps_par%A(II)  =param(j+4)
		j=j+3
		end do
 
		MZ=(NbZ*3)+2
		DO II = 1, Nb
		mumps_par%RHS(II)=param(II+MZ)
		END DO
	END IF
 
 
!  Call package for solution
      mumps_par%JOB = 6!combines the actions of calls with JOB=1 (analyse), JOB=2(factorisation), and JOB=3(solution)
      CALL DMUMPS(mumps_par)
      IF (mumps_par%INFOG(1).LT.0) THEN
       WRITE(6,'(A,A,I6,A,I9)') " ERROR RETURN: ", "  mumps_par%INFOG(1)= ", mumps_par%INFOG(1), &
            "  mumps_par%INFOG(2)= ", mumps_par%INFOG(2) 
       GOTO 500
      END IF
!  Solution has been assembled on the host
      IF ( mumps_par%MYID .eq. 0 ) THEN
        WRITE(6, * ) ' Solution is ',(mumps_par%RHS(II),II=1,Nb) !afficher à l'écran la solution
      END IF
!  Deallocate user data
      IF ( mumps_par%MYID .eq. 0 )THEN
        DEALLOCATE( mumps_par%IRN )! can be deallocated after analysis
        DEALLOCATE( mumps_par%JCN )! can be deallocated after analysis
        DEALLOCATE( mumps_par%A   )
        DEALLOCATE( mumps_par%RHS )
      END IF
!  Destroy the instance (deallocate internal data structures)
      mumps_par%JOB = -2 !termine l'exemple: toutes les data de la structure sont déallocated (sauf pour celles mises par l'utilsateur): on le mets au dernier appel de Dmumps
      CALL DMUMPS(mumps_par)
      IF (mumps_par%INFOG(1).LT.0) THEN
       WRITE(6,'(A,A,I6,A,I9)') " ERROR RETURN: ", "  mumps_par%INFOG(1)= ", mumps_par%INFOG(1), &
            "  mumps_par%INFOG(2)= ", mumps_par%INFOG(2) 
       GOTO 500
      END IF
!end do !fin boucle temps
 500  CALL MPI_FINALIZE(IERR)
      STOP
  	close (61)
END

Par contre, quand j'essaye de faire appel à ce package dans une subroutine, il plante (erreur: ERROR RETURN: mumps_par%INFOG(1)= -2 mumps_par%INFOG(2)= 0
: çàd qu'il n'arrive pas à lire Nb et NbZ.

le code est le suivant:

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
 
PROGRAM MUMPS_TEST
implicit none
     	INTEGER               Nb,NbZ,k,l,i ,nblignes,jj   
     	double precision       Q(1:5,1:5) 
   	double precision       BB(1:5), param(1:43)
	double precision       :: valeur
 
 
  	open (unit = 60,file = 'entrees.txt', status='unknown') ! fichier des input qui doit contenir :N,NZ, les éléments de la matrice A et du vecteur B
 	open (unit = 61,file = 'param.txt', status='unknown') 
!remplissage du fichier entrees.txt
 
Nb=5
NbZ=12
 
do k=1,Nb 
do l=1,Nb
Q(k,l)=0
enddo
enddo
!2) remplissage des éléments diagonaux
Q(1,2)= 3.0
Q(2,3)=-3.0
Q(4,3)= 2.0
Q(5,5)= 1.0
Q(2,1)= 3.0
Q(1,1)= 2.0
Q(5,2)= 4.0
Q(3,4)= 2.0
Q(2,5)= 6.0
Q(3,2)= -1.0
Q(1,3)= 4.0
Q(3,3)= 1.0 
 
!vecteur B (2 ème membre du système linéaire) tel que Ax=B
BB(1)=20.0
BB(2)=24.0
BB(3)=9.0
BB(4)=6.0
BB(5)=13.0
 
 
!AFFICHAGE 
write(60,'(I1)') Nb   !affichage du rang de Q
write(60,'(I2)') NbZ  !affichage du nbre des nonzéro de A sur 4 longueurs pour éviter les espaces en début de ligne du fichier input_1_iter (même si cet estpace n'est pas gênant)
 
! affichage de Q
do k=1,5 !je parcous A jusqu'à la ligne MN
	do l=1,5 !je parcours toutes les colonnes
 
		if (Q(k,l)/=0) then
		write(60,100) k,l,Q(k,l)
		100 FORMAT(I1,/,I1,/,f4.1) !format réel réparti sur 12 longueurs avec 10 chiffres après la virgule 
		end if
 
	end do
end do
 
!affichage du vecteur B
do i=1,5
	write(60,105)BB(i)
	105 FORMAT(f4.1)
end do
 
 close(60)
 
 
! lecture du fichier entrees.txt et attribution de ses valeurs au tableau param
open (60,file='entrees.txt',action='read') ! on ouvre le fichier (,status='old')
nblignes=43 !nblignes=2+(NZ*3)+(dimension du vecteur B)=2+(12*3)+5=43
 
   do jj=1,nblignes
   ! on lit une ligne du fichier
   read(60,*) valeur
  ! on stocke la donnee lue dans la case correspondante du tableau   
  param(jj)=valeur
write (61, *) param(jj)
  end do
 close(60)
 
!!!Début mumps    
  call submump
  !close (61)
end program MUMPS_TEST
 
 
subroutine submump
implicit none
INCLUDE 'mpif.h'
INCLUDE 'dmumps_struc.h'
TYPE (DMUMPS_STRUC) mumps_par
integer IERR, II, Nb, NbZ,j,MZ
double precision       param(1:43)
 
 
CALL MPI_INIT(IERR)
!do k=1,5 !boucle dans le temps
 
 
!les composantes: COMM, SYM et PAR doivent être déf depuis que vous avez commencé etinis par l'utilisateur avant a call with JOB=-1
! Define a communicator for the package.
      mumps_par%COMM = MPI_COMM_WORLD !mis avant l'initialisation (JOB=-1)
!  Initialize an instance of the package
!  for L U factorization (sym = 0, with working host)
 
      mumps_par%JOB = -1 !initialise l'exemple par des valeurs par défaults (comme ICNTL) avant tout call mumps
      mumps_par%SYM = 0! matrice asymétrique
      mumps_par%PAR = 1 !PAR=1 obligatoire dans la version séquentielle
      CALL DMUMPS(mumps_par)
      !open (unit = 17,file = 'input_simpletest_real')
      IF (mumps_par%INFOG(1).LT.0) THEN !informations de tous les processeurs
       WRITE(6,'(A,A,I6,A,I9)') " ERROR RETURN: ","  mumps_par%INFOG(1)= ", mumps_par%INFOG(1), &
            "  mumps_par%INFOG(2)= ", mumps_par%INFOG(2) 
       GOTO 500
      END IF
!  Define problem on the host (processor 0)
	IF ( mumps_par%MYID .eq. 0 ) THEN
!définir une matrice centralized assembled avant la factorisation (Job=2) et l'analyse (job=1)
 
 
 
        mumps_par%N=Nb
        mumps_par%NZ=NbZ
        ALLOCATE( mumps_par%IRN ( NbZ ) )
        ALLOCATE( mumps_par%JCN ( NbZ ) )
        ALLOCATE( mumps_par%A( NbZ ) )
        ALLOCATE( mumps_par%RHS ( Nb ) )
 
		j=1
		DO II = 1, NbZ
		mumps_par%IRN(II)=param(j+2)
		mumps_par%JCN(II)=param(j+3)
		mumps_par%A(II)  =param(j+4)
		j=j+3
		end do
 
		MZ=(NbZ*3)+2
		DO II = 1, Nb
		mumps_par%RHS(II)=param(II+MZ)
		END DO
	END IF
 
 
!  Call package for solution
      mumps_par%JOB = 6!combines the actions of calls with JOB=1 (analyse), JOB=2(factorisation), and JOB=3(solution)
      CALL DMUMPS(mumps_par)
      IF (mumps_par%INFOG(1).LT.0) THEN
       WRITE(6,'(A,A,I6,A,I9)') " ERROR RETURN: ", "  mumps_par%INFOG(1)= ", mumps_par%INFOG(1), &
            "  mumps_par%INFOG(2)= ", mumps_par%INFOG(2) 
       GOTO 500
      END IF
!  Solution has been assembled on the host
      IF ( mumps_par%MYID .eq. 0 ) THEN
        WRITE(6, * ) ' Solution is ',(mumps_par%RHS(II),II=1,Nb) !afficher à l'écran la solution
      END IF
!  Deallocate user data
      IF ( mumps_par%MYID .eq. 0 )THEN
        DEALLOCATE( mumps_par%IRN )! can be deallocated after analysis
        DEALLOCATE( mumps_par%JCN )! can be deallocated after analysis
        DEALLOCATE( mumps_par%A   )
        DEALLOCATE( mumps_par%RHS )
      END IF
!  Destroy the instance (deallocate internal data structures)
      mumps_par%JOB = -2 !termine l'exemple: toutes les data de la structure sont déallocated (sauf pour celles mises par l'utilsateur): on le mets au dernier appel de Dmumps
      CALL DMUMPS(mumps_par)
      IF (mumps_par%INFOG(1).LT.0) THEN
       WRITE(6,'(A,A,I6,A,I9)') " ERROR RETURN: ", "  mumps_par%INFOG(1)= ", mumps_par%INFOG(1), &
            "  mumps_par%INFOG(2)= ", mumps_par%INFOG(2) 
       GOTO 500
      END IF
!end do !fin boucle temps
 500  CALL MPI_FINALIZE(IERR)
      STOP
  	!close (61)
END

des idées de solution?
Merci,

Bonjour,

En effet, les valeurs ne sont pas passées à la routine submump : la manière dont est écrit le programme, la routine submump est externe et ne voit pas les variables qui sont définies dans le program principal. Dans le cas précédent, les variables sont définies et utilisées au même endroit, donc ça fonctionne.

Les possibilités sont:
soit préférablement passer les variables utiles en argument de la routine submump(data1,etc...).
soit définir les variables dans un module contenant submump (possiblement program) du style (mais je déconseille fortement pour cette utilisation, la routine submump se veut générale j'ai l'impression et ça va la spécialiser trop fort et rendre à mon sens la lisibilité plus compliquée.)

Code:

1 2 3 4 5 6 7 8 module/program x declare and initialize data [...] contains subroutine usedata [...] end subroutine usedata end module/program
Après, je préfère aussi alerter sur le fait que la manière dont est implémenté le côté MPI du programme est dangereux.
En effet, le concept du parallélisme MPI n'est pas du fork & join, c'est vraiment un programme qui est lancé sur autant de machines à la fois qui vont communiquer (à la demande du code) pour savoir quoi faire. Les initialisations MPI_INIT et FINALIZE ne peuvent être appelées qu'une seule fois dans un programme implémentant des directives MPI, et elles ne délimitent pas une zone parallèle et séquentielle : toutes les unités suivent le programme en entier, la zone ainsi délimitée définit seulement à partir de quand et jusqu'où les unités vont devoir et pouvoir communiquer.
Appeler MPI_INIT et FINALIZE dans des routines séparées pour les spécialiser (par exemple avec des initialisations complémentaires liées aux rangs des machines à la répartition de charge ou autre), OK, mais là, c'est je pense à proscrire.
Ici: Mettre les INIT et FINALIZE dans une même sous-routine signifie que la sous-routine ne peut être appelée qu'une seule fois.
A part cela, il faut aussi voir qu'il va sans doute avoir un conflit dans ce qui est avant et après la sousroutine : toutes les unités vont éxécuter le programme dans sa totalité. Y compris l'écriture des fichiers qui va se faire en multiples (et/ou qui va planter je ne sais pas). Bref, quelques petits trucs à repenser.
-> Donc pour tout ce qui est parallélisme, ce qui est fait au moment l'est je pense dans une mauvaise direction même si vous pouvez arriver à faire in fine ce que vous voulez, ce n'est pas un usage canonique.

Cordialement,
xflr6

Merci xflr6 pour votre aide,
Le problème est résolu
:)