multpot.f90

基于linux操作下的fortran快速付立变换的程序

        subroutine multpot(icplexwf,icplex,includelast,nd1,nd2,n2,lot,n1dfft,pot,zw)
        implicit real*8 (a-h,o-z)
        dimension zw(2,lot,n2),pot(icplex*nd1,nd2)

        if(icplexwf==1)then

         if(icplex==2)then

          write(6,*)' multpot : icplexwf=1 and icplex=2, bug '
          stop
      
         else

!         TO BE SPEEDED UP : should use the same trick as Stefan
          if(includelast==1)then 
           do i2=1,n2
            do j=1,n1dfft
             zw(1,j,i2)=zw(1,j,i2)*pot(2*j-1,i2)
             zw(2,j,i2)=zw(2,j,i2)*pot(2*j  ,i2)
            enddo
           enddo
          else
           do i2=1,n2
            do j=1,n1dfft-1
             zw(1,j,i2)=zw(1,j,i2)*pot(2*j-1,i2)
             zw(2,j,i2)=zw(2,j,i2)*pot(2*j  ,i2)
            enddo
             zw(1,n1dfft,i2)=zw(1,n1dfft,i2)*pot(2*n1dfft-1,i2)
           enddo
          endif

         endif

        else if(icplexwf==2)then

         if(icplex==1)then

          do 10, i2=1,n2-1,2
          do 10, j=1,n1dfft
          zw(1,j,i2+0)=zw(1,j,i2+0)*pot(j,i2+0)
          zw(2,j,i2+0)=zw(2,j,i2+0)*pot(j,i2+0)
          zw(1,j,i2+1)=zw(1,j,i2+1)*pot(j,i2+1)
          zw(2,j,i2+1)=zw(2,j,i2+1)*pot(j,i2+1)
10        continue
          if(2*(n2/2)/=n2)then
          do 20, j=1,n1dfft
          zw(1,j,n2  )=zw(1,j,n2  )*pot(j,n2  )
          zw(2,j,n2  )=zw(2,j,n2  )*pot(j,n2  )
20        continue
          endif

         else

          do 30, i2=1,n2-1,2
          do 30, j=1,n1dfft
          zw(1,j,i2+0)=zw(1,j,i2+0)*pot(2*j-1,i2+0)&
&                     -zw(2,j,i2+0)*pot(2*j-0,i2+0)
          zw(2,j,i2+0)=zw(2,j,i2+0)*pot(2*j-1,i2+0)&
&                     +zw(1,j,i2+0)*pot(2*j-0,i2+0)
          zw(1,j,i2+1)=zw(1,j,i2+1)*pot(2*j-1,i2+1)&
&                     -zw(2,j,i2+1)*pot(2*j-0,i2+1)
          zw(2,j,i2+1)=zw(2,j,i2+1)*pot(2*j-1,i2+1)&
&                     +zw(1,j,i2+1)*pot(2*j-0,i2+1)
30        continue
          if(2*(n2/2)/=n2)then
          do 40, j=1,n1dfft
          zw(1,j,n2  )=zw(1,j,n2  )*pot(2*j-1,n2  )&
&                     -zw(2,j,n2  )*pot(2*j-0,n2  )
          zw(2,j,n2  )=zw(2,j,n2  )*pot(2*j-1,n2  )&
&                     +zw(1,j,n2  )*pot(2*j-0,n2  )
40        continue
          endif

         endif

        endif

        return
        end