dox/html/fftmod_8f90_source.html

00001 !     =============
00002 !     MODULE FFTMOD
00003 !     =============
00004
00005       module fftmod
00006       parameter(NRES = 12)
00007       integer :: nallowed(NRES) = (/16,32,48,64,96,128,256,384,512,1024,2048,4096/)
00008 !     T3    - N16   : 8-2
00009 !     T10   - N32   : 8-2-2
00010 !     T15   - N48   : 8-3-2
00011 !     T21   - N64   : 8-4-2
00012 !     T31   - N96   : 8-4-3
00013 !     T42   - N128  : 8-4-4
00014 !     T85   - N256  : 8-4-4-2
00015 !     T127  - N384  : 8-4-4-3
00016 !     T170  - N512  : 8-4-4-4
00017 !     T341  - N1024 : 8-4-4-4-2
00018 !     T682  - N2048 : 8-4-4-4-4
00019 !     T1365 - N4096 : 8-4-4-4-4-2
00020
00021       integer :: lastn = 0
00022       real,allocatable :: trigs(:)
00023       end module fftmod
00024
00025 !     ================
00026 !     SUBROUTINE GP2FC
00027 !     ================
00028
00029       subroutine gp2fc(a,n,lot)
00030       use fftmod
00031       real a(n,lot)
00032
00033       if (n /= lastn) then
00034          if (allocated(trigs)) deallocate(trigs)
00035          allocate(trigs(n))
00036          lastn = n
00037          call fftini(n)
00038       endif
00039
00040       call dfft8(a,a,n,lot)
00041       la = n / 8
00042       do while (la >= 4)
00043          call dfft4(a,trigs,n,lot,la)
00044       enddo
00045
00046       if (la == 3) then
00047          do l = 1 , lot
00048             call dfft3(a(1,l),trigs,n)
00049          enddo
00050       endif
00051
00052       if (la == 2) then
00053          do l = 1 , lot
00054             call dfft2(a(1,l),trigs,n)
00055          enddo
00056       endif
00057       return
00058       end subroutine gp2fc
00059
00060 !     ================
00061 !     SUBROUTINE FC2GP
00062 !     ================
00063
00064       subroutine fc2gp(a,n,lot)
00065       use fftmod
00066       real a(n,lot)
00067
00068       if (n /= lastn) then
00069          if (allocated(trigs)) deallocate(trigs)
00070          allocate(trigs(n))
00071          lastn = n
00072          call fftini(n)
00073       endif
00074
00075       nf = n/8
00076       do while (nf >= 4)
00077          nf = nf/4
00078       enddo
00079       la = 1
00080       if (nf == 2) call ifft2(a,trigs,n,lot,la)
00081       if (nf == 3) call ifft3(a,trigs,n,lot,la)
00082       do while (la < n/8)
00083          call ifft4(a,trigs,n,lot,la)
00084       enddo
00085       call ifft8(a,a,n,lot)
00086       return
00087       end subroutine fc2gp
00088
00089 !     =================
00090 !     SUBROUTINE FFTINI
00091 !     =================
00092
00093       subroutine fftini(n)
00094       use fftmod
00095       logical labort
00096
00097 !     check for allowed values of n
00098
00099       labort = .true.
00100       do j = 1 , NRES
00101          if (n == nallowed(j)) labort = .false.
00102       enddo
00103
00104       if (labort) then
00105          write (*,*) '*** FFT does not support n = ',n,' ***'
00106          write (*,*) 'Following resolutions may be used:'
00107          write (*,*) '----------------------------------'
00108          do j = 1 , NRES
00109             write (*,1000) nallowed(j), nallowed(j)/2, nallowed(j)/3
00110          enddo
00111          stop
00112       endif
00113  1000 format(' NLON=',I5,'  NLAT=',I5,'  NTRU=',I5)
00114
00115       del = 4.0 * asin(1.0) / n
00116       do k=0,n/2-1
00117         angle = k * del
00118         trigs(2*k+1) = cos(angle)
00119         trigs(2*k+2) = sin(angle)
00120       enddo
00121       return
00122       end subroutine fftini
00123
00124 !     ================
00125 !     SUBROUTINE DFFT2
00126 !     ================
00127
00128       subroutine dfft2(a,trigs,n)
00129       dimension a(n),c(n),trigs(n)
00130
00131       c(1) = a(1) + a(2)
00132       c(2) = 0.0
00133
00134       ja = 3
00135       jb = n - 1
00136
00137       do i=3,n-5,4
00138          c1 = trigs(ja  )
00139          s1 = trigs(ja+1)
00140          a1p3 = c1 * a(i+1) + s1 * a(i+3)
00141          a3m1 = c1 * a(i+3) - s1 * a(i+1)
00142          c(ja  ) = a(i) + a1p3
00143          c(jb  ) = a(i) - a1p3
00144          c(ja+1) = a3m1 + a(i+2)
00145          c(jb+1) = a3m1 - a(i+2)
00146          ja = ja + 2
00147          jb = jb - 2
00148       enddo
00149
00150       c(ja  ) =  a(n-1)
00151       c(ja+1) = -a(n  )
00152
00153       a = c
00154       return
00155       end subroutine dfft2
00156
00157 !     ================
00158 !     SUBROUTINE DFFT3
00159 !     ================
00160
00161       subroutine dfft3(a,trigs,n)
00162       parameter(SIN60 = 0.866025403784438D0)
00163       dimension a(n),c(n),trigs(n)
00164
00165       ja = 1              !  1
00166       jb = 2 * (n/3)  + 1 ! 65
00167       jc = jb             ! 65
00168
00169       c(ja  ) = a(1) + a(2) + a(3)
00170       c(ja+1) = 0.0
00171       c(jb  ) = a(1) - 0.5 * (a(2) + a(3))
00172       c(jb+1) =      SIN60 * (a(3) - a(2))
00173
00174       ja = 3         !  3, 5, 7, ... ,31
00175       jb = jb + 2    ! 67,69,71, ... ,95
00176       jc = jc - 2    ! 63,61,59, ... ,35
00177
00178       do i = 4 , n-8 , 6 ! 88
00179          c1 = trigs(ja  )
00180          s1 = trigs(ja+1)
00181          c2 = trigs(ja+ja-1)
00182          s2 = trigs(ja+ja  )
00183          a1 = (c1*a(i+1)+s1*a(i+4))+(c2*a(i+2)+s2*a(i+5))
00184          b1 = (c1*a(i+4)-s1*a(i+1))+(c2*a(i+5)-s2*a(i+2))
00185          a2 = a(i  ) - 0.5 * a1
00186          b2 = a(i+3) - 0.5 * b1
00187          a3 = SIN60*((c1*a(i+1)+s1*a(i+4))-(c2*a(i+2)+s2*a(i+5)))
00188          b3 = SIN60*((c1*a(i+4)-s1*A(i+1))-(c2*a(i+5)-s2*a(i+2)))
00189          c(ja  ) = a(i  ) + a1
00190          c(ja+1) = a(i+3) + b1
00191          c(jb  ) = a2 + b3
00192          c(jb+1) = b2 - a3
00193          c(jc  ) = a2 - b3
00194          c(jc+1) =-b2 - a3
00195          ja = ja + 2
00196          jb = jb + 2
00197          jc = jc - 2
00198       enddo
00199
00200       if (ja <= jc) then ! ja=33  jc=33
00201          c(ja  ) = a(n-2) + 0.5 * (a(n-1) - a(n)) ! 33
00202          c(ja+1) =       -SIN60 * (a(n-1) + a(n)) ! 34
00203       endif
00204       a(:) = c(:)
00205       return
00206       end subroutine dfft3
00207
00208 !     ================
00209 !     SUBROUTINE DFFT4
00210 !     ================
00211
00212       subroutine dfft4(a,trigs,n,lot,la)
00213       dimension a(n,lot),c(n,lot),trigs(n)
00214       la = la / 4
00215
00216       i1 = la
00217       i2 = la + i1
00218       i3 = la + i2
00219       i4 = la + i3
00220       i5 = la + i4
00221       i6 = la + i5
00222       i7 = la + i6
00223
00224       j1 = n/2 - la
00225       j2 = n - la
00226       j3 = j1
00227       j5 = j1 + la
00228
00229       do i=1,la
00230          do l=1,lot
00231          a0p2 = a(i   ,l) + a(i2+i,l)
00232          a1p3 = a(i1+i,l) + a(i3+i,l)
00233          c(   i,l) = a0p2 + a1p3
00234          c(j2+i,l) = a0p2 - a1p3
00235          c(j1+i,l) = a(   i,l) - a(i2+i,l)
00236          c(j5+i,l) = a(i3+i,l) - a(i1+i,l)
00237          enddo
00238       enddo
00239
00240       jink = 2 * la
00241       j0 = la
00242       j1 = j1 + jink
00243       j2 = j2 - jink
00244       j3 = j3 - jink
00245       j4 = j0 + la
00246       j5 = j1 + la
00247       j6 = j2 + la
00248       j7 = j3 + la
00249
00250       ibase=4*la
00251
00252       do 450 k=la,(n-4)/8,la
00253          kb=k+k
00254          kc=kb+kb
00255          kd=kc+kb
00256          c1=trigs(kb+1)
00257          s1=trigs(kb+2)
00258          c2=trigs(kc+1)
00259          s2=trigs(kc+2)
00260          c3=trigs(kd+1)
00261          s3=trigs(kd+2)
00262
00263          i=ibase+1
00264          do j=1,la
00265             do l=1,lot
00266             a1p5 = c1 * a(i1+i,l) + s1 * a(i5+i,l)
00267             a2p6 = c2 * a(i2+i,l) + s2 * a(i6+i,l)
00268             a3p7 = c3 * a(i3+i,l) + s3 * a(i7+i,l)
00269             a5m1 = c1 * a(i5+i,l) - s1 * a(i1+i,l)
00270             a6m2 = c2 * a(i6+i,l) - s2 * a(i2+i,l)
00271             a7m3 = c3 * a(i7+i,l) - s3 * a(i3+i,l)
00272             a0 = a(i,l) + a2p6
00273             a2 = a(i,l) - a2p6
00274             a1 = a1p5 + a3p7
00275             a3 = a3p7 - a1p5
00276             b0 = a(i4+i,l) + a6m2
00277             b2 = a(i4+i,l) - a6m2
00278             b1 = a5m1 + a7m3
00279             b3 = a5m1 - a7m3
00280             c(j0+j,l) = a0+a1
00281             c(j2+j,l) = a0-a1
00282             c(j4+j,l) = b0+b1
00283             c(j6+j,l) = b1-b0
00284             c(j1+j,l) = a2+b3
00285             c(j3+j,l) = a2-b3
00286             c(j5+j,l) = a3+b2
00287             c(j7+j,l) = a3-b2
00288             enddo
00289             i=i+1
00290          enddo
00291
00292         ibase=ibase+8*la
00293         j0 = j0 + jink
00294         j1 = j1 + jink
00295         j2 = j2 - jink
00296         j3 = j3 - jink
00297         j4 = j0 + la
00298         j5 = j1 + la
00299         j6 = j2 + la
00300         j7 = j3 + la
00301 450   continue
00302       if (j1 <= j2) then
00303          sin45=sqrt(0.5)
00304          i=ibase+1
00305          do j=1,la
00306             do l=1,lot
00307             a1p3 = sin45 * (a(i1+i,l) + a(i3+i,l))
00308             a1m3 = sin45 * (a(i1+i,l) - a(i3+i,l))
00309             c(j0+j,l) =  a(   i,l) + a1m3
00310             c(j1+j,l) =  a(   i,l) - a1m3
00311             c(j4+j,l) = -a(i2+i,l) - a1p3
00312             c(j5+j,l) =  a(i2+i,l) - a1p3
00313             enddo
00314             i=i+1
00315          enddo
00316       endif
00317       if (la == 1) then
00318          do l=1,lot
00319          a(1,l) = c(1,l)
00320          a(2,l) = 0.0
00321          a(3:n,l) = c(2:n-1,l)
00322          enddo
00323       else
00324          a = c
00325       endif
00326       return
00327       end subroutine dfft4
00328
00329 !     ================
00330 !     SUBROUTINE DFFT8
00331 !     ================
00332
00333       subroutine dfft8(a,c,n,lot)
00334       real a(n*lot),c(n*lot)
00335       la = n / 8
00336       z  = 1.0 / n
00337       zsin45 = z * sqrt(0.5)
00338
00339       do i=0,la*lot-1
00340          i0 = (i/la) * n + mod(i,la) + 1
00341          i1 = i0 + la
00342          i2 = i1 + la
00343          i3 = i2 + la
00344          i4 = i3 + la
00345          i5 = i4 + la
00346          i6 = i5 + la
00347          i7 = i6 + la
00348
00349          a0p4 =  a(i0) + a(i4)
00350          a1p5 =  a(i1) + a(i5)
00351          a2p6 =  a(i2) + a(i6)
00352          a3p7 =  a(i3) + a(i7)
00353          a5m1 =  a(i5) - a(i1)
00354          a7m3 =  a(i7) - a(i3)
00355          a0m4 = (a(i0) - a(i4)) * z
00356          a6m2 = (a(i6) - a(i2)) * z
00357
00358          a0p4p2p6 = a0p4 + a2p6
00359          a1p5p3p7 = a1p5 + a3p7
00360          a7m3p5m1 = (a7m3 + a5m1) * zsin45
00361          a7m3m5m1 = (a7m3 - a5m1) * zsin45
00362
00363          c(i0) = z * (a0p4p2p6 + a1p5p3p7)
00364          c(i7) = z * (a0p4p2p6 - a1p5p3p7)
00365          c(i3) = z * (a0p4 - a2p6)
00366          c(i4) = z * (a3p7 - a1p5)
00367          c(i1) = a0m4 + a7m3m5m1
00368          c(i5) = a0m4 - a7m3m5m1
00369          c(i2) = a7m3p5m1 + a6m2
00370          c(i6) = a7m3p5m1 - a6m2
00371       enddo
00372       return
00373       end subroutine dfft8
00374
00375 !     ================
00376 !     SUBROUTINE IFFT4
00377 !     ================
00378
00379       subroutine ifft4(c,trigs,n,lot,la)
00380       dimension a(n,lot),c(n,lot),trigs(n)
00381
00382       if (la == 1) then
00383          a(1,:) = 0.5 * c(1,:)
00384          a(n,:) = 0.0
00385          a(2:n-1,:) = c(3:n,:)
00386       else
00387          a = c
00388       endif
00389
00390       m=n/4
00391       kstop=(n-4)/8
00392
00393       i1 = n/2 - la
00394       i2 = n   - la
00395       i5 = i1  + la
00396
00397       j1 = la
00398       j2 = la+j1
00399       j3 = la+j2
00400       j4 = la+j3
00401       j5 = la+j4
00402       j6 = la+j5
00403       j7 = la+j6
00404
00405       do i=1,la
00406       do l=1,lot
00407          c(   i,l) = a(i,l) + a(i2+i,l) + a(i1+i,l)
00408          c(j1+i,l) = a(i,l) - a(i2+i,l) - a(i5+i,l)
00409          c(j2+i,l) = a(i,l) + a(i2+i,l) - a(i1+i,l)
00410          c(j3+i,l) = a(i,l) - a(i2+i,l) + a(i5+i,l)
00411       enddo
00412       enddo
00413
00414       iink  = 2 * la
00415       jbase = 4 * la + 1
00416       i0    = la
00417       i1    = i0 + n/2
00418       i2    = n - 3 * la
00419       i3    = i2 - n/2
00420       i4    = i0 + la
00421       i5    = i1 + la
00422       i6    = i2 + la
00423       i7    = i3 + la
00424
00425       do 450 k=la,kstop,la
00426         kb=k+k
00427         kc=kb+kb
00428         kd=kc+kb
00429         c1=trigs(kb+1)
00430         s1=trigs(kb+2)
00431         c2=trigs(kc+1)
00432         s2=trigs(kc+2)
00433         c3=trigs(kd+1)
00434         s3=trigs(kd+2)
00435         do i = 1 , la
00436            j = jbase
00437            do l=1,lot
00438            a0p2 = a(i0+i,l) + a(i2+i,l)
00439            a0m2 = a(i0+i,l) - a(i2+i,l)
00440            a1p3 = a(i1+i,l) + a(i3+i,l)
00441            a1m3 = a(i1+i,l) - a(i3+i,l)
00442            a4p6 = a(i4+i,l) + a(i6+i,l)
00443            a4m6 = a(i4+i,l) - a(i6+i,l)
00444            a5p7 = a(i5+i,l) + a(i7+i,l)
00445            a5m7 = a(i5+i,l) - a(i7+i,l)
00446
00447            a0p2m1p3 = a0p2 - a1p3
00448            a4m6m5m7 = a4m6 - a5m7
00449
00450            c(   j,l) = a0p2 + a1p3
00451            c(j4+j,l) = a4m6 + a5m7
00452            c(j2+j,l) = c2 * a0p2m1p3 - s2 * a4m6m5m7
00453            c(j6+j,l) = s2 * a0p2m1p3 + c2 * a4m6m5m7
00454            c(j1+j,l) = c1*(a0m2-a5p7)-s1*(a4p6+a1m3)
00455            c(j5+j,l) = s1*(a0m2-a5p7)+c1*(a4p6+a1m3)
00456            c(j3+j,l) = c3*(a0m2+a5p7)-s3*(a4p6-a1m3)
00457            c(j7+j,l) = s3*(a0m2+a5p7)+c3*(a4p6-a1m3)
00458            enddo
00459            jbase=jbase+1
00460         enddo
00461         i0 = i0 + iink
00462         i1 = i1 + iink
00463         i2 = i2 - iink
00464         i3 = i3 - iink
00465         i4 = i4 + iink
00466         i5 = i5 + iink
00467         i6 = i6 - iink
00468         i7 = i7 - iink
00469         jbase=jbase+7*la
00470 450     continue
00471
00472       if (i1 <= i2) then
00473          sin45=sqrt(0.5)
00474          do i=1,la
00475             j=jbase
00476             do l=1,lot
00477             c(   j,l)=a(i0+i,l)+a(i1+i,l)
00478             c(j1+j,l)=sin45*((a(i0+i,l)-a(i1+i,l))-(a(la+i0+i,l)+a(la+i1+i,l)))
00479             c(j2+j,l)=a(la+i1+i,l)-a(la+i0+i,l)
00480             c(j3+j,l)=-sin45*((a(i0+i,l)-a(i1+i,l))+(a(la+i0+i,l)+a(la+i1+i,l)))
00481             enddo
00482             jbase=jbase+1
00483          enddo
00484       endif
00485       la = la * 4
00486       return
00487       end subroutine ifft4
00488
00489 !     ================
00490 !     SUBROUTINE IFFT2
00491 !     ================
00492
00493       subroutine ifft2(a,trigs,n,lot,la)
00494       dimension a(n,lot),c(n,lot),trigs(n)
00495
00496       c(1,:) = 0.5 * a(1,:)
00497       c(2,:) = c(1,:)
00498
00499       ia    =   3
00500       ib    = n-1
00501
00502       do j = 3 , n-5 , 4
00503          c1 = trigs(ia  )
00504          s1 = trigs(ia+1)
00505          do l=1,lot
00506          amb = a(ia  ,l) - a(ib  ,l)
00507          apb = a(ia+1,l) + a(ib+1,l)
00508          c(j  ,l) = a(ia  ,l) + a(ib  ,l)
00509          c(j+2,l) = a(ia+1,l) - a(ib+1,l)
00510          c(j+1,l) = c1 * amb - s1 * apb
00511          c(j+3,l) = s1 * amb + c1 * apb
00512          enddo
00513          ia = ia + 2
00514          ib = ib - 2
00515       enddo
00516       c(n-1,:) =  a(ia  ,:)
00517       c(n  ,:) = -a(ia+1,:)
00518
00519       a(:,:) = c(:,:)
00520       la = 2
00521       return
00522       end subroutine ifft2
00523
00524 !     ================
00525 !     SUBROUTINE IFFT3
00526 !     ================
00527
00528       subroutine ifft3(a,trigs,n,lot,la)
00529       dimension a(n,lot),c(n,lot),trigs(n)
00530       parameter(SIN60 = 0.866025403784438D0)
00531
00532       ib = 2 * (n/3) + 1
00533
00534       c(1,:) = 0.5 * a(1,:) + a(ib,:)
00535       c(2,:) = 0.5 * a(1,:) - 0.5 * a(ib,:) - SIN60 * a(ib+1,:)
00536       c(3,:) = 0.5 * a(1,:) - 0.5 * a(ib,:) + SIN60 * a(ib+1,:)
00537
00538       ia = 3
00539       ic = ib - 2
00540       ib = ib + 2
00541
00542       do j = 4 , n-8 , 6
00543          c1 = trigs(ia  )
00544          s1 = trigs(ia+1)
00545          c2 = trigs(ia+ia-1)
00546          s2 = trigs(ia+ia  )
00547
00548          do l = 1 , lot
00549             hbpc = a(ia  ,l) - 0.5 * (a(ib  ,l) + a(ic  ,l))
00550             hbmc = a(ia+1,l) - 0.5 * (a(ib+1,l) - a(ic+1,l))
00551             sbmc = SIN60 * (a(ib  ,l) - a(ic  ,l))
00552             sbpc = SIN60 * (a(ib+1,l) + a(ic+1,l))
00553
00554             c(j  ,l) = a(ia  ,l) + a(ib  ,l) + a(ic  ,l)
00555             c(j+3,l) = a(ia+1,l) + a(ib+1,l) - a(ic+1,l)
00556             c(j+1,l) = c1 * (hbpc-sbpc) - s1 * (hbmc+sbmc)
00557             c(j+4,l) = s1 * (hbpc-sbpc) + c1 * (hbmc+sbmc)
00558             c(j+2,l) = c2 * (hbpc+sbpc) - s2 * (hbmc-sbmc)
00559             c(j+5,l) = s2 * (hbpc+sbpc) + c2 * (hbmc-sbmc)
00560          enddo
00561          ia = ia + 2
00562          ib = ib + 2
00563          ic = ic - 2
00564       enddo
00565
00566       c(n-2,:) = a(ia,:)
00567       c(n-1,:) =   0.5 * a(ia,:) - SIN60 * a(ia+1,:)
00568       c(n  ,:) = - 0.5 * a(ia,:) - SIN60 * a(ia+1,:)
00569
00570       a(:,:)  = c(:,:)
00571       la = 3
00572       return
00573       end subroutine ifft3
00574
00575 !     ================
00576 !     SUBROUTINE IFFT8
00577 !     ================
00578
00579       subroutine ifft8(a,c,n,lot)
00580       parameter(SQRT2 = 1.414213562373095D0)
00581       dimension a(n*lot),c(n*lot)
00582       la = n / 8
00583
00584       do i=0,la*lot-1
00585          i0 = (i/la) * n + mod(i,la) + 1
00586          i1 = i0 + la
00587          i2 = i1 + la
00588          i3 = i2 + la
00589          i4 = i3 + la
00590          i5 = i4 + la
00591          i6 = i5 + la
00592          i7 = i6 + la
00593
00594          a0p7 = a(i0) + a(i7)
00595          a0m7 = a(i0) - a(i7)
00596          a1p5 = a(i1) + a(i5)
00597          a1m5 = a(i1) - a(i5)
00598          a2p6 = a(i2) + a(i6)
00599          a2m6 = a(i2) - a(i6)
00600
00601          a0p7p3   = a0p7 + a(i3)
00602          a0p7m3   = a0p7 - a(i3)
00603          a0m7p4   = 2.0 * (a0m7 + a(i4))
00604          a0m7m4   = 2.0 * (a0m7 - a(i4))
00605          a1m5p2p6 = SQRT2 * (a1m5 + a2p6)
00606          a1m5m2p6 = SQRT2 * (a1m5 - a2p6)
00607
00608          c(i0)  = 2.0 * (a0p7p3 + a1p5)
00609          c(i2)  = 2.0 * (a0p7m3 - a2m6)
00610          c(i4)  = 2.0 * (a0p7p3 - a1p5)
00611          c(i6)  = 2.0 * (a0p7m3 + a2m6)
00612
00613          c(i1)  = a0m7m4 + a1m5m2p6
00614          c(i3)  = a0m7p4 - a1m5p2p6
00615          c(i5)  = a0m7m4 - a1m5m2p6
00616          c(i7)  = a0m7p4 + a1m5p2p6
00617       enddo
00618       return
00619       end