! File: scaled2d.f90 ! Public domain 2004 James Van Buskirk module mykinds implicit none integer, parameter :: dp = selected_real_kind(15,300) integer, parameter :: qp_preferred = selected_real_kind(30,1000) integer, parameter :: qp = (1+sign(1,qp_preferred))/2*qp_preferred+ & (1-sign(1,qp_preferred))/2*dp end module mykinds module avl_mod use mykinds implicit none type avl_node real(qp) data ! Quadruple precision value represented integer index ! Index in list of constants type(avl_node), pointer :: left => NULL() type(avl_node), pointer :: right => NULL() integer :: balance = 0 ! -1 = heavy left; +1 = heavy right end type avl_node interface operator(<) module procedure lessn end interface operator(<) interface operator(>) module procedure greatern end interface operator(>) integer, parameter :: ulps = 100 private ulps type absorber type(avl_node), pointer :: root => NULL() integer :: num_nodes = 0 end type absorber contains function lessn(a, b) logical lessn real(qp), intent(in) :: a type(avl_node), intent(in) :: b lessn = abs(a) < abs(b%data)-ulps*spacing(b%data) end function lessn function greatern(a, b) logical greatern real(qp), intent(in) :: a type(avl_node), intent(in) :: b greatern = abs(a) > abs(b%data)+ulps*spacing(b%data) end function greatern subroutine digest(food,body,poop) real(qp), intent(in) :: food type(absorber) body integer, intent(out) :: poop integer balance if(.NOT.associated(body%root)) then allocate(body%root) body%num_nodes = body%num_nodes+1 poop = body%num_nodes body%root%data = abs(food) body%root%index = poop else call avl_insert(food,body,body%root,poop,balance) end if end subroutine digest recursive subroutine avl_insert(x,body,root,index,balance) real(qp), intent(in) :: x type(absorber) body type(avl_node), pointer :: root integer, intent(out) :: index integer, intent(out) :: balance type(avl_node), pointer :: rotate_temp if(x < root) then if(.NOT.associated(root%left)) then allocate(root%left) body%num_nodes = body%num_nodes+1 index = body%num_nodes root%left%data = abs(x) root%left%index = index balance = 1-root%balance root%balance = root%balance-1 else call avl_insert(x,body,root%left,index,balance) if(balance == 1) then if(root%balance == -1) then if(root%left%balance == 1) then ! double rotate call force_left(root%left) end if rotate_temp => root%left root%left => rotate_temp%right rotate_temp%right => root root => rotate_temp root%balance = 0 root%right%balance = 0 balance = 0 else if(root%balance == 0) then root%balance = -1 ! balance = 1 else ! root%balance == 1 root%balance = 0 balance = 0 end if end if end if else if(x > root) then if(.NOT.associated(root%right)) then allocate(root%right) body%num_nodes = body%num_nodes+1 index = body%num_nodes root%right%data = abs(x) root%right%index = index balance = 1+root%balance root%balance = root%balance+1 else call avl_insert(x,body,root%right,index,balance) if(balance == 1) then if(root%balance == 1) then if(root%right%balance == -1) then ! double rotate call force_right(root%right) end if rotate_temp => root%right root%right => rotate_temp%left rotate_temp%left => root root => rotate_temp root%balance = 0 root%left%balance = 0 balance = 0 else if(root%balance == 0) then root%balance = 1 ! balance = 1 else ! root%balance == -1 root%balance = 0 balance = 0 end if end if end if else ! x == root index = root%index balance = 0 end if end subroutine avl_insert recursive subroutine force_left(root) type(avl_node), pointer :: root type(avl_node), pointer :: rotate_temp if(root%right%balance == -1) then call force_right(root%right) end if rotate_temp => root%right root%right => rotate_temp%left rotate_temp%left => root root => rotate_temp root%left%balance = -root%balance root%balance = -1 end subroutine force_left recursive subroutine force_right(root) type(avl_node), pointer :: root type(avl_node), pointer :: rotate_temp if(root%left%balance == 1) then call force_left(root%left) end if rotate_temp => root%left root%left => rotate_temp%right rotate_temp%right => root root => rotate_temp root%right%balance = -root%balance root%balance = 1 end subroutine force_right subroutine list2array(body,array) type(absorber) body real(qp), intent(out) :: array(body%num_nodes) if(associated(body%root)) then call clean(body%root,array) end if end subroutine list2array recursive subroutine clean(root,array) type(avl_node), pointer :: root real(qp) array(*) if(associated(root%left)) then call clean(root%left,array) end if array(root%index) = root%data if(associated(root%right)) then call clean(root%right,array) end if deallocate(root) end subroutine clean end module avl_mod module correct use mykinds use avl_mod implicit none ! type(absorber) body type(absorber), save :: body type get_factor_stuff real(qp), pointer :: factor(:) end type get_factor_stuff type(get_factor_stuff), allocatable :: get_factor_data(:) type odd_tail_stuff integer, pointer :: matrixc(:,:) integer, pointer :: matrixs(:,:) end type odd_tail_stuff type(odd_tail_stuff), allocatable :: odd_tail_data(:) type scaled_odd_tail_stuff integer, pointer :: matrixc(:) integer, pointer :: matrixs(:) integer, pointer :: matrix1c(:,:) integer, pointer :: matrix1s(:,:) integer, pointer :: diag(:) end type scaled_odd_tail_stuff type(scaled_odd_tail_stuff), allocatable :: scaled_odd_tail_data(:) type correct_stuff integer rsq2 integer c1 integer c3 end type correct_stuff type(correct_stuff) correct_data real(dp), allocatable :: const_array(:) type var_stuff integer depth integer index character(2) label end type var_stuff contains ! Subroutine get_factor ! Computes scale_factor(n) recursive subroutine get_factor(factor, n, old_factor) integer, intent(in) :: n real(qp), intent(out) :: factor(2**(n-4)) real(qp), intent(in) :: old_factor(2**(n-6)) integer i integer k real(qp) pi pi = 4*atan(1.0_qp) if(n <= 3) then factor = 1 else factor = (/(cos(2*pi*i/2**n),i=1,2**(n-3),2)/) if(n > 5) then k = 2**(n-6) factor(0*k+1:0*k+k) = factor(0*k+1:0*k+k)*old_factor factor(1*k+1:1*k+k) = factor(1*k+1:1*k+k)*old_factor(k:1:-1) factor(2*k+1:2*k+k) = factor(2*k+1:2*k+k)*old_factor factor(3*k+1:3*k+k) = factor(3*k+1:3*k+k)*old_factor(k:1:-1) end if end if end subroutine get_factor ! Subroutine prepare_correct(n) subroutine prepare_correct(n) ! type(absorber) body integer, intent(in) :: n real(qp) constant real(qp) pi if(n > 2) then pi = 4*atan(1.0_qp) constant = cos(pi/4) call digest(constant, body, correct_data%rsq2) end if if(n > 3) then constant = cos(pi/8) call digest(constant, body, correct_data%c1) constant = cos(3*pi/8) call digest(constant, body, correct_data%c3) end if end subroutine prepare_correct ! Subroutine scaled_odd_tail_correct ! Calculates scaled_odd_tail(n)*h recursive subroutine scaled_odd_tail_correct(h,n) integer, intent(in) :: n ! real(dp), intent(inout) :: h(2**(n-1)) type(var_stuff), intent(inout) :: h(2**(n-1)) ! real(dp), parameter :: rsq2 = 0.70710678118654752440084436210485_dp ! real(dp) xc ! real(dp) xs type(var_stuff) h1(2**(n-1)), h2(2**(n-1)) integer max_depth integer i select case(n) case(1) case(2) case(3) max_depth = maxval(h%depth) h1(1) = var_stuff(max_depth+1, h(1)%index, h(1)%label) write(10,'(3a,i0,5a)') ' ', write_var(h1(1)),' = d',correct_data%rsq2, & '*(',write_var(h(2)),'-',write_var(h(4)),')' ! xc = const_array(correct_data%rsq2)*(h(2)-h(4)) h1(2) = var_stuff(max_depth+1, h(2)%index, h(2)%label) write(10,'(3a,i0,5a)') ' ', write_var(h1(2)),' = d',correct_data%rsq2, & '*(',write_var(h(2)),'+',write_var(h(4)),')' ! xs = const_array(correct_data%rsq2)*(h(2)+h(4)) do i = 1, 2**(n-1) h2(i) = var_stuff(max_depth+2,h(i)%index,h(i)%label) end do write(10,'(6a)') ' ',write_var(h2(1)),' = ', & write_var(h(1)),'+',write_var(h1(1)) write(10,'(6a)') ' ',write_var(h2(2)),' = ', & write_var(h(1)),'-',write_var(h1(1)) write(10,'(6a)') ' ',write_var(h2(3)),' = ', & write_var(h1(2)),'+',write_var(h(3)) write(10,'(6a)') ' ',write_var(h2(4)),' = ', & write_var(h1(2)),'-',write_var(h(3)) ! h = (/h(1)+xc,h(1)-xc,xs+h(3),xs-h(3)/) do i = 1, 2**(n-1) h(i) = var_stuff(max_depth+2,h2(i)%index,h2(i)%label) end do case default stop ' bigger n than 3 in scaled_odd_tail_correct' end select end subroutine scaled_odd_tail_correct ! Subroutine odd_tail_correct ! Computes odd_tail(n)*h using the theoretical formula subroutine odd_tail_correct(h,n) integer, intent(in) :: n ! real(dp), intent(inout) :: h(2**(n-1)) type(var_stuff), intent(inout) :: h(2**(n-1)) ! real(dp) xc ! real(dp) xs ! real(dp) x3_5, x1_7, x35, x17, y3, y4, y5, y6 ! real(dp) x0, x2_6, x4, x26, y1, y2, y7, y8 type(var_stuff) h1(2**(n-1)), h2(2**(n-1)) ! real(dp), parameter :: rsq2 = 0.70710678118654752440084436210485_dp ! real(dp), parameter :: c1 = 0.92387953251128675612818318939679_dp ! real(dp), parameter :: c3 = 0.3826834323650897717284599840304_dp integer max_depth integer i select case(n) case(1) case(2) case(3) max_depth = maxval(h%depth) h1(1) = var_stuff(max_depth+1, h(1)%index, h(1)%label) write(10,'(3a,i0,5a)') ' ', write_var(h1(1)),' = d',correct_data%rsq2, & '*(',write_var(h(2)),'-',write_var(h(4)),')' ! xc = const_array(correct_data%rsq2)*(h(2)-h(4)) h1(2) = var_stuff(max_depth+1, h(2)%index, h(2)%label) write(10,'(3a,i0,5a)') ' ', write_var(h1(2)),' = d',correct_data%rsq2, & '*(',write_var(h(2)),'+',write_var(h(4)),')' ! xs = const_array(correct_data%rsq2)*(h(2)+h(4)) do i = 1, 2**(n-1) h2(i) = var_stuff(max_depth+2,h(i)%index,h(i)%label) end do write(10,'(6a)') ' ',write_var(h2(1)),' = ', & write_var(h(1)),'+',write_var(h1(1)) write(10,'(6a)') ' ',write_var(h2(2)),' = ', & write_var(h(1)),'-',write_var(h1(1)) write(10,'(6a)') ' ',write_var(h2(3)),' = ', & write_var(h1(2)),'+',write_var(h(3)) write(10,'(6a)') ' ',write_var(h2(4)),' = ', & write_var(h1(2)),'-',write_var(h(3)) ! h = (/h(1)+xc,h(1)-xc,xs+h(3),xs-h(3)/) do i = 1, 2**(n-1) h(i) = var_stuff(max_depth+2,h2(i)%index,h2(i)%label) end do case(4) max_depth = maxval(h%depth) do i = 1, 2**(n-1) h1(i) = var_stuff(max_depth+1,h(i)%index,h(i)%label) h2(i) = var_stuff(max_depth+2,h(i)%index,h(i)%label) end do write(10,'(6a)') ' ',write_var(h1(2)),' = ', & write_var(h(2)),'+',write_var(h(8)) ! x17 = h(2)+h(8) write(10,'(6a)') ' ',write_var(h1(6)),' = ', & write_var(h(2)),'-',write_var(h(8)) ! x1_7 = h(2)-h(8) write(10,'(6a)') ' ',write_var(h1(4)),' = ', & write_var(h(4)),'+',write_var(h(6)) ! x35 = h(4)+h(6) write(10,'(6a)') ' ',write_var(h1(8)),' = ', & write_var(h(4)),'-',write_var(h(6)) ! x3_5 = h(4)-h(6) write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(h2(3)),' = d',& correct_data%c1,'*',write_var(h1(6)),'+d', & correct_data%c3,'*',write_var(h1(8)) ! y3 = correct_data%c1*x1_7+correct_data%c3*x3_5 write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(h2(4)),' = d',& correct_data%c3,'*',write_var(h1(6)),'-d', & correct_data%c1,'*',write_var(h1(8)) ! y4 = correct_data%c3*x1_7-correct_data%c1*x3_5 write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(h2(5)),' = d',& correct_data%c3,'*',write_var(h1(2)),'+d', & correct_data%c1,'*',write_var(h1(4)) ! y5 = correct_data%c3*x17+correct_data%c1*x35 write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(h2(6)),' = d',& correct_data%c1,'*',write_var(h1(2)),'-d', & correct_data%c3,'*',write_var(h1(4)) ! y6 = correct_data%c1*x17-correct_data%c3*x35 ! x0 = h(1) write(10,'(3a,i0,5a)') ' ',write_var(h1(7)),' = d', & correct_data%rsq2,'*(',write_var(h(3)), & '-', write_var(h(7)), ')' ! x2_6 = correct_data%rsq2*(h(3)-h(7)) ! x4 = h(5) write(10,'(3a,i0,5a)') ' ',write_var(h1(3)),' = d', & correct_data%rsq2,'*(',write_var(h(3)), & '+', write_var(h(7)), ')' ! x26 = correct_data%rsq2*(h(3)+h(7)) write(10,'(6a)') ' ',write_var(h2(1)),' = ', & write_var(h(1)),'+',write_var(h1(7)) ! y1 = x0+x2_6 write(10,'(6a)') ' ',write_var(h2(2)),' = ', & write_var(h(1)),'-',write_var(h1(7)) ! y2 = x0-x2_6 write(10,'(6a)') ' ',write_var(h2(7)),' = ', & write_var(h1(3)),'+',write_var(h(5)) ! y7 = x26+x4 write(10,'(6a)') ' ',write_var(h2(8)),' = ', & write_var(h1(3)),'-',write_var(h(5)) ! y8 = x26-x4 do i = 1, 2**(n-1) h(i) = var_stuff(max_depth+3, h(i)%index, h(i)%label) end do write(10,'(6a)') ' ',write_var(h(1)),' = ',& write_var(h2(1)),'+',write_var(h2(3)) write(10,'(6a)') ' ',write_var(h(2)),' = ',& write_var(h2(2)),'+',write_var(h2(4)) write(10,'(6a)') ' ',write_var(h(3)),' = ',& write_var(h2(2)),'-',write_var(h2(4)) write(10,'(6a)') ' ',write_var(h(4)),' = ',& write_var(h2(1)),'-',write_var(h2(3)) write(10,'(6a)') ' ',write_var(h(5)),' = ',& write_var(h2(5)),'+',write_var(h2(7)) write(10,'(6a)') ' ',write_var(h(6)),' = ',& write_var(h2(6)),'+',write_var(h2(8)) write(10,'(6a)') ' ',write_var(h(7)),' = ',& write_var(h2(6)),'-',write_var(h2(8)) write(10,'(6a)') ' ',write_var(h(8)),' = ',& write_var(h2(5)),'-',write_var(h2(7)) ! h = (/y1+y3,y2+y4,y2-y4,y1-y3,y5+y7,y6+y8,y6-y8,y5-y7/) case default stop ' bigger n than 4 in odd_tail_correct' end select end subroutine odd_tail_correct ! Subroutine odd_tail_gen(h,n) ! Compute odd_tail constants and opcounts subroutine odd_tail_gen(stuff,n) integer, intent(in) :: n type(odd_tail_stuff), intent(out) :: stuff integer i integer j integer k real(qp) pi integer :: ulps = 100 real(qp) matrix(2**(n-3),2**(n-3)) pi = 4*atan(1.0_qp) allocate(stuff%matrixc(2,2**(n-3))) allocate(stuff%matrixs(2,2**(n-3))) if(n > 4) then k = 2**(n-4) do i = 1, k matrix(1,i) = cos(2*pi*(2*i-1)/2**n) matrix(1,2*k+1-i) = sin(2*pi*(2*i-1)/2**n) matrix(2,i) = -sin(2*pi*(2*i-1)/2**n) matrix(2,2*k+1-i) = cos(2*pi*(2*i-1)/2**n) end do k = 2**(n-6) do j = 1, k matrix(:,0*k+j) = matrix(:,0*k+j)* & get_factor_data(n-2)%factor(j) matrix(:,8*k+1-j) = matrix(:,8*k+1-j)* & get_factor_data(n-2)%factor(j) matrix(:,1*k+j) = matrix(:,1*k+j)* & get_factor_data(n-2)%factor(k+1-j) matrix(:,7*k+1-j) = matrix(:,7*k+1-j)* & get_factor_data(n-2)%factor(k+1-j) matrix(:,2*k+j) = matrix(:,2*k+j)* & get_factor_data(n-2)%factor(j) matrix(:,6*k+1-j) = matrix(:,6*k+1-j)* & get_factor_data(n-2)%factor(j) matrix(:,3*k+j) = matrix(:,3*k+j)* & get_factor_data(n-2)%factor(k+1-j) matrix(:,5*k+1-j) = matrix(:,5*k+1-j)* & get_factor_data(n-2)%factor(k+1-j) end do do i = 1, 2**(n-3) do j = 1, 2 call digest(abs(matrix(j,i)), body, stuff%matrixc(j,i)) end do end do k = 2**(n-4) do i = 1, k matrix(1,i) = sin(2*pi*(2*i-1)/2**n) matrix(1,2*k+1-i) = cos(2*pi*(2*i-1)/2**n) matrix(2,i) = cos(2*pi*(2*i-1)/2**n) matrix(2,2*k+1-i) = -sin(2*pi*(2*i-1)/2**n) end do k = 2**(n-6) do j = 1, k matrix(:,0*k+j) = matrix(:,0*k+j)* & get_factor_data(n-2)%factor(j) matrix(:,8*k+1-j) = matrix(:,8*k+1-j)* & get_factor_data(n-2)%factor(j) matrix(:,1*k+j) = matrix(:,1*k+j)* & get_factor_data(n-2)%factor(k+1-j) matrix(:,7*k+1-j) = matrix(:,7*k+1-j)* & get_factor_data(n-2)%factor(k+1-j) matrix(:,2*k+j) = matrix(:,2*k+j)* & get_factor_data(n-2)%factor(j) matrix(:,6*k+1-j) = matrix(:,6*k+1-j)* & get_factor_data(n-2)%factor(j) matrix(:,3*k+j) = matrix(:,3*k+j)* & get_factor_data(n-2)%factor(k+1-j) matrix(:,5*k+1-j) = matrix(:,5*k+1-j)* & get_factor_data(n-2)%factor(k+1-j) end do do i = 1, 2**(n-3) do j = 1, 2 call digest(abs(matrix(j,i)), body, stuff%matrixs(j,i)) end do end do end if end subroutine odd_tail_gen ! Subroutine odd_tail(h,n) ! Compute action of odd tail matrix of order 2**(n-1) recursive subroutine odd_tail(h,stuff,n) integer, intent(in) :: n type(odd_tail_stuff), intent(in) :: stuff type(var_stuff), intent(inout) :: h(2**(n-1)) type(var_stuff) hc(2**(n-3)), hc1(2**(n-3)) type(var_stuff) hs(2**(n-3)), hs1(2**(n-3)) type(var_stuff) h2(2**(n-2)) integer i integer j integer k integer M integer max_depth if(n <= 4) then call odd_tail_correct(h,n) else max_depth = maxval(h%depth) M = 2**(n-3) do i = 1, M h2(i) = h(2*i-1) hc(i) = var_stuff(max_depth+1,h(2*i)%index,h(2*i)%label) hs(i) = var_stuff(max_depth+1,h(2*(i+M))%index,h(2*(i+M))%label) h2(i+M) = h(2*(i+M)-1) end do do i = 1, M write(10,'(6a)') ' ',write_var(hc(i)),' = ', & write_var(h(4*i-2)),'-',write_var(h(4*(M+1-i))) end do ! hc = h(2::4)-h(2**(n-1):1:-4) call scaled_odd_tail(hc,scaled_odd_tail_data(n-2),n-2) max_depth = maxval(hc%depth) do i = 1, M hc1(i) = var_stuff(max_depth+1,hc(i)%index,hc(i)%label) end do M = 2**(n-4) do i = 1, M write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(hc1(i)),' = d', & stuff%matrixc(1,i),'*',write_var(hc(i)),'-d', & stuff%matrixc(2,i),'*',write_var(hc(M+i)) end do do i = M+1,2*M write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(hc1(i)),' = d', & stuff%matrixc(1,i),'*',write_var(hc(2*M+1-i)),'+d', & stuff%matrixc(2,i),'*',write_var(hc(3*M+1-i)) end do ! hc = (/(/(stuff%matrixc(1,i)*hc(i)- & ! stuff%matrixc(2,i)*hc(M+i),i=1,M)/), & ! (/(stuff%matrixc(1,i)*hc(2*M+1-i)+ & ! stuff%matrixc(2,i)*hc(3*M+1-i),i=M+1,2*M)/)/) M = 2**(n-3) do i = 1, M write(10,'(6a)') ' ',write_var(hs(i)),' = ', & write_var(h(4*i-2)),'+',write_var(h(4*(M+1-i))) end do ! hs = h(2::4)+h(2**(n-1):1:-4) call scaled_odd_tail(hs,scaled_odd_tail_data(n-2),n-2) max_depth = maxval(hs%depth) do i = 1, M hs1(i) = var_stuff(max_depth+1,hs(i)%index,hs(i)%label) end do M = 2**(n-4) do i = 1, M write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(hs1(i)),' = d', & stuff%matrixs(1,i),'*',write_var(hs(i)),'+d', & stuff%matrixs(2,i),'*',write_var(hs(M+i)) end do do i = M+1,2*M write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(hs1(i)),' = d', & stuff%matrixs(1,i),'*',write_var(hs(2*M+1-i)),'-d', & stuff%matrixs(2,i),'*',write_var(hs(3*M+1-i)) end do ! hs = (/(/(stuff%matrixs(1,i)*hs(i)+ & ! stuff%matrixs(2,i)*hs(M+i),i=1,M)/), & ! (/(stuff%matrixs(1,i)*hs(2*M+1-i)- & ! stuff%matrixs(2,i)*hs(3*M+1-i),i=M+1,2*M)/)/) ! h2 = h(::2) call odd_tail(h2,odd_tail_data(n-1),n-1) k = 2**(n-3) max_depth = max(maxval(hc1%depth),maxval(hs1%depth),maxval(h2%depth)) do i = 1, k h(2*i-1) = var_stuff(max_depth+1,h2(i)%index,h2(i)%label) h(2*i) = var_stuff(max_depth+1,hc1(i)%index,hc1(i)%label) h(2*(i+k)) = var_stuff(max_depth+1,hs1(i)%index,hs1(i)%label) h(2*(i+k)-1) = var_stuff(max_depth+1,h2(i+k)%index,h2(i+k)%label) end do do i = 1, k write(10,'(6a)') ' ',write_var(h(i)),' = ', & write_var(h2(i)),'+',write_var(hc1(i)) end do ! h(0*k+1:1*k) = h2(:k)+hc do i = 1, k write(10,'(6a)') ' ',write_var(h(k+i)),' = ', & write_var(h2(k+1-i)),'-',write_var(hc1(k+1-i)) end do ! h(1*k+1:2*k) = h2(k:1:-1)-hc(k:1:-1) do i = 1, k write(10,'(6a)') ' ',write_var(h(2*k+i)),' = ', & write_var(hs1(i)),'+',write_var(h2(k+i)) end do ! h(2*k+1:3*k) = hs+h2(k+1:) do i = 1, k write(10,'(6a)') ' ',write_var(h(3*k+i)),' = ', & write_var(hs1(k+1-i)),'-',write_var(h2(2*k+1-i)) end do ! h(3*k+1:4*k) = hs(k:1:-1)-h2(2*k:k+1:-1) end if end subroutine odd_tail ! Subroutine scaled_odd_tail_gen ! Generates constants and opcounts for scaled_odd_tail subroutine scaled_odd_tail_gen(stuff,n) integer, intent(in) :: n type(scaled_odd_tail_stuff), intent(out) :: stuff real(qp) factor(2**(n-4)) integer i integer j integer k real(qp) pi integer :: ulps = 100 real(qp) matrix(2**(n-3)) real(qp) matrix1(min(2,2**(n-4)),2**(n-4)) real(qp) diag(2**(n-3)) pi = 4*atan(1.0_qp) allocate(stuff%matrixc(2**(n-3))) allocate(stuff%matrixs(2**(n-3))) allocate(stuff%matrix1c(min(2,2**(n-4)),2**(n-4))) allocate(stuff%matrix1s(min(2,2**(n-4)),2**(n-4))) allocate(stuff%diag(2**(n-3))) if(n == 3) then ! Insert garbage values: should never be used stuff%matrixc = 0 stuff%matrixs = 0 stuff%diag = 0 else if(n > 3) then factor = get_factor_data(n)%factor ! matrix = 0 k = 2**(n-4) do i = 1, k matrix(2*k+1-i) = tan(2*pi*(2*i-1)/2**n) matrix(i) = -tan(2*pi*(2*i-1)/2**n) end do do i = 1, 2**(n-3) call digest(abs(matrix(i)), body, stuff%matrixc(i)) end do ! ! matrix = 0 k = 2**(n-4) do i = 1, k matrix(i) = tan(2*pi*(2*i-1)/2**n) matrix(2*k+1-i) = -tan(2*pi*(2*i-1)/2**n) end do do i = 1, 2**(n-3) call digest(abs(matrix(i)), body, stuff%matrixs(i)) end do ! ! matrix1 = 0 k = 2**(n-5) matrix1(1,1) = cos(pi/4) do i = 1, k matrix1(1,i) = cos(2*pi*(2*i-1)/2**(n-1)) matrix1(1,2*k+1-i) = sin(2*pi*(2*i-1)/2**(n-1)) matrix1(2,i) = -sin(2*pi*(2*i-1)/2**(n-1)) matrix1(2,2*k+1-i) = cos(2*pi*(2*i-1)/2**(n-1)) end do k = 2**(n-7) do j = 1, k matrix1(:,0*k+j) = matrix1(:,0*k+j)* & get_factor_data(n-3)%factor(j) matrix1(:,8*k+1-j) = matrix1(:,8*k+1-j)* & get_factor_data(n-3)%factor(j) matrix1(:,1*k+j) = matrix1(:,1*k+j)* & get_factor_data(n-3)%factor(k+1-j) matrix1(:,7*k+1-j) = matrix1(:,7*k+1-j)* & get_factor_data(n-3)%factor(k+1-j) matrix1(:,2*k+j) = matrix1(:,2*k+j)* & get_factor_data(n-3)%factor(j) matrix1(:,6*k+1-j) = matrix1(:,6*k+1-j)* & get_factor_data(n-3)%factor(j) matrix1(:,3*k+j) = matrix1(:,3*k+j)* & get_factor_data(n-3)%factor(k+1-j) matrix1(:,5*k+1-j) = matrix1(:,5*k+1-j)* & get_factor_data(n-3)%factor(k+1-j) end do k = 2**(n-4) do j = 1, k matrix1(:,0*k+j) = matrix1(:,0*k+j)/factor(j) end do do i = 1, 2**(n-4) do j = 1, min(2,2**(n-4)) call digest(abs(matrix1(j,i)), body, stuff%matrix1c(j,i)) end do end do ! matrix1 = 0 k = 2**(n-5) matrix1(1,1) = cos(pi/4) do i = 1, k matrix1(1,i) = sin(2*pi*(2*i-1)/2**(n-1)) matrix1(1,2*k+1-i) = cos(2*pi*(2*i-1)/2**(n-1)) matrix1(2,i) = cos(2*pi*(2*i-1)/2**(n-1)) matrix1(2,2*k+1-i) = -sin(2*pi*(2*i-1)/2**(n-1)) end do k = 2**(n-7) do j = 1, k matrix1(:,0*k+j) = matrix1(:,0*k+j)* & get_factor_data(n-3)%factor(j) matrix1(:,8*k+1-j) = matrix1(:,8*k+1-j)* & get_factor_data(n-3)%factor(j) matrix1(:,1*k+j) = matrix1(:,1*k+j)* & get_factor_data(n-3)%factor(k+1-j) matrix1(:,7*k+1-j) = matrix1(:,7*k+1-j)* & get_factor_data(n-3)%factor(k+1-j) matrix1(:,2*k+j) = matrix1(:,2*k+j)* & get_factor_data(n-3)%factor(j) matrix1(:,6*k+1-j) = matrix1(:,6*k+1-j)* & get_factor_data(n-3)%factor(j) matrix1(:,3*k+j) = matrix1(:,3*k+j)* & get_factor_data(n-3)%factor(k+1-j) matrix1(:,5*k+1-j) = matrix1(:,5*k+1-j)* & get_factor_data(n-3)%factor(k+1-j) end do k = 2**(n-4) do j = 1, k matrix1(:,0*k+j) = matrix1(:,0*k+j)/factor(j) end do do i = 1, 2**(n-4) do j = 1, min(2,2**(n-4)) call digest(abs(matrix1(j,i)), body, stuff%matrix1s(j,i)) end do end do k = 2**(n-6) diag = 1 do j = 1, k diag(0*k+j) = get_factor_data(n-2)%factor(j) diag(1*k+j) = get_factor_data(n-2)%factor(k+1-j) diag(2*k+j) = get_factor_data(n-2)%factor(j) diag(3*k+j) = get_factor_data(n-2)%factor(k+1-j) diag(4*k+j) = get_factor_data(n-2)%factor(j) diag(5*k+j) = get_factor_data(n-2)%factor(k+1-j) diag(6*k+j) = get_factor_data(n-2)%factor(j) diag(7*k+j) = get_factor_data(n-2)%factor(k+1-j) end do k = 2**(n-4) do j = 1, k diag(0*k+j) = diag(0*k+j)/factor(j) diag(1*k+j) = diag(1*k+j)/factor(j) end do do j = 1, 2**(n-3) call digest(abs(diag(j)), body, stuff%diag(j)) end do end if end subroutine scaled_odd_tail_gen ! Subroutine scaled_odd_tail ! Computes action of scaled odd tail matrix of order 2**(n-1) recursive subroutine scaled_odd_tail(h,stuff,n) integer, intent(in) :: n type(var_stuff), intent(inout) :: h(2**(n-1)) type(scaled_odd_tail_stuff), intent(in) :: stuff type(var_stuff) hc(2**(n-3)), hcx(2**(n-3)) type(var_stuff) hs(2**(n-3)), hsx(2**(n-3)) type(var_stuff) h2(2**(n-2)) type(var_stuff) hc1(2**(n-4)), hc1x(2**(n-4)) type(var_stuff) hs1(2**(n-4)), hs1x(2**(n-4)) type(var_stuff) h21(2**(n-3)), h21x(2**(n-3)) integer i integer j integer k integer M integer max_depth if(n <= 3) then call odd_tail_correct(h,n) else max_depth = maxval(h%depth) M = 2**(n-3) do i = 1, M h2(i) = h(2*i-1) hc(i) = var_stuff(max_depth+1,h(2*i)%index,h(2*i)%label) hs(i) = var_stuff(max_depth+1,h(2*(i+M))%index,h(2*(i+M))%label) h2(i+M) = h(2*(i+M)-1) end do do i = 1, M write(10,'(6a)') ' ',write_var(hc(i)),' = ', & write_var(h(4*i-2)),'-',write_var(h(4*(M+1-i))) end do ! hc = h(2::4)-h(2**(n-1):1:-4) call scaled_odd_tail(hc,scaled_odd_tail_data(n-2),n-2) max_depth = maxval(hc%depth) do i = 1, M hcx(i) = var_stuff(max_depth+1,hc(i)%index,hc(i)%label) end do M = 2**(n-4) do i = 1, M write(10,'(5a,i0,2a)') ' ',write_var(hcx(i)),' = ', & write_var(hc(i)),'-d', & stuff%matrixc(i),'*',write_var(hc(M+i)) end do do i = M+1,2*M write(10,'(3a,i0,4a)') ' ',write_var(hcx(i)),' = d', & stuff%matrixc(i),'*',write_var(hc(2*M+1-i)),'+', & write_var(hc(3*M+1-i)) end do ! hc = (/(/(hc(i)-stuff%matrixc(i)*hc(M+i),i=1,M)/), & ! (/(stuff%matrixc(i)*hc(2*M+1-i)+hc(3*M+1-i),i=M+1,2*M)/)/) M = 2**(n-3) do i = 1, M write(10,'(6a)') ' ',write_var(hs(i)),' = ', & write_var(h(4*i-2)),'+',write_var(h(4*(M+1-i))) end do ! hs = h(2::4)+h(2**(n-1):1:-4) call scaled_odd_tail(hs,scaled_odd_tail_data(n-2),n-2) max_depth = maxval(hs%depth) do i = 1, M hsx(i) = var_stuff(max_depth+1,hs(i)%index,hs(i)%label) end do M = 2**(n-4) do i = 1, M write(10,'(3a,i0,4a)') ' ',write_var(hsx(i)),' = d', & stuff%matrixs(i),'*',write_var(hs(i)),'+', & write_var(hs(M+i)) end do do i = M+1,2*M write(10,'(5a,i0,2a)') ' ',write_var(hsx(i)),' = ', & write_var(hs(2*M+1-i)),'-d', & stuff%matrixs(i),'*',write_var(hs(3*M+1-i)) end do ! hs = (/(/(stuff%matrixs(i)*hs(i)+hs(M+i),i=1,M)/), & ! (/(hs(2*M+1-i)-stuff%matrixs(i)*hs(3*M+1-i),i=M+1,2*M)/)/) ! h2 = h(::2) max_depth = maxval(h2%depth) M = 2**(n-4) do i = 1, M h21(i) = h2(2*i-1) hc1(i) = var_stuff(max_depth+1,h2(2*i)%index,h2(2*i)%label) hs1(i) = var_stuff(max_depth+1,h2(2*(i+M))%index,h2(2*(i+M))%label) h21(i+M) = h2(2*(i+M)-1) end do do i = 1, M write(10,'(6a)') ' ',write_var(hc1(i)),' = ', & write_var(h2(4*i-2)),'-',write_var(h2(4*(M+1-i))) end do ! hc1 = h2(2::4)-h2(2**(n-2):1:-4) call scaled_odd_tail(hc1,scaled_odd_tail_data(n-3),n-3) max_depth = maxval(hc1%depth) do i = 1, M hc1x(i) = var_stuff(max_depth+1,hc1(i)%index,hc1(i)%label) end do M = 2**(n-5) if(n == 4) then write(10,'(3a,i0,2a)') ' ',write_var(hc1x(1)),' = d', & stuff%matrix1c(1,1),'*',write_var(hc1(1)) else do i = 1, M write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(hc1x(i)),' = d', & stuff%matrix1c(1,i),'*',write_var(hc1(i)),'-d', & stuff%matrix1c(2,i),'*',write_var(hc1(M+i)) end do do i = M+1,2*M write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(hc1x(i)),' = d', & stuff%matrix1c(1,i),'*',write_var(hc1(2*M+1-i)),'+d', & stuff%matrix1c(2,i),'*',write_var(hc1(3*M+1-i)) end do end if ! if(n == 4) then ! hc1(1) = stuff%matrix1c(1,1)*hc1(1) ! else ! hc1 = (/(/(stuff%matrix1c(1,i)*hc1(i)- & ! stuff%matrix1c(2,i)*hc1(M+i),i=1,M)/), & ! (/(stuff%matrix1c(1,i)*hc1(2*M+1-i)+ & ! stuff%matrix1c(2,i)*hc1(3*M+1-i),i=M+1,2*M)/)/) ! end if M = 2**(n-4) do i = 1, M write(10,'(6a)') ' ',write_var(hs1(i)),' = ', & write_var(h2(4*i-2)),'+',write_var(h2(4*(M+1-i))) end do ! hs1 = h2(2::4)+h2(2**(n-2):1:-4) call scaled_odd_tail(hs1,scaled_odd_tail_data(n-3),n-3) max_depth = maxval(hs1%depth) do i = 1, M hs1x(i) = var_stuff(max_depth+1,hs1(i)%index,hs1(i)%label) end do M = 2**(n-5) if(n == 4) then write(10,'(3a,i0,2a)') ' ',write_var(hs1x(1)),' = d', & stuff%matrix1s(1,1),'*',write_var(hs1(1)) else do i = 1, M write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(hs1x(i)),' = d', & stuff%matrix1s(1,i),'*',write_var(hs1(i)),'+d', & stuff%matrix1s(2,i),'*',write_var(hs1(M+i)) end do do i = M+1,2*M write(10,'(3a,i0,3a,i0,2a)') ' ',write_var(hs1x(i)),' = d', & stuff%matrix1s(1,i),'*',write_var(hs1(2*M+1-i)),'-d', & stuff%matrix1s(2,i),'*',write_var(hs1(3*M+1-i)) end do end if ! if(n == 4) then ! hs1(1) = stuff%matrix1s(1,1)*hs1(1) ! else ! hs1 = (/(/(stuff%matrix1s(1,i)*hs1(i)+ & ! stuff%matrix1s(2,i)*hs1(M+i),i=1,M)/), & ! (/(stuff%matrix1s(1,i)*hs1(2*M+1-i)- & ! stuff%matrix1s(2,i)*hs1(3*M+1-i),i=M+1,2*M)/)/) ! end if ! h21 = h2(::2) call scaled_odd_tail(h21,scaled_odd_tail_data(n-2),n-2) max_depth = maxval(h21%depth) M = 2**(n-3) do i = 1, M h21x(i) = var_stuff(max_depth+1,h21(i)%index,h21(i)%label) end do do i = 1, M write(10,'(3a,i0,2a)') ' ',write_var(h21x(i)),' = d', & stuff%diag(i),'*',write_var(h21(i)) end do ! h21 = stuff%diag*h21 k = 2**(n-4) max_depth = max(maxval(hc1x%depth),maxval(hs1x%depth),maxval(h21x%depth)) do i = 1, k h2(2*i-1) = var_stuff(max_depth+1,h21x(i)%index,h21x(i)%label) h2(2*i) = var_stuff(max_depth+1,hc1x(i)%index,hc1x(i)%label) h2(2*(i+k)) = var_stuff(max_depth+1,hs1x(i)%index,hs1x(i)%label) h2(2*(i+k)-1) = var_stuff(max_depth+1,h21x(i+k)%index,h21x(i+k)%label) end do do i = 1, k write(10,'(6a)') ' ',write_var(h2(i)),' = ', & write_var(h21x(i)),'+',write_var(hc1x(i)) end do ! h2(0*k+1:1*k) = h21(:k)+hc1 do i = 1, k write(10,'(6a)') ' ',write_var(h2(k+i)),' = ', & write_var(h21x(k+1-i)),'-',write_var(hc1x(k+1-i)) end do ! h2(1*k+1:2*k) = h21(k:1:-1)-hc1(k:1:-1) do i = 1, k write(10,'(6a)') ' ',write_var(h2(2*k+i)),' = ', & write_var(hs1x(i)),'+',write_var(h21x(k+i)) end do ! h2(2*k+1:3*k) = hs1+h21(k+1:) do i = 1, k write(10,'(6a)') ' ',write_var(h2(3*k+i)),' = ', & write_var(hs1x(k+1-i)),'-',write_var(h21x(2*k+1-i)) end do ! h2(3*k+1:4*k) = hs1(k:1:-1)-h21(2*k:k+1:-1) k = 2**(n-3) max_depth = max(maxval(hcx%depth),maxval(hsx%depth),maxval(h2%depth)) do i = 1, k h(2*i-1) = var_stuff(max_depth+1,h2(i)%index,h2(i)%label) h(2*i) = var_stuff(max_depth+1,hcx(i)%index,hcx(i)%label) h(2*(i+k)) = var_stuff(max_depth+1,hsx(i)%index,hsx(i)%label) h(2*(i+k)-1) = var_stuff(max_depth+1,h2(i+k)%index,h2(i+k)%label) end do do i = 1, k write(10,'(6a)') ' ',write_var(h(i)),' = ', & write_var(h2(i)),'+',write_var(hcx(i)) end do ! h(0*k+1:1*k) = h2(:k)+hc do i = 1, k write(10,'(6a)') ' ',write_var(h(k+i)),' = ', & write_var(h2(k+1-i)),'-',write_var(hcx(k+1-i)) end do ! h(1*k+1:2*k) = h2(k:1:-1)-hc(k:1:-1) do i = 1, k write(10,'(6a)') ' ',write_var(h(2*k+i)),' = ', & write_var(hsx(i)),'+',write_var(h2(k+i)) end do ! h(2*k+1:3*k) = hs+h2(k+1:) do i = 1, k write(10,'(6a)') ' ',write_var(h(3*k+i)),' = ', & write_var(hsx(k+1-i)),'-',write_var(h2(2*k+1-i)) end do ! h(3*k+1:4*k) = hs(k:1:-1)-h2(2*k:k+1:-1) end if end subroutine scaled_odd_tail ! Subroutine real_fft_correct ! Calculates real-half complex dft of order 2**n from definition subroutine real_fft_correct(h,n) integer, intent(in) :: n ! real(dp), intent(inout) :: h(0:2**n-1) type(var_stuff), intent(inout) :: h(0:2**n-1) type(var_stuff) h1(0:2**n-1) type(var_stuff) h2(0:2**n-1) integer i, j ! real(dp) pi integer M integer max_depth ! pi = 4*atan(1.0_dp) M = 2**n ! h = (/(/(sum((/(h(j)*cos(2*pi*mod(i*j,M)/M),j=0,M-1)/)),i=0,M/2)/), & ! (/(sum((/(h(j)*sin(2*pi*mod(i*j,M)/M),j=0,M-1)/)),i=1,(M-1)/2)/)/) select case(n) case(1) max_depth = maxval(h%depth) do i = 0, 2**n-1 h1(i) = var_stuff(max_depth+1,h(i)%index,h(i)%label) end do if(max_depth == 0) then write(10,'(6a)') ' ',write_var(h1(0)),' = ', & write_input(h(0)),'+',write_input(h(1)) write(10,'(6a)') ' ',write_var(h1(1)),' = ', & write_input(h(0)),'-',write_input(h(1)) else write(10,'(6a)') ' ',write_var(h1(0)),' = ', & write_var(h(0)),'+',write_var(h(1)) write(10,'(6a)') ' ',write_var(h1(1)),' = ', & write_var(h(0)),'-',write_var(h(1)) end if do i = 0, 2**n-1 h(i) = h1(i) end do case(2) !write(*,*) 'in real_fft_correct, max_depth = ', max_depth max_depth = maxval(h%depth) do i = 0, 2**n-1 h1(i) = var_stuff(max_depth+1,h(i)%index,h(i)%label) h2(i) = var_stuff(max_depth+2,h(i)%index,h(i)%label) end do if(max_depth == 0) then write(10,'(6a)') ' ',write_var(h1(0)),' = ', & write_input(h(0)),'+',write_input(h(2)) write(10,'(6a)') ' ',write_var(h1(2)),' = ', & write_input(h(1)),'+',write_input(h(3)) write(10,'(6a)') ' ',write_var(h2(0)),' = ', & write_var(h1(0)),'+',write_var(h1(2)) write(10,'(6a)') ' ',write_var(h2(1)),' = ', & write_input(h(0)),'-',write_input(h(2)) write(10,'(6a)') ' ',write_var(h2(2)),' = ', & write_var(h1(0)),'-',write_var(h1(2)) write(10,'(6a)') ' ',write_var(h2(3)),' = ', & write_input(h(1)),'-',write_input(h(3)) else write(10,'(6a)') ' ',write_var(h1(0)),' = ', & write_var(h(0)),'+',write_var(h(2)) write(10,'(6a)') ' ',write_var(h1(2)),' = ', & write_var(h(1)),'+',write_var(h(3)) write(10,'(6a)') ' ',write_var(h2(0)),' = ', & write_var(h1(0)),'+',write_var(h1(2)) write(10,'(6a)') ' ',write_var(h2(1)),' = ', & write_var(h(0)),'-',write_var(h(2)) write(10,'(6a)') ' ',write_var(h2(2)),' = ', & write_var(h1(0)),'-',write_var(h1(2)) write(10,'(6a)') ' ',write_var(h2(3)),' = ', & write_var(h(1)),'-',write_var(h(3)) end if do i = 0, 2**n-1 h(i) = h2(i) end do ! case(3) case default stop 'bigger n than 3 in real_fft_correct' end select end subroutine real_fft_correct ! Subroutine real_fft ! Calculates real-half complex fft of order 2**n recursive subroutine real_fft(h,n) integer, intent(in) :: n ! real(dp), intent(inout) :: h(0:2**n-1) type(var_stuff), intent(inout) :: h(0:2**n-1) type(var_stuff) h0(0:2**(n-1)-1) type(var_stuff) h1(0:2**(n-1)-1) integer max_depth integer M integer i if(n < 3) then call real_fft_correct(h,n) else max_depth = maxval(h%depth) if(max_depth == 0) then M = 2**(n-1) do i = 0, 2**(n-1)-1 h0(i) = var_stuff(max_depth+1, h(i)%index, h(i)%label) write(10,'(6a)') ' ',write_var(h0(i)),' = ', & write_input(h(i)),'+',write_input(h(i+M)) end do ! h0 = h(0:2**(n-1)-1)+h(2**(n-1):) call real_fft(h0,n-1) do i = 0, 2**(n-1)-1 h1(i) = var_stuff(max_depth+1, h(i+M)%index, h(i+M)%label) write(10,'(6a)') ' ',write_var(h1(i)),' = ', & write_input(h(i)),'-',write_input(h(i+M)) end do ! h1 = h(0:2**(n-1)-1)-h(2**(n-1):) else M = 2**(n-1) do i = 0, 2**(n-1)-1 h0(i) = var_stuff(max_depth+1, h(i)%index, h(i)%label) write(10,'(6a)') ' ',write_var(h0(i)),' = ', & write_var(h(i)),'+',write_var(h(i+M)) end do ! h0 = h(0:2**(n-1)-1)+h(2**(n-1):) call real_fft(h0,n-1) do i = 0, 2**(n-1)-1 h1(i) = var_stuff(max_depth+1, h(i+M)%index, h(i+M)%label) write(10,'(6a)') ' ',write_var(h1(i)),' = ', & write_var(h(i)),'-',write_var(h(i+M)) end do ! h1 = h(0:2**(n-1)-1)-h(2**(n-1):) end if call odd_tail(h1,odd_tail_data(n),n) do i = 0, M-1 h(2*i) = h0(i) h(2*i+1) = h1(i) end do ! h(0::2) = h0 ! h(1::2) = h1 end if end subroutine real_fft ! Subroutine complex_fft ! Computes full complex fft of order 2**n ! subroutine complex_fft(h,n) subroutine complex_fft(n) integer, intent(in) :: n ! type(var_stuff), intent(out) :: h(0:2**n-1) type(var_stuff) hreal(0:2**n-1) type(var_stuff) himag(0:2**n-1) integer M integer i M = 2**n ! hreal = real(h) do i = 0, 2**n-1 hreal(i) = var_stuff(0,i,'xr') end do call real_fft(hreal,n) ! himag = aimag(h) do i = 0, 2**n-1 himag(i) = var_stuff(0,i,'xi') end do call real_fft(himag,n) ! h(0) = cmplx(hreal(0),himag(0),dp) write(10,'(5a)') ' h(0) = cmplx(', & write_var(hreal(0)),',',write_var(himag(0)),',wp)' do i = 1, (M-1)/2 write(10,'(a,i0,9a)') ' h(',i,') = cmplx(', & write_var(hreal(i)),'+',write_var(himag(M/2+i)),',', & write_var(himag(i)),'-',write_var(hreal(M/2+i)),',wp)' ! h(i) = cmplx(hreal(i)-himag(M/2+i),hreal(M/2+i)+himag(i),dp) write(10,'(a,i0,9a)') ' h(',M-i,') = cmplx(', & write_var(hreal(i)),'-',write_var(himag(M/2+i)),',', & write_var(hreal(M/2+i)),'+',write_var(himag(i)),',wp)' ! h(M-i) = cmplx(hreal(i)+himag(M/2+i),-hreal(M/2+i)+himag(i),dp) end do write(10,'(a,i0,5a)') ' h(',M/2,') = cmplx(', & write_var(hreal(M/2)),',',write_var(himag(M/2)),',wp)' ! h(M/2) = cmplx(hreal(M/2),himag(M/2),dp) end subroutine complex_fft function write_var(var) type(var_stuff), intent(in) :: var character(write_var_len(var)) write_var write(write_var,'(a,i0,a,i0)') trim(adjustl(var%label)), & var%depth, '_', var%index end function write_var pure function write_var_len(var) integer write_var_len type(var_stuff), intent(in) :: var character(80) test_string write(test_string,'(a,i0,a,i0)') trim(adjustl(var%label)), & var%depth, '_', var%index write_var_len = len_trim(test_string) end function write_var_len function write_input(var) type(var_stuff), intent(in) :: var character(write_input_len(var)) write_input if(scan(var%label,'iI') == 0) then write(write_input,'(a,i0,a)') 'real(h(',var%index,'))' else write(write_input,'(a,i0,a)') 'aimag(h(',var%index,'))' end if end function write_input pure function write_input_len(var) type(var_stuff), intent(in) :: var integer write_input_len character(80) test_string if(scan(var%label,'iI') == 0) then write(test_string,'(a,i0,a)') 'real(h(',var%index,'))' else write(test_string,'(a,i0,a)') 'aimag(h(',var%index,'))' end if write_input_len = len_trim(test_string) end function write_input_len end module correct program test use mykinds use correct use avl_mod implicit none integer n ! complex(dp), allocatable :: h1(:) ! complex(dp), allocatable :: h2(:) ! real(dp), allocatable :: h1(:) ! real(dp), allocatable :: h2(:) ! type(var_stuff), allocatable :: h2(:) integer i integer j integer k integer :: ulps = 100 integer M real(dp) pi real(dp) t0, t1 real(qp), allocatable :: qconst_array(:) character(80) filename character(8) date character(10) time integer An, Mn call cpu_time(t0) pi = 4*atan(1.0_dp) write(*,'(a)',advance='no') ' Enter the exponent of the problem:> ' read(*,*) n write(filename,'(a,i0,a)') 'fft', 2**n, '.f90' open(10,file=filename,status='replace') !************************************* write(10,'(a,a)') '! File: ', trim(filename) call date_and_time(date, time) write(10,'(12a)') '! Generated by scaled2d.f90 ', date(5:6), & '/', date(7:8), '/', date(1:4), ' ', time(1:2), & ':', time(3:4), ':', time(4:10) write(10,'(a)') 'module mykinds' write(10,'(a)') ' implicit none' write(10,'(a)') ' integer, parameter :: sp = selected_real_kind(6,30)' write(10,'(a)') ' integer, parameter :: dp = selected_real_kind(15,300)' write(10,'(a)') 'end module mykinds' write(10,'(a)') '' write(10,'(a)') 'module dp_mod' write(10,'(a)') ' use mykinds, only: wp => dp' write(10,'(a)') ' implicit none' write(10,'(a)') ' private wp' write(10,'(a)') ' contains' write(10,'(a,i0,a)') 'include ''fft',2**n,'.i90''' write(10,'(a)') 'end module dp_mod' write(10,'(a)') '' write(10,'(a)') 'program test' write(10,'(a)') ' use mykinds' write(10,'(a)') ' use dp_mod' write(10,'(a)') ' implicit none' write(10,'(a,i0)') ' integer, parameter :: N = ',2**n write(10,'(a)') ' complex(dp) h1(0:N-1)' write(10,'(a)') ' complex(dp) h2(0:N-1)' write(10,'(a)') ' complex(dp), parameter :: ri(2) = (/(1,0),(0,1)/)' write(10,'(a)') ' integer i, j, k, L' write(10,'(a)') ' real(dp) pi' write(10,'(a)') ' real(dp) err, maxerr' write(10,'(a)') '' write(10,'(a)') ' pi = 4*atan(1.0_dp)' write(10,'(a)') ' maxerr = 0' write(10,'(a)') ' do i = 0, N-1' write(10,'(a)') ' do j = 1, 2' write(10,'(a)') ' h2 = 0' write(10,'(a)') ' h2(i) = ri(j)' ! write(10,'(a)') ' h1 = (/(sum((/(exp(-2*pi*(0,1)*mod(k*L,N)/N)*'// & ! 'h2(L),L=0,N-1)/)),k=0,N-1)/)' write(10,'(a)') ' h1 = (/(h2(i)*exp(-2*pi*(0,1)*mod(i*k,N)/N),k=0,N-1)/)' write(10,'(a,i0,a)') ' call fft',2**n,'(h2)' write(10,'(a)') ' err = max(maxval(abs(real(h1-h2))),maxval(abs(aimag(h1-h2))))' write(10,'(a)') ' write(*,*) i, j, err' write(10,'(a)') ' maxerr = max(err, maxerr)' write(10,'(a)') ' end do' write(10,'(a)') ' end do' write(10,'(a)') ' write(*,*)' write(10,'(a)') ' write(*,*) maxerr' write(10,'(a)') 'end program test' !************************************* close(10) write(filename,'(a,i0,a)') 'fft', 2**n, '.i90' open(10,file=filename,status='replace') write(10,'(a,a)') '! File: ', trim(filename) call date_and_time(date, time) write(10,'(12a)') '! Generated by scaled2d.f90 ', date(5:6), & '/', date(7:8), '/', date(1:4), ' ', time(1:2), & ':', time(3:4), ':', time(4:10) write(10,'(a,i0,a)') 'subroutine fft', 2**n, '(h)' An = (8*3*n*2**n-16*2**n-2*(-1)**n)/9+2 Mn = (3*10*n*2**n-76*2**n-2*3*n*(-1)**n+22*(-1)**n)/27-2*n+6 write(10,'(a,i0,a,i0,a)') '! ',An,' adds, ',Mn,' muls' write(10,'(a)') ' implicit real(wp) (x,y)' write(10,'(a,i0,a)') ' complex(wp), intent(inout) :: h(0:', 2**n-1, ')' allocate(get_factor_data(n-2)) do i = 1, n-2 allocate(get_factor_data(i)%factor(2**(i-4))) if(i >= 4) then call get_factor(get_factor_data(i)%factor,i,get_factor_data(i-2)%factor) end if end do call prepare_correct(n) allocate(odd_tail_data(n)) do i = 1, n call odd_tail_gen(odd_tail_data(i),i) end do allocate(scaled_odd_tail_data(n-2)) do i = 1, n-2 call scaled_odd_tail_gen(scaled_odd_tail_data(i),i) end do allocate(qconst_array(body%num_nodes)) call list2array(body, qconst_array) allocate(const_array(size(qconst_array))) const_array = qconst_array M = 2**n write(*,'()') ! allocate(h1(M)) ! allocate(h2(M)) call cpu_time(t1) do i = 1, size(qconst_array) write(10,'(a,i0,a,f0.34,a)') ' real(wp), parameter :: d', & i, ' = ', qconst_array(i), '_wp' end do deallocate(qconst_array) write(*,'(a,f0.3)') ' Prep time = ', t1-t0 call cpu_time(t0) ! h1 = (/(h1(i)*exp(-2*pi*(0,1)*mod((i-1)*j,M)/M),j=0,M-1)/) ! h1 = (/(/(h1(i)*cos(2*pi*mod((i-1)*j,M)/M),j=0,M/2)/), & ! (/(h1(i)*sin(2*pi*mod((i-1)*j,M)/M),j=1,M/2-1)/)/) call complex_fft(n) ! call real_fft(h2,n) call cpu_time(t1) write(*,'(a,f0.3)') ' Execution time = ', t1-t0 write(10,'(a,i0)') 'end subroutine fft', 2**n end program test