Discussion :

[Johnson25]

bonjour; avant toute chose le programme que je vais donner fonctionne parfaitement. Il calcule la transformée Fourier discrète en deux dimensions. Je l'ai testé avec une matrice 4x4 unitaire. Cela me donne la matrice :

Code :

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

1
2
3
$ ./testforonematrix
              (0.99999994039535522,0.0000000000000000)               (0.0000000000000000,0.0000000000000000)
               (0.0000000000000000,0.0000000000000000)               (0.0000000000000000,0.0000000000000000)

Normalement cela devrait être un 1.000000000 en première composante ( j'ai donc une erreur de 10^-7 ). Au début je pensais que le problème venait de la simple précision du programme, j'ai donc changé en double précision. Pourtant j'ai toujours une erreur de 10^-7.
D'où pourrait venir l'erreur?
En vous remerciant de votre aide future

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
 
module mdfourier
 
implicit none
 
contains
 
  subroutine ft2d(n1,n2,cp,isig)
    implicit none
    integer :: n1,n2,isig,i1,i2
    complex*8 :: cp(n1,n2),cw(n2)
    external fork
 
    do i2= 1,n2
      call fork(n1,cp(1,i2),isig)
    end do
 
    do i1= 1,n1
      do i2= 1,n2
        cw(i2)= cp(i1,i2)
      end do
      call fork(n2,cw,isig)
      do i2= 1,n2
        cp(i1,i2)= cw(i2)
      end do
    end do
   do i1=1,n1
print*,((1.D0/n1)*cp(:,i1))
end do
    return
  end subroutine ft2d
 
end module mdfourier
 
  subroutine fork(lx,cx,isig)
 
    implicit none
    integer :: lx,isig,i,j,l,m,istep
    complex*8 :: cx(lx),carg,cexp,cw,ctemp
    real*8 :: pii,sc
    pii= 3.141592653589793
 
    j= 1
    sc= sqrt(1./lx)
    do i= 1,lx
      if(i <= j) then
        ctemp= cx(j)*sc
        cx(j)= cx(i)*sc
        cx(i)= ctemp
      end if
      m= lx/2
      do while (m >= 1 .and. j > m)
        j= j-m
        m= m/2
      end do
      j= j+m
    end do
    l= 1
    do while (l < lx)
      istep= 2*l
      do m= 1,l
        carg= cmplx(0.,1.)*(pii*isig*(m-1))/l
        cw= cexp(carg)
        do i= m,lx,istep
          ctemp= cw*cx(i+l)
          cx(i+l)= cx(i)-ctemp
          cx(i)= cx(i)+ctemp
        end do
      end do
      l= istep
    end do
    return
  end subroutine fork
 
 
 
 
 
 
 
 
PROGRAM fourier
  use mdfourier
  IMPLICIT NONE
  complex*8, dimension(2,2) :: A
  complex*8 :: II=(0.D0,1.D0)
  integer :: i,j 
 
do j=1,2  
do i=1,2
A(i,j)=1.D0;
end do 
end do 
 
 
call ft2d(2,2,A,-1)
 
 
END PROGRAM fourier

[Sylvain Bergeron]

Il y a au moins une erreur ici :

Code :

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

1
2
3
 
real*8 :: pii
pii= 3.141592653589793

pii est bien real 8, mais tu lui attribues une constante real 4. Quand bien même tu lui ajouterais 100 décimales, c'est une constante real 4. Pour que ça fonctionne, tu dois dire explicitement que la constante est real 8 : 3.141592653589793d0 ou 3.141592653589793_8 ou par la méthode standard avec une constante comme rkind8 = selected_real... et pii = 3.141592653589793_rkind8