Discussion :

[lili2002]

Bonjour,
Je veux paralléliser la methode du gradient conjugué en mpi avec fortran, mais je ne sais pas comment proceder, je suis debutant en mpi.
J'ai déja écrit le code sequentiel en fortran 90.

Merci de bien vouloir m'aider a ecrire le code parallele.

[Jatukam]

http://www.idris.fr/data/cours/paral...choix_doc.html

[lili2002]

Merci, je suis déja allé sur ce site. Je n'ai pas trouvé un truc qui peut m'aider.

Voici le code que j'ai essayé de faire mas ça ne donne pas la bonne solution.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
!--------------------------------------------------------------------------
!programme parallele : méthode de gradient conjugué                       !
!--------------------------------------------------------------------------
program gradconjpara
 
 
!######################################################################
 
 
 
!-----------------------Déclaration des variables--------------------------
 
   implicit none
   include "mpif.h"
 
!----------------------Déclaration des entiers réelles---------------------
 
   real*8 :: eps,alpha,beta,dx,dy,dt,Lx,Ly,D1,T,x,y,error,c,c1,c2
   integer :: i,j,m,nx,ny,nt,n,l,ns
 
 
!----------------------Déclaration des tableaux---------------------------
 
   real*8, dimension(:  , :    )   , allocatable ::A,A1
   real*8, dimension(  :   ), allocatable :: k,B,r,r1,d,w
 
!--------------------------mpiv------------------------------------------ 
!--------------------Déclaration variables parallèle----------------------
 
   integer :: me,statinfo,np,nyloc,i1,in,dim,dim1,j1,jn
   real*8  :: somme,somme1,somme2
   integer, parameter :: Msgtag=50
   integer, dimension(MPI_STATUS_SIZE) :: status
 
 
 
!----------------------Lecture fichier données----------------------------
 
   open(11,file='data',status='old',form='formatted')
 
    read(11,*) nx    
    read(11,*) ny     
    read(11,*) nt        
    read(11,*)Lx  
    read(11,*)Ly
    read(11,*)T   
    read(11,*)D1    
    read(11,*)eps   
    read(11,*)m
 
 
!-------------------------Maillage---------------------------------------
 
   dx= real(Lx)/(nx+1)
   dy= real(Ly)/(ny+1)
   dt=real(T)/(nt+1) 
   c= (2.*D1)/(dx*dx) + (2.*D1)/(dy*dy) 
   c1=D1/(dx*dx) 
   c2= D1/(dy*dy)
   ny=4
   nx=5
   ns=nx*ny
 
!----------------------main-----------------------------------------------
!----------------------mpi-----------------------------------------------
 
    call  MPI_INIT(statinfo)
    call  MPI_COMM_SIZE(MPI_COMM_WORLD,np,statinfo)
    call  MPI_COMM_RANK(MPI_COMM_WORLD, me, statinfo)
    call  repartitioncharge(me, ns, np, i1, in)
    !print*,i1,in
!--------------------Allocation dynamique des tableaux-------------------   
 
    allocate(A(1:ns,1:ns))
    allocate(B(1:ns))
    allocate(k(1:ns))
    allocate(r(1:ns)) 
    allocate(d(1:ns))
    allocate(w(1:ns))
    allocate(r1(1:ns))
    allocate(A1(i1:in,i1:in+nx))
 
 
 
 
!--------------------Initialisation----------------------------------------
   if (np/=ny ) stop 'Numbers of processors must be nx!!!' 
 
   A=0.
   k=0.
   B=1.
 
!----------------------Appel subroutine-----------------------------------
 
   call    mat(c,c1,c2,nx,ny,A)
   call    gradientconj(i1,in,nx,ny,m,eps,b,k,A)
   do j= 0,np-1
       call  repartitioncharge(j, ns, np, j1, jn)
       dim1=jn-j1+1
       call MPI_BCAST(k(j1:jn),dim1, MPI_real8 ,j, MPI_COMM_WORLD ,statinfo)
   end do
   if (me==0) then
     print*,k
   end if
!--------------------main-------------------------------------------------
 
!--------------------Désallouement dynamique------------------------------
 
    deallocate(A)
    deallocate(b)
    deallocate(k)
    deallocate(r) 
    deallocate(d)
    deallocate(w)
    deallocate(r1)  
 
 call  MPI_FINALIZE( statinfo)
   contains  
!--------------------Définition de la matrice A----------------------------
 
 
     subroutine mat(c,c1,c2,nx,ny,A)
        implicit none
        real*8, dimension(1:ny*nx,1:ny*nx), intent(inout) :: A
        real*8, intent(in) :: c,c1,c2
        integer, intent(in):: nx, ny
        integer ::i,j   
 
 
        do i=1,nx*ny
              a(i,i)=c
 
 
            if( i .gt. Nx ) a(i,i-Nx) = -c2
 
 
            if( i .lt. (nx*ny-nx+1) )  a(i,i+nx) = -c2
 
 
        end do
 
        do i=2,nx*ny
            a(i,i-1)=-c1
        end do
 
        do i=1, nx*ny-1
 
           a(i,i+1)=-c1
        end do
 
        do i=1,ny-1
          a(i*nx+1,i*nx)=0
          a(i*nx,i*nx +1)=0
        end do   
 
    end subroutine mat
 
 
 
!--------------------------------------------------------------------------
!subroutine :parallélisation de la  méthode de gradient conjugué          !
!--------------------------------------------------------------------------
 
 
   subroutine gradientconj(i1,in,nx,ny,m,eps,b,k,A)
        implicit none
        integer, intent(in) ::m,i1,in,nx,ny
        integer :: ns, l,i,j1,jn,dim,dim1 
        real*8, dimension(1:nx*ny,1:nx*ny), intent(in) :: A
        real*8, dimension(1:nx*ny), intent(inout) :: k
        real*8, dimension(1:nx*ny), intent(in) :: b       
        real*8, intent(in) :: eps
        real*8, dimension(1:nx*ny) :: d,r,w
        real*8 :: alpha, beta,somme
 
        ns=ny*nx    
        l=0
        if (i1==1)then
           r(i1:in)=matmul(A(i1:in,i1:in+nx),k(i1:in+nx))-B(i1:in)
        elseif (in==ns) then
           r(i1:in)=matmul(A(i1:in,i1-nx:in),k(i1-nx:in))-B(i1:in)
        else
           r(i1:in)=matmul(A(i1:in,i1-nx:in+nx),k(i1-nx:in+nx))-B(i1:in)
        end if
 
        d(i1:in)=r(i1:in)
 
 
        do while ( 0<=l .and. l<=m .and. sqrt(dot_product(r,r))>eps)
              do j= 0,np-1
                call  repartitioncharge(j, ns, np, j1, jn)
                dim1=jn-j1+1
                call MPI_BCAST(d(j1:jn),dim1, MPI_real8 ,j, MPI_COMM_WORLD ,statinfo)
             end do
 
             if (i1==1)then
                w(i1:in)=matmul(A(i1:in,i1:in+nx),d(i1:in+nx))
             elseif (in==ns) then
                w(i1:in)=matmul(A(i1:in,i1-nx:in),d(i1-nx:in))
             else
9                w(i1:in)=matmul(A(i1:in,i1-nx:in+nx),d(i1-nx:in+nx))
             end if
 
 
 
             alpha=(dot_product(d(i1:in),r(i1:in)))/(dot_product(d(i1:in),w(i1:in)))
             call MPI_ALLREDUCE (alpha,somme,1, MPI_REAL , MPI_SUM , MPI_COMM_WORLD ,statinfo)
 
             k(i1:in)=k(i1:in)-alpha*d(i1:in)
 
 
             beta=1./(dot_product(r(i1:in),r(i1:in)))
             call MPI_ALLREDUCE (beta,somme,1, MPI_REAL8 , MPI_SUM , MPI_COMM_WORLD ,statinfo)
 
             r(i1:in)=r(i1:in)-alpha*w(i1:in)
 
             beta=beta*(dot_product(r(i1:in),r(i1:in)))
             call MPI_ALLREDUCE (beta, somme,1, MPI_REAL8 , MPI_SUM , MPI_COMM_WORLD ,statinfo)
 
             d(i1:in)=r(i1:in)+beta*d(i1:in)
 
             l=l+1
 
        end do
   end subroutine gradientconj
!---------------------------------charge-------------------------------------------------------------------
 
   subroutine repartitioncharge(me, ns, np, i1, in)
          implicit none
 
         integer,  intent(in) ::  ns, np, me
         integer,  intent(out) :: i1, in
         integer :: r
 
         r=ns- np*(ns/np) 
 
 
         if (me<r) then
            i1=me*(ns/np+1) +1
            in = (me+1)*(ns/np+1)+1 
 
         else
             i1=me*ns/np +r+1
             in = i1 + ns/np-1 
         end if
       end subroutine repartitioncharge
 
end program gradconjpara

fichier data

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
3   !nx
3   !ny
6   !nt
1.0 !Lx
1.0 !Lx
1.0 !T
1.0 !D1
1.d-10 !eps
10000 !m

la solution doit varier entre 0 et 1.
Merci de bien vouloir m'aider!!!!!

[plumedesiles]

Salut,

Tu peux essayer sur google scholar avec conjugate gradient mpi, tu trouveras surement quelque chose.

Quelle version du CG utilises tu ? Steepest Descent
Une ref compléte sur le suejt http://www.cs.cmu.edu/~quake-papers/...e-gradient.pdf

La parallélisation du gradient conjugué repose essentiellement sur la parallélisation du produit matrice vecteur (creux généralement).

Chacun de tes process va devoir calculer une ligne du produit matrice vecteur, puis tu peux faire un gather à la fin pour calculer la norme et la faire vérifier par rapport à ton erreur.
Voilà comment en gros je procéderais, je me suis jamais réllement penché sur le problème.

Peux tu donner la version de l'algo que tu veux implémenté ? On peut y réfléchir ensemble et je peux te proposer une solution

[lili2002]

Bonjour,
Merci beaucoup de ton aide.

Voici le code sequentiel (gradient conjugué) en fortran 90 que je veux paralléliser(MPI).

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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
!--------------------------------------------------------------------------
!subroutine : méthode de gradient conjugué                                 !
!--------------------------------------------------------------------------
 
 
   subroutine gradientconj(nx,ny,m,eps,b,k,A)
        implicit none
        real*8, dimension(1:nx*ny,1:nx*ny), intent(in) :: A
        real*8, dimension(1:nx*ny), intent(inout) :: k
        real*8, dimension(1:nx*ny), intent(in) :: b       
        real*8, intent(in) :: eps
        integer, intent(in) :: nx,ny,m
        real*8, dimension(1:nx*ny) :: d,r,w,r1
        real*8 :: alpha, beta
        integer :: l    
 
        l=0
        r=matmul(A,k)-b
        d=r
        do while ( 0<=l .and. l<=m .and. sqrt(dot_product(r,r))>eps)
            w=(matmul(A,d))
            alpha=(dot_product(d,r))/(dot_product(d,w))
            k=k-alpha*d
            r1=r-alpha*w
            beta=(dot_product(r1,r1))/(dot_product(r,r))
            d=r1+beta*d
            l=l+1
            r=r1               
        end do
   end subroutine gradientconj
end program gradient

[plumedesiles]

Il faut paralléliser tes matmul et tes dot_product

Première chose répartir les lignes de ta matrice A sur chaque processeurs.
Par exemple 10 lignes par processeur.
Chaque processeur sera en charge de calculer le résultat pour ces 10 lignes
Pour tes matmul qui sont en faite des produits matrice vecteur (qui doivent être creux tu peux confirmer?) Ax = y
Chaque processeur calcule 10 lignes de A, le vecteur x existe sur chaque processeur pour éviter les communications. De toute manière c'est pas le plus gros consommateur en mémoire.
Tu te retrouves donc à faire sur chaque processeur un matmul, en fait un produit matrice vecteur, de petite dimension => Aij * Xj = Yi
Ton résultat Y est donc répartir par groupe de 10lignes sur chaque processeur

As tu compris celà et peux ton aller plus loin ou dois je être plus clair ?
Si tu as compris j'explique la prochaine fois le dot_product comment on peut procéder.

[lili2002]

oui j' ai compris merci.
Expliques moi le dot-product parce que quand je l'ajoute la solution me donne NaN.
A plus