| 1 | ! Transformee de Fourier reelle-complexe 3d |
|---|
| 2 | |
|---|
| 3 | program tjmscfft3d |
|---|
| 4 | |
|---|
| 5 | implicit none |
|---|
| 6 | |
|---|
| 7 | integer, parameter :: n = 8 |
|---|
| 8 | integer, parameter :: m = 16 |
|---|
| 9 | integer, parameter :: l = 4 |
|---|
| 10 | real, dimension(0:n-1,0:m-1,0:l-1) :: x, xx |
|---|
| 11 | complex, dimension(0:n/2,0:m-1,0:l-1) :: y |
|---|
| 12 | |
|---|
| 13 | integer, parameter :: ntable = 100+2*(n+m+l) |
|---|
| 14 | real, dimension(0:ntable-1) :: table |
|---|
| 15 | |
|---|
| 16 | ! Les routines jm ont besoin de 2*2*(n/2+1)*m*l. D'ou tronconnage. |
|---|
| 17 | integer, parameter :: nwork = 512*max(n,m,l) |
|---|
| 18 | real, dimension(0:nwork-1) :: work |
|---|
| 19 | |
|---|
| 20 | integer :: isign |
|---|
| 21 | real :: scale |
|---|
| 22 | integer :: isys |
|---|
| 23 | integer :: i, j, k |
|---|
| 24 | integer :: iz, j1, k1 |
|---|
| 25 | real :: twopi |
|---|
| 26 | complex :: s |
|---|
| 27 | |
|---|
| 28 | twopi = 2 * acos(real(-1)) |
|---|
| 29 | |
|---|
| 30 | ! On prepare le tableau d'entree |
|---|
| 31 | call random_number( x ) |
|---|
| 32 | xx = x |
|---|
| 33 | |
|---|
| 34 | isys = 0 |
|---|
| 35 | scale = 1./real(n) |
|---|
| 36 | |
|---|
| 37 | isign = 0 |
|---|
| 38 | call scfft3d(isign,n,m,l,scale,x,n,m,y,n/2+1,m,table,work,isys) |
|---|
| 39 | isign = 1 |
|---|
| 40 | print *,'jmscfft3d ',n,m,l,isign,scale |
|---|
| 41 | call scfft3d(isign,n,m,l,scale,x,n,m,y,n/2+1,m,table,work,isys) |
|---|
| 42 | |
|---|
| 43 | ! On imprime le tableau de sortie |
|---|
| 44 | open(10,file='temp1',status='unknown',form='formatted') |
|---|
| 45 | write(10,'(2e25.12)') y |
|---|
| 46 | |
|---|
| 47 | ! Ce qu'il faut trouver |
|---|
| 48 | open(11,file='temp2',status='unknown',form='formatted') |
|---|
| 49 | ! On reprepare le tableau d'entree |
|---|
| 50 | x = xx |
|---|
| 51 | ! Et on calcule |
|---|
| 52 | do k = 0,l-1 |
|---|
| 53 | do j = 0,m-1 |
|---|
| 54 | do i = 0,n/2 |
|---|
| 55 | s = 0 |
|---|
| 56 | do k1 = 0,l-1 |
|---|
| 57 | do j1 = 0,m-1 |
|---|
| 58 | do iz = 0,n-1 |
|---|
| 59 | s = s + cmplx( & |
|---|
| 60 | & cos(twopi*i*iz/real(n)), & |
|---|
| 61 | & isign*sin(twopi*i*iz/real(n))) & |
|---|
| 62 | & * cmplx( & |
|---|
| 63 | & cos(twopi*j*j1/real(m)), & |
|---|
| 64 | & isign*sin(twopi*j*j1/real(m))) & |
|---|
| 65 | & * cmplx( & |
|---|
| 66 | & cos(twopi*k*k1/real(l)), & |
|---|
| 67 | & isign*sin(twopi*k*k1/real(l))) & |
|---|
| 68 | & * x(iz,j1,k1) |
|---|
| 69 | end do |
|---|
| 70 | end do |
|---|
| 71 | end do |
|---|
| 72 | write(11,'(2e25.12)') s*scale |
|---|
| 73 | end do |
|---|
| 74 | end do |
|---|
| 75 | end do |
|---|
| 76 | |
|---|
| 77 | end program tjmscfft3d |
|---|
