Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer | :: | Ido | ||||
integer | :: | l1 | ||||
real(kind=rk) | :: | Cc | ||||
real(kind=rk) | :: | Ch | ||||
real(kind=rk) | :: | Wa1 | ||||
real(kind=rk) | :: | Wa2 | ||||
real(kind=rk) | :: | Wa3 | ||||
real(kind=rk) | :: | Wa4 |
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
real(kind=rk), | public | :: | ci2 | ||||
real(kind=rk), | public | :: | ci3 | ||||
real(kind=rk), | public | :: | ci4 | ||||
real(kind=rk), | public | :: | ci5 | ||||
real(kind=rk), | public | :: | cr2 | ||||
real(kind=rk), | public | :: | cr3 | ||||
real(kind=rk), | public | :: | cr4 | ||||
real(kind=rk), | public | :: | cr5 | ||||
real(kind=rk), | public | :: | di2 | ||||
real(kind=rk), | public | :: | di3 | ||||
real(kind=rk), | public | :: | di4 | ||||
real(kind=rk), | public | :: | di5 | ||||
real(kind=rk), | public | :: | dr2 | ||||
real(kind=rk), | public | :: | dr3 | ||||
real(kind=rk), | public | :: | dr4 | ||||
real(kind=rk), | public | :: | dr5 | ||||
integer, | public | :: | i | ||||
integer, | public | :: | ic | ||||
integer, | public | :: | idp2 | ||||
integer, | public | :: | k | ||||
real(kind=rk), | public, | parameter | :: | pi | = | acos(-1.0_rk) | |
real(kind=rk), | public, | parameter | :: | ti11 | = | sin(2.0_rk*pi/5.0_rk) | |
real(kind=rk), | public, | parameter | :: | ti12 | = | sin(4.0_rk*pi/5.0_rk) | |
real(kind=rk), | public | :: | ti2 | ||||
real(kind=rk), | public | :: | ti3 | ||||
real(kind=rk), | public | :: | ti4 | ||||
real(kind=rk), | public | :: | ti5 | ||||
real(kind=rk), | public, | parameter | :: | tr11 | = | cos(2.0_rk*pi/5.0_rk) | |
real(kind=rk), | public, | parameter | :: | tr12 | = | cos(4.0_rk*pi/5.0_rk) | |
real(kind=rk), | public | :: | tr2 | ||||
real(kind=rk), | public | :: | tr3 | ||||
real(kind=rk), | public | :: | tr4 | ||||
real(kind=rk), | public | :: | tr5 |
subroutine radb5(Ido,l1,Cc,Ch,Wa1,Wa2,Wa3,Wa4)
use fftpack_kind
implicit none
real(rk) :: Cc , Ch , ci2 , ci3 , ci4 , ci5 , cr2 , cr3 , &
cr4 , cr5 , di2 , di3 , di4 , di5 , dr2 , dr3 , &
dr4 , dr5
real(rk) :: ti2 , ti3 , ti4 , ti5 , tr2 , tr3, &
tr4 , tr5 , Wa1 , Wa2 , Wa3 , Wa4
integer :: i , ic , Ido , idp2 , k , l1
dimension Cc(Ido,5,l1) , Ch(Ido,l1,5) , Wa1(*) , Wa2(*) , Wa3(*), &
Wa4(*)
real(rk),parameter :: pi = acos(-1.0_rk)
real(rk),parameter :: tr11 = cos(2.0_rk * pi / 5.0_rk)
real(rk),parameter :: ti11 = sin(2.0_rk * pi / 5.0_rk)
real(rk),parameter :: tr12 = cos(4.0_rk * pi / 5.0_rk)
real(rk),parameter :: ti12 = sin(4.0_rk * pi / 5.0_rk)
do k = 1 , l1
ti5 = Cc(1,3,k) + Cc(1,3,k)
ti4 = Cc(1,5,k) + Cc(1,5,k)
tr2 = Cc(Ido,2,k) + Cc(Ido,2,k)
tr3 = Cc(Ido,4,k) + Cc(Ido,4,k)
Ch(1,k,1) = Cc(1,1,k) + tr2 + tr3
cr2 = Cc(1,1,k) + tr11*tr2 + tr12*tr3
cr3 = Cc(1,1,k) + tr12*tr2 + tr11*tr3
ci5 = ti11*ti5 + ti12*ti4
ci4 = ti12*ti5 - ti11*ti4
Ch(1,k,2) = cr2 - ci5
Ch(1,k,3) = cr3 - ci4
Ch(1,k,4) = cr3 + ci4
Ch(1,k,5) = cr2 + ci5
enddo
if ( Ido==1 ) return
idp2 = Ido + 2
do k = 1 , l1
do i = 3 , Ido , 2
ic = idp2 - i
ti5 = Cc(i,3,k) + Cc(ic,2,k)
ti2 = Cc(i,3,k) - Cc(ic,2,k)
ti4 = Cc(i,5,k) + Cc(ic,4,k)
ti3 = Cc(i,5,k) - Cc(ic,4,k)
tr5 = Cc(i-1,3,k) - Cc(ic-1,2,k)
tr2 = Cc(i-1,3,k) + Cc(ic-1,2,k)
tr4 = Cc(i-1,5,k) - Cc(ic-1,4,k)
tr3 = Cc(i-1,5,k) + Cc(ic-1,4,k)
Ch(i-1,k,1) = Cc(i-1,1,k) + tr2 + tr3
Ch(i,k,1) = Cc(i,1,k) + ti2 + ti3
cr2 = Cc(i-1,1,k) + tr11*tr2 + tr12*tr3
ci2 = Cc(i,1,k) + tr11*ti2 + tr12*ti3
cr3 = Cc(i-1,1,k) + tr12*tr2 + tr11*tr3
ci3 = Cc(i,1,k) + tr12*ti2 + tr11*ti3
cr5 = ti11*tr5 + ti12*tr4
ci5 = ti11*ti5 + ti12*ti4
cr4 = ti12*tr5 - ti11*tr4
ci4 = ti12*ti5 - ti11*ti4
dr3 = cr3 - ci4
dr4 = cr3 + ci4
di3 = ci3 + cr4
di4 = ci3 - cr4
dr5 = cr2 + ci5
dr2 = cr2 - ci5
di5 = ci2 - cr5
di2 = ci2 + cr5
Ch(i-1,k,2) = Wa1(i-2)*dr2 - Wa1(i-1)*di2
Ch(i,k,2) = Wa1(i-2)*di2 + Wa1(i-1)*dr2
Ch(i-1,k,3) = Wa2(i-2)*dr3 - Wa2(i-1)*di3
Ch(i,k,3) = Wa2(i-2)*di3 + Wa2(i-1)*dr3
Ch(i-1,k,4) = Wa3(i-2)*dr4 - Wa3(i-1)*di4
Ch(i,k,4) = Wa3(i-2)*di4 + Wa3(i-1)*dr4
Ch(i-1,k,5) = Wa4(i-2)*dr5 - Wa4(i-1)*di5
Ch(i,k,5) = Wa4(i-2)*di5 + Wa4(i-1)*dr5
enddo
enddo
end subroutine radb5