Discussion :

[djocin]

Bonjour

Dans le but de programmer l'algorithme du gradient conjugué pour résoudre un système du type Ax=b, je souhaite utiliser le produit matrice-vecteur sous forme additive, i.e que je ne veux pas écrire explicitement la matrice A et le vecteur x, pour ensuite en faire le produit, mais je veux écrire une subroutine qui me calcul directement le produit Ax.

Voici la subroutine que j'ai pour l'instant faite :

Code :

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

1
2
3
4
5
6
7
8
9
10
11
12
13
 
subroutine Mlt(u,v,nx,ny)
implicit none
real*8,dimension(:,:),allocatable::u,v
integer::i,j,nx,ny
 
do i=1,nx+1
   do j=1,ny+1
      v(i,j)=v(i,j)+u(i,j)
   enddo
enddo
 
end subroutine

Le vecteur d'entrée est u (écrit en fait sous forme d'une matrice), et la matrice de sortie est v (c'est donc Ax).

Le problème est que cette écriture ne marche pas vraiment !!

Avez-vous une idée s'il vous plaît ?

Merci.

[Ladgalen]

bonjour

D'un point de vue programmation il y a plusieurs problèmes :

* tes matrices u et v sont allocatables mais ne sont pas allouées.
* dans une subroutine (à moins qu'elle soit dans un module) les tableaux ne peuvent pas être dynamique (pas allocatable)
* tu écrits v = v + u mais v n'est pas initialisé

D'un point de vue algorithme je ne comprends pas ce que tu veux faire. D'autre part si tu veux juste résoudre A x = b, si A est inversible ce n'est pas la peine de faire une méthode du gradient conjugué. Tu est obligé ?

Bon courage

[genteur slayer]

Lorsque l'on cherche à résoudre un problème du type AX=B le produit AX représente le membre non constant d'une équation (ou d'un système)

donc si tu ne veut pas écrire A, mais AX, il te FAUT coder l'équation à la place (ce qui te fournira directement AX)

prenons l'exemple de l'équation de poisson: <laplacien>U=f on va dire en 3D et en différence finie pour faire simple du coups, l'équation discrétisé pour i,j,k au centre (il faudra faire évidement un traitement particulier pour les conditions de bords mais je n'en dirai rien aujourd'hui) on suppose dx, dy, dz les pas de discrétisations et U=(u, v ,w) f=(fx,fy,fz).
(u(i+1,j,k)-2*u(i,j,k)+u(i-1,j,k))/dx+(u(i,j+1,k)-2*u(i,j,k)+u(i,j+1,k))/dy+(u(i,j,k+1)-2*u(i,j,k)+u(i,j,k-1))/dz=fx
et pareil pour v et w
on définit donc AX comme un tableau à 4 dims: i,j,k et trois composantes d'où une fonction (tant qu'à faire):

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
 
function mlt(u,v,w,nx,ny,nz,dx,dy,dz)
  implicit none
  real (kind=8),dimension (1:,1:,1:),intent(in) :: u,v,w
  real (kind=8),intent(in) :: dx,dy,dz
  integer,intent(in) :: nx,ny,nz
  real(kind=8),dimension(1:nx,1:ny,1:nz,1:3) :: mlt
 integer :: i,j,k
  forall(i=2:nx-1,j=2:ny-1,z=2:nz-1) !attention, je ne traite pas les bords, il faut le faire à part
    mlt(i,j,k,1)=(u(i+1,j,k)-2.d0*u(i,j,k)+u(i-1,j,k))/dx+(u(i,j+1,k)-2.d0*u(i,j,k)+u(i,j+1,k))/dy+(u(i,j,k+1)-2.d0*u(i,j,k)+u(i,j,k-1))/dz
    mlt(i,j,k,2)=(v(i+1,j,k)-2.d0*v(i,j,k)+v(i-1,j,k))/dx+(v(i,j+1,k)-2.d0*v(i,j,k)+v(i,j+1,k))/dy+(v(i,j,k+1)-2.d0*v(i,j,k)+v(i,j,k-1))/dz
    mlt(i,j,k,3)=(w(i+1,j,k)-2.d0*w(i,j,k)+w(i-1,j,k))/dx+(w(i,j+1,k)-2.d0*w(i,j,k)+w(i,j+1,k))/dy+(w(i,j,k+1)-2.d0*w(i,j,k)+u(i,j,k-1))/dz
  end forall
  return
end function

dans cet exemple on est obligé de passé les taille de tableau dans les paramètre car il faut donner explicitement la taille du tableau de sortie mlt, sinon, l'utilisation de la fonction size permet d'obtenir les tailles de tes tableau: nx=size(u,1);ny=size(u,2);nz=size(u,3)

[Sylvain Bergeron]

...on est obligé de passé les taille de tableau dans les paramètre car il faut donner explicitement la taille du tableau de sortie mlt, sinon, l'utilisation de la fonction size permet d'obtenir les tailles de tes tableau: nx=size(u,1);ny=size(u,2);nz=size(u,3)

Il est permis d'utiliser des fonctions d'information dans les déclarations, de sorte que le code suivant est légal :

Code :

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

1
2
3
4
5
6
7
8
 
function mlt(u,v,w,dx,dy,dz)
  implicit none
  real (kind=8),dimension (1:,1:,1:),intent(in) :: u,v,w
  real (kind=8),intent(in) :: dx,dy,dz
  real(kind=8),dimension(1:size(u,1),1:size(u,2),1:size(u,3),1:3) :: mlt
 integer :: i,j,k
...

[genteur slayer]

ah oui? et bien c'est encore mieux: moins de paramètres d'entrée!!!