! File: scaled2h.f90 ! Public domain 2004 James Van Buskirk ! PURPOSE: Generates full complex FFT algorithms of order 2**n. ! Prompts the user for the exponent of the problem, and if the ! user wants to generate FFT code for problems of order 2**n, he ! enters n at the prompt. Two files are then created: ! 1) fft###.i90 contains a subroutine fft###(h) that converts a ! complex vector h of order ### = 2**n into its DFT: ! h(0:2**n-1) = & ! (/(sum((/(exp(-2*pi*(0,1)*i*j/2**n),j=0,2**n-1)/)),i=0,2**n-1)/) ! 2) fft###.f90 contains a test program that shows how one has ! to INCLUDE the above file in a module that defines the ! parameter wp to be a kind type parameter suitable for real ! types, and how to call fft### to convert a vector into its ! DFT. ! If everything works OK, the fft###.f90 file should compile ! to a program that prints out 2*2**n lines of output that are ! near EPSILON(1.0d0) and a final line of output that is the ! maximum of the previous 2*2**n lines. ! *************************************************************** ! module mykinds ! Sets up kind type parameter for double precision and ! (if available) quadruple precision. If there is no ! quad precision type, double precision will be used. ! *************************************************************** 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 ! Contains tools and derived types to store randomly generated ! data in an AVL tree and convert the tree to an array. ! NOTE ABOUT RE-USE: This is not a general-purpose AVL tree ! module: it was designed to work on only positive data, and ! there is no deletion procedure nor would it be simple to ! implement one for this module without adding an extra data ! structure to handle the indices. ! *************************************************************** module avl_mod use mykinds implicit none !................................................................ ! type avl_node ! A node in the AVL tree like any avl node, but also has ! an index that tells us where the key will end up in ! the eventual array. !................................................................ type avl_node real(qp) key ! Actual (positive) data integer index ! Index into ultimate array integer :: balance = 0 ! -1 = heavy left, +1 = heavey right type(avl_node), pointer :: left => NULL() ! Pointer to left child type(avl_node), pointer :: right => NULL() ! Pointer to right child end type avl_node !................................................................ !................................................................ ! type avl_bag ! Contains an AVL tree root and a counter of its number of nodes !................................................................ type avl_bag type(avl_node), pointer :: root => NULL() ! Root of AVL tree integer :: num_nodes = 0 ! Node counter end type avl_bag !................................................................ !................................................................ ! Interface operator(<) ! Allows us to fuzzy compare a real(qp) with an avl_node !................................................................ interface operator(<) module procedure lessn end interface operator(<) !................................................................ !................................................................ ! Interface operator(>) ! Allows us to fuzzy compare a real(qp) with an avl_node !................................................................ interface operator(>) module procedure greatern end interface operator(>) !................................................................ ! Fuzzy compares above will be up to 100 ulps integer, parameter :: ulps = 100 ! Module entities private by default private ! Need the below public avl_node ! Need the above ! Only public entities: public avl_bag, key2index, bag2array contains !................................................................ ! function lessn(key,root) ! Performs fuzzy compare of real(qp) key with type(avl_tree) root !................................................................ pure function lessn(key,root) real(qp), intent(in) :: key ! Left operand of comparison type(avl_node), intent(in) :: root ! Right operand of comparison logical lessn ! Result of comparison lessn = abs(key) < abs(root%key)-ulps*spacing(root%key) end function lessn !................................................................ !................................................................ ! function greatern(key,root) ! Performs fuzzy compare of real(qp) key with type(avl_tree) root !................................................................ pure function greatern(key,root) real(qp), intent(in) :: key ! Left operand of comparison type(avl_node), intent(in) :: root ! Right operand of comparison logical greatern ! Result of comparison greatern = abs(key) > abs(root%key)+ulps*spacing(root%key) end function greatern !................................................................ !................................................................ ! function key2index(key,bag,index) ! Stores key as avl_node in bag and returns index of key ! in eventual array. Does not store duplicate copies. !................................................................ subroutine key2index(key,bag,index) real(qp), intent(in) :: key ! Key to be compared and stored type(avl_bag), intent(inout) :: bag ! Where the key will be stored integer, intent(out) :: index ! Index of incoming key logical changed ! Required for avl_insert; ignored index = bag%num_nodes if(associated(bag%root)) then ! We have an AVL tree to search call avl_insert(key,bag%root,index,changed) bag%num_nodes = max(bag%num_nodes, index) else ! Start new AVL tree bag%num_nodes = 1 allocate(bag%root) bag%root%key = key index = 1 bag%root%index = index end if end subroutine key2index !................................................................ !................................................................ ! subroutine avl_insert(key,root,index,changed) ! Searched AVL tree starting at root for key. If found, returns ! index of matching key and clears changed. If not found, creates ! new node containing key and next index. Sets change accordingly. !................................................................ recursive subroutine avl_insert(key,root,index,changed) real(qp), intent(in) :: key ! Key to search for in AVL tree type(avl_node), pointer :: root ! Root of AVL tree (must be associated) integer, intent(inout) :: index ! On input, last valid index ! On output, index of key logical, intent(out) :: changed ! Set if tree grew longer if(key < root) then if(associated(root%left)) then call avl_insert(key,root%left,index,changed) ! After insertion, if tree grew longer we must recalculate ! balance factor and rotate to maintain AVL condition if ! necessary. if(changed) then if(root%balance == -1) then if(root%left%balance == 1) then call rot2right(root) else ! root%left%balance <= 0 call rot1right(root) end if changed = .FALSE. else if(root%balance == 0) then root%balance = -1 else ! root%balance == 1 root%balance = 0 changed = .FALSE. end if end if else ! .NOT.associated(root%left) ! If there is an empty subtree to search, we know the key ! goes here. index = index+1 allocate(root%left) root%left%key = key root%left%index = index changed = root%balance == 0 root%balance = root%balance-1 end if else if(key > root) then if(associated(root%right)) then call avl_insert(key,root%right,index,changed) ! After insertion, if tree grew longer we must recalculate ! balance factor and rotate to maintain AVL condition if ! necessary. if(changed) then if(root%balance == -1) then root%balance = 0 changed = .FALSE. else if(root%balance == 0) then root%balance = 1 else ! root%balance == 1 if(root%right%balance == -1) then call rot2left(root) else ! root%right%balance <= 0 call rot1left(root) end if changed = .FALSE. end if end if else ! .NOT.associated(root%right) ! If there is an empty subtree to search, we know the key ! goes here. index = index+1 allocate(root%right) root%right%key = key root%right%index = index changed = root%balance == 0 root%balance = root%balance+1 end if else ! key == root ! If we found the key, return the matching index. changed = .FALSE. index = root%index end if end subroutine avl_insert !................................................................ !................................................................ ! Subroutine rot2right(root) ! Performs a double right rotation of the AVL tree at root !................................................................ subroutine rot2right(root) type(avl_node), pointer :: root ! Root of AVL tree type(avl_node), pointer :: rotate_temp ! Temporary pointer to new root ! Create new connections rotate_temp => root%left%right root%left%right => rotate_temp%left rotate_temp%left => root%left root%left => rotate_temp%right rotate_temp%right => root root => rotate_temp ! Recalculate balance factors root%right%balance = max(-root%balance, 0) root%left%balance = -(root%balance+root%right%balance) root%balance = 0 end subroutine rot2right !................................................................ !................................................................ ! Subroutine rot2left(root) ! Performs a double left rotation of the AVL tree at root !................................................................ subroutine rot2left(root) type(avl_node), pointer :: root ! Root of AVL tree type(avl_node), pointer :: rotate_temp ! Temporary pointer to new root ! Create new connections rotate_temp => root%right%left root%right%left => rotate_temp%right rotate_temp%right => root%right root%right => rotate_temp%left rotate_temp%left => root root => rotate_temp ! Recalculate balance factors root%left%balance = min(-root%balance, 0) root%right%balance = -(root%balance+root%left%balance) root%balance = 0 end subroutine rot2left !................................................................ !................................................................ ! Subroutine rot1right(root) ! Performs a single right rotation of the AVL tree at root !................................................................ subroutine rot1right(root) type(avl_node), pointer :: root ! Root of AVL tree type(avl_node), pointer :: rotate_temp ! Temporary pointer to new root ! Create new connections rotate_temp => root%left root%left => rotate_temp%right rotate_temp%right => root root => rotate_temp ! Recalculate balance factors root%balance = root%balance+1 root%right%balance = -root%balance end subroutine rot1right !................................................................ !................................................................ ! Subroutine rot1left(root) ! Performs a single left rotation of the AVL tree at root !................................................................ subroutine rot1left(root) type(avl_node), pointer :: root ! Root of AVL tree type(avl_node), pointer :: rotate_temp ! Temporary pointer to new root ! Create new connections rotate_temp => root%right root%right => rotate_temp%left rotate_temp%left => root root => rotate_temp ! Recalculate balance factors root%balance = root%balance-1 root%left%balance = -root%balance end subroutine rot1left !................................................................ !................................................................ ! subroutine bag2array(bag,array) ! Destroys AVL tree rooted in bag and stores the keys at ! the appropriate indices in array !................................................................ subroutine bag2array(bag, array) type(avl_bag), intent(inout) :: bag ! Has tree to be converted real(qp), intent(out) :: array(bag%num_nodes) ! Resultant array of keys ! Set up call to first level of avl_detroy if(associated(bag%root)) then call avl_destroy(bag%root,array) bag%num_nodes = 0 end if end subroutine bag2array !................................................................ !................................................................ ! subroutine avl_destroy(root,array) ! Destroys AVL tree rooted at root and stores keys at the ! appropriate indices in array. !................................................................ recursive subroutine avl_destroy(root,array) type(avl_node), pointer :: root ! Root of AVL tree to destroy real(qp), intent(inout) :: array(*) ! Array to store contents of tree ! Set up call to next level on left if(associated(root%left)) then call avl_destroy(root%left,array) end if ! Set up call to next level on right if(associated(root%right)) then call avl_destroy(root%right,array) end if ! Store current key in array in destroy current node array(root%index) = root%key deallocate(root) end subroutine avl_destroy !................................................................ end module avl_mod ! *************************************************************** ! *************************************************************** ! module code_gen ! Contains tools and derived types to generate full complex DFT code ! *************************************************************** module code_gen use mykinds use avl_mod implicit none !................................................................ ! type get_factor_stuff ! Contains a pointer to a quad-precision array !................................................................ type get_factor_stuff real(qp), pointer :: factor(:) ! Pointer to quad array for current order end type get_factor_stuff !................................................................ !................................................................ ! type odd_tail_stuff ! Contains pointers to 2-d arrays for sine and cosine array entries ! Entries will contain the indices of the constants !................................................................ type odd_tail_stuff integer, pointer :: matrixc(:,:) ! Matrix of cosine coefficient indices integer, pointer :: matrixs(:,:) ! Matrix of sine coefficient indices end type odd_tail_stuff !................................................................ !................................................................ ! type scaled_odd_tail_stuff ! Contains pointers to arrays of tangents, sine and cosine entries, ! and the last diagonal. Entries will contain the indices of the constants !................................................................ type scaled_odd_tail_stuff integer, pointer :: matrixc(:) ! Array of cosine tangent indices integer, pointer :: matrixs(:) ! Array of sine tangent indices integer, pointer :: matrix1c(:,:) ! Matrix of cosine coefficient indices integer, pointer :: matrix1s(:,:) ! Matrix of sin coefficient indices integer, pointer :: diag(:) ! Array of last diagonal indices end type scaled_odd_tail_stuff !................................................................ !................................................................ ! type correct_stuff ! Contains data to compute odd_tail(8). Entries ! will contain the indices of the constants. !................................................................ type correct_stuff integer rsq2 ! Index of 1/sqrt(2) end type correct_stuff !................................................................ !................................................................ ! type var_stuff ! Contains latitude, longtude, and label of current variable. ! Note LF95 pads at the end to fit into arrays !................................................................ type var_stuff integer depth ! Number of LOC executed to reach variable integer index ! Index of variable among current set character(2) label ! xr = real even, yr = real odd ! xi = imaginary even, yi = imaginary odd ! 'y' labels not yet implemented end type var_stuff !................................................................ contains !................................................................ ! Subroutine get_factor(factor,n,old_factor) ! Computes scale factor thrown out by scaled_odd_tail(n) !................................................................ recursive subroutine get_factor(factor, n, old_factor) integer, intent(in) :: n ! Order of the problem real(qp), intent(out) :: factor(2**(n-4)) ! Scale factor real(qp), intent(in) :: old_factor(2**(n-6)) ! Previous scale factor integer i ! Multiple of fundamental angle integer k ! Current block in factor 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,body,correct_data) ! Prepares constants for odd_tail algorithms leading to degenerate ! recursion !................................................................ subroutine prepare_correct(n, body, correct_data) type(avl_bag), intent(inout) :: body ! Holds any keys created integer, intent(in) :: n ! Order of overall problem type(correct_stuff) correct_data ! Holds indices of keys created real(qp) constant ! Current constant created real(qp) pi ! Pi if(n > 2) then pi = 4*atan(1.0_qp) constant = cos(pi/4) call key2index(constant, body, correct_data%rsq2) end if end subroutine prepare_correct !................................................................ !................................................................ ! Subroutine odd_tail_correct(h,n,correct_data,iunit) ! Computes odd_tail(n)*h using the theoretical formula for ! orders too degenerate for recursion and immotalizes the ! result in code !................................................................ subroutine odd_tail_correct(h,n,correct_data,iunit) integer, intent(in) :: n ! Order of the problem type(var_stuff), intent(inout) :: h(2**(n-1)) ! Array of variable mnemonics type(correct_stuff) correct_data ! Data for degenerate problem integer, intent(in) :: iunit ! Write code to this unit type(var_stuff) h1(2**(n-1)), h2(2**(n-1)) ! Temporary variable mnemonics integer max_depth ! Previous maximum depth integer i ! General index 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) call WX(' '//write_var(h1(1))//' = '//write_const(correct_data%rsq2)// & '*('//write_var(h(2))//'-'//write_var(h(4))//')',iunit) ! xc = const_array(correct_data%rsq2)*(h(2)-h(4)) h1(2) = var_stuff(max_depth+1, h(2)%index, h(2)%label) call WX(' '//write_var(h1(2))//' = '//write_const(correct_data%rsq2)// & '*('//write_var(h(2))//'+'//write_var(h(4))//')',iunit) ! 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 call WX(' '//write_var(h2(1))//' = '// & write_var(h(1))//'+'//write_var(h1(1)),iunit) call WX(' '//write_var(h2(2))//' = '// & write_var(h(1))//'-'//write_var(h1(1)),iunit) call WX(' '//write_var(h2(3))//' = '// & write_var(h1(2))//'+'//write_var(h(3)),iunit) call WX(' '//write_var(h2(4))//' = '// & write_var(h1(2))//'-'//write_var(h(3)),iunit) ! 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 4 in odd_tail_correct' end select end subroutine odd_tail_correct !................................................................ !................................................................ ! Subroutine odd_tail_gen(stuff,n,body,get_factor_data) ! Compute odd_tail constants and places their indices in arrays !................................................................ subroutine odd_tail_gen(stuff,n,body,get_factor_data) integer, intent(in) :: n ! Order of the problem type(avl_bag) body ! AVL tree to hold constants type(odd_tail_stuff), intent(out) :: stuff ! Indices of data for ! current order type(get_factor_stuff), intent(in) :: get_factor_data(n-2) ! Factors ! to absorb integer i ! Current row integer j ! Column select integer k ! Current blocksize real(qp) pi ! Pi real(qp) matrix(2,2**(n-3)) ! Array of actual data 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 i = 1, k matrix(:,0*k+i) = matrix(:,0*k+i)* & get_factor_data(n-2)%factor(i) matrix(:,8*k+1-i) = matrix(:,8*k+1-i)* & get_factor_data(n-2)%factor(i) matrix(:,1*k+i) = matrix(:,1*k+i)* & get_factor_data(n-2)%factor(k+1-i) matrix(:,7*k+1-i) = matrix(:,7*k+1-i)* & get_factor_data(n-2)%factor(k+1-i) matrix(:,2*k+i) = matrix(:,2*k+i)* & get_factor_data(n-2)%factor(i) matrix(:,6*k+1-i) = matrix(:,6*k+1-i)* & get_factor_data(n-2)%factor(i) matrix(:,3*k+i) = matrix(:,3*k+i)* & get_factor_data(n-2)%factor(k+1-i) matrix(:,5*k+1-i) = matrix(:,5*k+1-i)* & get_factor_data(n-2)%factor(k+1-i) end do do i = 1, 2**(n-3) do j = 1, 2 call key2index(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 i = 1, k matrix(:,0*k+i) = matrix(:,0*k+i)* & get_factor_data(n-2)%factor(i) matrix(:,8*k+1-i) = matrix(:,8*k+1-i)* & get_factor_data(n-2)%factor(i) matrix(:,1*k+i) = matrix(:,1*k+i)* & get_factor_data(n-2)%factor(k+1-i) matrix(:,7*k+1-i) = matrix(:,7*k+1-i)* & get_factor_data(n-2)%factor(k+1-i) matrix(:,2*k+i) = matrix(:,2*k+i)* & get_factor_data(n-2)%factor(i) matrix(:,6*k+1-i) = matrix(:,6*k+1-i)* & get_factor_data(n-2)%factor(i) matrix(:,3*k+i) = matrix(:,3*k+i)* & get_factor_data(n-2)%factor(k+1-i) matrix(:,5*k+1-i) = matrix(:,5*k+1-i)* & get_factor_data(n-2)%factor(k+1-i) end do do i = 1, 2**(n-3) do j = 1, 2 call key2index(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,correct_data, & ! odd_tail_data,scaled_odd_tail_data,iunit) ! Compute action of odd tail matrix of order 2**(n-1) ! and immortalize in code !................................................................ recursive subroutine odd_tail(h,n,correct_data, & odd_tail_data,scaled_odd_tail_data,iunit) integer, intent(in) :: n ! Order of the problem type(var_stuff), intent(inout) :: h(2**(n-1)) ! Array of variable mnemonics type(correct_stuff) correct_data ! Indices of data for ! degenerate orders type(odd_tail_stuff), intent(in) :: odd_tail_data(n) ! Indices of data ! For the problem type(scaled_odd_tail_stuff) scaled_odd_tail_data(n-2) ! Indices of data for ! scaled version of problem integer, intent(in) :: iunit ! Write code to this unit type(var_stuff) hc(2**(n-3)), hc1(2**(n-3)) ! Temporary varaible mnemonics type(var_stuff) hs(2**(n-3)), hs1(2**(n-3)) ! Temporary varaible mnemonics type(var_stuff) h2(2**(n-2)) ! Temporary varaible mnemonics integer i ! Row number integer k ! Blocksize integer M ! Blocksize integer max_depth ! Previous maximum depth if(n < 4) then call odd_tail_correct(h,n,correct_data,iunit) 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 call WX(' '//write_var(hc(i))//' = '// & write_var(h(4*i-2))//'-'//write_var(h(4*(M+1-i))),iunit) end do ! hc = h(2::4)-h(2**(n-1):1:-4) call scaled_odd_tail(hc,n-2,correct_data,scaled_odd_tail_data,iunit) 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 call WX(' '//write_var(hc1(i))//' = '// & write_const(odd_tail_data(n)%matrixc(1,i))//'*'// & write_var(hc(i))//'-'//write_const(odd_tail_data(n)%matrixc(2,i)) & //'*'//write_var(hc(M+i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hc1(i))//' = '// & write_const(odd_tail_data(n)%matrixc(1,i))//'*'// & write_var(hc(2*M+1-i))//'+'// & write_const(odd_tail_data(n)%matrixc(2,i))//'*'// & write_var(hc(3*M+1-i)),iunit) end do ! hc = (/(/(odd_tail_data(n)%matrixc(1,i)*hc(i)- & ! odd_tal_data(n)%matrixc(2,i)*hc(M+i),i=1,M)/), & ! (/(odd_tail_data(n)%matrixc(1,i)*hc(2*M+1-i)+ & ! odd_tail_data(n)%matrixc(2,i)*hc(3*M+1-i),i=M+1,2*M)/)/) M = 2**(n-3) do i = 1, M call WX(' '//write_var(hs(i))//' = '// & write_var(h(4*i-2))//'+'//write_var(h(4*(M+1-i))),iunit) end do ! hs = h(2::4)+h(2**(n-1):1:-4) call scaled_odd_tail(hs,n-2,correct_data,scaled_odd_tail_data,iunit) 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 call WX(' '//write_var(hs1(i))//' = '// & write_const(odd_tail_data(n)%matrixs(1,i))//'*'// & write_var(hs(i))//'+'// & write_const(odd_tail_data(n)%matrixs(2,i))//'*'// & write_var(hs(M+i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hs1(i))//' = '// & write_const(odd_tail_data(n)%matrixs(1,i))//'*'// & write_var(hs(2*M+1-i))//'-'// & write_const(odd_tail_data(n)%matrixs(2,i))//'*'// & write_var(hs(3*M+1-i)),iunit) end do ! hs = (/(/(odd_tail_data(n)%matrixs(1,i)*hs(i)+ & ! odd_tail_data(n)%matrixs(2,i)*hs(M+i),i=1,M)/), & ! (/(odd_tail_data(n)%matrixs(1,i)*hs(2*M+1-i)- & ! odd_tail_data(n)%matrixs(2,i)*hs(3*M+1-i),i=M+1,2*M)/)/) ! h2 = h(::2) call odd_tail(h2,n-1,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) 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 call WX(' '//write_var(h(i))//' = '// & write_var(h2(i))//'+'//write_var(hc1(i)),iunit) end do ! h(0*k+1:1*k) = h2(:k)+hc do i = 1, k call WX(' '//write_var(h(k+i))//' = '// & write_var(h2(k+1-i))//'-'//write_var(hc1(k+1-i)),iunit) end do ! h(1*k+1:2*k) = h2(k:1:-1)-hc(k:1:-1) do i = 1, k call WX(' '//write_var(h(2*k+i))//' = '// & write_var(hs1(i))//'+'//write_var(h2(k+i)),iunit) end do ! h(2*k+1:3*k) = hs+h2(k+1:) do i = 1, k call WX(' '//write_var(h(3*k+i))//' = '// & write_var(hs1(k+1-i))//'-'//write_var(h2(2*k+1-i)),iunit) 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(stuff,n,body,get_factor_data) ! Generates constants for scaled_odd_tail and store their indices !................................................................ subroutine scaled_odd_tail_gen(stuff,n,body,get_factor_data) integer, intent(in) :: n ! Order of the problem type(avl_bag) body ! Holds actual data type(scaled_odd_tail_stuff), intent(out) :: stuff ! Holds ! Indices of coefficients type(get_factor_stuff), intent(in) :: get_factor_data(n) ! Factors ! thrown out by scaled_odd_tail real(qp) factor(2**(n-4)) ! Factors passed through integer i ! Row number integer j ! Column select integer k ! Blocksize real(qp) pi ! Pi real(qp) matrix(2**(n-3)) ! Actual array of coefficients real(qp) matrix1(min(2,2**(n-4)),2**(n-4)) ! Actual matrix of coefficients real(qp) diag(2**(n-3)) ! Actual array of coefficients 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 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 key2index(abs(matrix(i)), body, stuff%matrixc(i)) end do 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 key2index(abs(matrix(i)), body, stuff%matrixs(i)) end do 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 i = 1, k matrix1(:,0*k+i) = matrix1(:,0*k+i)* & get_factor_data(n-3)%factor(i) matrix1(:,8*k+1-i) = matrix1(:,8*k+1-i)* & get_factor_data(n-3)%factor(i) matrix1(:,1*k+i) = matrix1(:,1*k+i)* & get_factor_data(n-3)%factor(k+1-i) matrix1(:,7*k+1-i) = matrix1(:,7*k+1-i)* & get_factor_data(n-3)%factor(k+1-i) matrix1(:,2*k+i) = matrix1(:,2*k+i)* & get_factor_data(n-3)%factor(i) matrix1(:,6*k+1-i) = matrix1(:,6*k+1-i)* & get_factor_data(n-3)%factor(i) matrix1(:,3*k+i) = matrix1(:,3*k+i)* & get_factor_data(n-3)%factor(k+1-i) matrix1(:,5*k+1-i) = matrix1(:,5*k+1-i)* & get_factor_data(n-3)%factor(k+1-i) end do k = 2**(n-4) do i = 1, k matrix1(:,0*k+i) = matrix1(:,0*k+i)/factor(i) end do do i = 1, 2**(n-4) do j = 1, min(2,2**(n-4)) call key2index(abs(matrix1(j,i)), body, stuff%matrix1c(j,i)) end do end do 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 i = 1, k matrix1(:,0*k+i) = matrix1(:,0*k+i)* & get_factor_data(n-3)%factor(i) matrix1(:,8*k+1-i) = matrix1(:,8*k+1-i)* & get_factor_data(n-3)%factor(i) matrix1(:,1*k+i) = matrix1(:,1*k+i)* & get_factor_data(n-3)%factor(k+1-i) matrix1(:,7*k+1-i) = matrix1(:,7*k+1-i)* & get_factor_data(n-3)%factor(k+1-i) matrix1(:,2*k+i) = matrix1(:,2*k+i)* & get_factor_data(n-3)%factor(i) matrix1(:,6*k+1-i) = matrix1(:,6*k+1-i)* & get_factor_data(n-3)%factor(i) matrix1(:,3*k+i) = matrix1(:,3*k+i)* & get_factor_data(n-3)%factor(k+1-i) matrix1(:,5*k+1-i) = matrix1(:,5*k+1-i)* & get_factor_data(n-3)%factor(k+1-i) end do k = 2**(n-4) do i = 1, k matrix1(:,0*k+i) = matrix1(:,0*k+i)/factor(i) end do do i = 1, 2**(n-4) do j = 1, min(2,2**(n-4)) call key2index(abs(matrix1(j,i)), body, stuff%matrix1s(j,i)) end do end do k = 2**(n-6) diag = 1 do i = 1, k diag(0*k+i) = get_factor_data(n-2)%factor(i) diag(1*k+i) = get_factor_data(n-2)%factor(k+1-i) diag(2*k+i) = get_factor_data(n-2)%factor(i) diag(3*k+i) = get_factor_data(n-2)%factor(k+1-i) diag(4*k+i) = get_factor_data(n-2)%factor(i) diag(5*k+i) = get_factor_data(n-2)%factor(k+1-i) diag(6*k+i) = get_factor_data(n-2)%factor(i) diag(7*k+i) = get_factor_data(n-2)%factor(k+1-i) end do k = 2**(n-4) do i = 1, k diag(0*k+i) = diag(0*k+i)/factor(i) diag(1*k+i) = diag(1*k+i)/factor(i) end do do i = 1, 2**(n-3) call key2index(abs(diag(i)), body, stuff%diag(i)) end do end if end subroutine scaled_odd_tail_gen !................................................................ ! Subroutine scaled_odd_tail(h,n,correct_data,scaled_odd_tail_data,iunit) ! Computes action of scaled odd tail matrix of order 2**(n-1) and ! immortalizes in code. !................................................................ recursive subroutine scaled_odd_tail(h,n,correct_data,scaled_odd_tail_data,iunit) integer, intent(in) :: n ! Order of the problem type(var_stuff), intent(inout) :: h(2**(n-1)) ! Array of variable mnemonics type(scaled_odd_tail_stuff), intent(in) :: scaled_odd_tail_data(n) ! Indices ! of data for scaled_odd_tail type(correct_stuff) correct_data ! Indices for degenerate problems integer, intent(in) :: iunit ! Write code to this unit type(var_stuff) hc(2**(n-3)), hcx(2**(n-3)) ! Temporary variable mnemonics type(var_stuff) hs(2**(n-3)), hsx(2**(n-3)) ! Temporary variable mnemonics type(var_stuff) h2(2**(n-2)) ! Temporary variable mnemonics type(var_stuff) hc1(2**(n-4)), hc1x(2**(n-4)) ! Temporary variable mnemonics type(var_stuff) hs1(2**(n-4)), hs1x(2**(n-4)) ! Temporary variable mnemonics type(var_stuff) h21(2**(n-3)), h21x(2**(n-3)) ! Temporary variable mnemonics integer i ! Row number integer k ! Blocksize integer M ! Blocksize integer max_depth ! Previous maximum depth if(n <= 3) then call odd_tail_correct(h,n,correct_data,iunit) 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 call WX(' '//write_var(hc(i))//' = '// & write_var(h(4*i-2))//'-'//write_var(h(4*(M+1-i))),iunit) end do ! hc = h(2::4)-h(2**(n-1):1:-4) call scaled_odd_tail(hc,n-2,correct_data,scaled_odd_tail_data,iunit) 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 call WX(' '//write_var(hcx(i))//' = '// & write_var(hc(i))//'-'// & write_const(scaled_odd_tail_data(n)%matrixc(i))//'*'// & write_var(hc(M+i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hcx(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrixc(i))//'*'// & write_var(hc(2*M+1-i))//'+'//write_var(hc(3*M+1-i)),iunit) end do ! hc = (/(/(hc(i)-scaled_odd_tail_data(n)%matrixc(i)*hc(M+i),i=1,M)/), & ! (/(scaled_odd_tail_data(n)%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 call WX(' '//write_var(hs(i))//' = '// & write_var(h(4*i-2))//'+'//write_var(h(4*(M+1-i))),iunit) end do ! hs = h(2::4)+h(2**(n-1):1:-4) call scaled_odd_tail(hs,n-2,correct_data,scaled_odd_tail_data,iunit) 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 call WX(' '//write_var(hsx(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrixs(i))//'*'// & write_var(hs(i))//'+'//write_var(hs(M+i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hsx(i))//' = '// & write_var(hs(2*M+1-i))//'-'// & write_const(scaled_odd_tail_data(n)%matrixs(i))//'*'// & write_var(hs(3*M+1-i)),iunit) end do ! hs = (/(/(scaled_odd_tail_data(n)%matrixs(i)*hs(i)+hs(M+i),i=1,M)/), & ! (/(hs(2*M+1-i)-scaled_odd_tail_data(n)%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 call WX(' '//write_var(hc1(i))//' = '// & write_var(h2(4*i-2))//'-'//write_var(h2(4*(M+1-i))),iunit) end do ! hc1 = h2(2::4)-h2(2**(n-2):1:-4) call scaled_odd_tail(hc1,n-3,correct_data,scaled_odd_tail_data,iunit) 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 call WX(' '//write_var(hc1x(1))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1c(1,1))//'*'// & write_var(hc1(1)),iunit) else do i = 1, M call WX(' '//write_var(hc1x(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1c(1,i))//'*'// & write_var(hc1(i))//'-'// & write_const(scaled_odd_tail_data(n)%matrix1c(2,i))// & '*'//write_var(hc1(M+i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hc1x(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1c(1,i))//'*'// & write_var(hc1(2*M+1-i))//'+'// & write_const(scaled_odd_tail_data(n)%matrix1c(2,i))//'*'// & write_var(hc1(3*M+1-i)),iunit) end do end if ! if(n == 4) then ! hc1(1) = scaled_odd_tail_data(n)%matrix1c(1,1)*hc1(1) ! else ! hc1 = (/(/(scaled_odd_tail_data(n)%matrix1c(1,i)*hc1(i)- & ! scaled_odd_tail_data(n)%matrix1c(2,i)*hc1(M+i),i=1,M)/), & ! (/(scaled_odd_tail_data(n)%matrix1c(1,i)*hc1(2*M+1-i)+ & ! scaled_odd_tail_data(n)%matrix1c(2,i)*hc1(3*M+1-i),i=M+1,2*M)/)/) ! end if M = 2**(n-4) do i = 1, M call WX(' '//write_var(hs1(i))//' = '// & write_var(h2(4*i-2))//'+'//write_var(h2(4*(M+1-i))),iunit) end do ! hs1 = h2(2::4)+h2(2**(n-2):1:-4) call scaled_odd_tail(hs1,n-3,correct_data,scaled_odd_tail_data,iunit) 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 call WX(' '//write_var(hs1x(1))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1s(1,1))//'*'// & write_var(hs1(1)),iunit) else do i = 1, M call WX(' '//write_var(hs1x(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1s(1,i))//'*'// & write_var(hs1(i))//'+'// & write_const(scaled_odd_tail_data(n)%matrix1s(2,i))//'*'// & write_var(hs1(M+i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hs1x(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1s(1,i))//'*'// & write_var(hs1(2*M+1-i))//'-'// & write_const(scaled_odd_tail_data(n)%matrix1s(2,i))//'*'// & write_var(hs1(3*M+1-i)),iunit) end do end if ! if(n == 4) then ! hs1(1) = scaled_odd_tail_data(n)%matrix1s(1,1)*hs1(1) ! else ! hs1 = (/(/(scaled_odd_tail_data(n)%matrix1s(1,i)*hs1(i)+ & ! scaled_odd_tail_data(n)%matrix1s(2,i)*hs1(M+i),i=1,M)/), & ! (/(scaled_odd_tail_data(n)%matrix1s(1,i)*hs1(2*M+1-i)- & ! scaled_odd_tail_data(n)%matrix1s(2,i)*hs1(3*M+1-i),i=M+1,2*M)/)/) ! end if ! h21 = h2(::2) call scaled_odd_tail(h21,n-2,correct_data,scaled_odd_tail_data,iunit) 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 call WX(' '//write_var(h21x(i))//' = '// & write_const(scaled_odd_tail_data(n)%diag(i))//'*'// & write_var(h21(i)),iunit) end do ! h21 = scaled_odd_tail_data(n)%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 call WX(' '//write_var(h2(i))//' = '// & write_var(h21x(i))//'+'//write_var(hc1x(i)),iunit) end do ! h2(0*k+1:1*k) = h21(:k)+hc1 do i = 1, k call WX(' '//write_var(h2(k+i))//' = '// & write_var(h21x(k+1-i))//'-'//write_var(hc1x(k+1-i)),iunit) end do ! h2(1*k+1:2*k) = h21(k:1:-1)-hc1(k:1:-1) do i = 1, k call WX(' '//write_var(h2(2*k+i))//' = '// & write_var(hs1x(i))//'+'//write_var(h21x(k+i)),iunit) end do ! h2(2*k+1:3*k) = hs1+h21(k+1:) do i = 1, k call WX(' '//write_var(h2(3*k+i))//' = '// & write_var(hs1x(k+1-i))//'-'//write_var(h21x(2*k+1-i)),iunit) 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 call WX(' '//write_var(h(i))//' = '// & write_var(h2(i))//'+'//write_var(hcx(i)),iunit) end do ! h(0*k+1:1*k) = h2(:k)+hc do i = 1, k call WX(' '//write_var(h(k+i))//' = '// & write_var(h2(k+1-i))//'-'//write_var(hcx(k+1-i)),iunit) end do ! h(1*k+1:2*k) = h2(k:1:-1)-hc(k:1:-1) do i = 1, k call WX(' '//write_var(h(2*k+i))//' = '// & write_var(hsx(i))//'+'//write_var(h2(k+i)),iunit) end do ! h(2*k+1:3*k) = hs+h2(k+1:) do i = 1, k call WX(' '//write_var(h(3*k+i))//' = '// & write_var(hsx(k+1-i))//'-'//write_var(h2(2*k+1-i)),iunit) 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(h,n,iunit) ! Calculates real-half complex dft of order 2**n from definition ! for degenerate cases and immortalizes the result in code !................................................................ subroutine real_fft_correct(h,n,iunit) integer, intent(in) :: n ! Order of the problem type(var_stuff), intent(inout) :: h(0:2**n-1) ! Array of variable mnemonics integer, intent(in) :: iunit ! Write code to this unit type(var_stuff) h1(0:2**n-1) ! Temporary variable mnemonics type(var_stuff) h2(0:2**n-1) ! Temporary variable mnemonics integer i ! Row index integer max_depth ! Depth of previous problem ! 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 call WX(' '//write_var(h1(0))//' = '// & write_input(h(0))//'+'//write_input(h(1)),iunit) call WX(' '//write_var(h1(1))//' = '// & write_input(h(0))//'-'//write_input(h(1)),iunit) else call WX(' '//write_var(h1(0))//' = '// & write_var(h(0))//'+'//write_var(h(1)),iunit) call WX(' '//write_var(h1(1))//' = '// & write_var(h(0))//'-'//write_var(h(1)),iunit) end if do i = 0, 2**n-1 h(i) = h1(i) end do case(2) 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 call WX(' '//write_var(h1(0))//' = '// & write_input(h(0))//'+'//write_input(h(2)),iunit) call WX(' '//write_var(h1(2))//' = '// & write_input(h(1))//'+'//write_input(h(3)),iunit) call WX(' '//write_var(h2(0))//' = '// & write_var(h1(0))//'+'//write_var(h1(2)),iunit) call WX(' '//write_var(h2(1))//' = '// & write_input(h(0))//'-'//write_input(h(2)),iunit) call WX(' '//write_var(h2(2))//' = '// & write_var(h1(0))//'-'//write_var(h1(2)),iunit) call WX(' '//write_var(h2(3))//' = '// & write_input(h(1))//'-'//write_input(h(3)),iunit) else call WX(' '//write_var(h1(0))//' = '// & write_var(h(0))//'+'//write_var(h(2)),iunit) call WX(' '//write_var(h1(2))//' = '// & write_var(h(1))//'+'//write_var(h(3)),iunit) call WX(' '//write_var(h2(0))//' = '// & write_var(h1(0))//'+'//write_var(h1(2)),iunit) call WX(' '//write_var(h2(1))//' = '// & write_var(h(0))//'-'//write_var(h(2)),iunit) call WX(' '//write_var(h2(2))//' = '// & write_var(h1(0))//'-'//write_var(h1(2)),iunit) call WX(' '//write_var(h2(3))//' = '// & write_var(h(1))//'-'//write_var(h(3)),iunit) end if do i = 0, 2**n-1 h(i) = h2(i) end do case default stop 'bigger n than 2 in real_fft_correct' end select end subroutine real_fft_correct !................................................................ !................................................................ ! Subroutine real_fft(h,n,correct_data,odd_tail_data, & ! scaled_odd_tail_data,iunit) ! Calculates real-half complex fft of order 2**n and ! immortalizes the result in code !................................................................ recursive subroutine real_fft(h,n,correct_data,odd_tail_data, & scaled_odd_tail_data,iunit) integer, intent(in) :: n ! Order of the problem type(var_stuff), intent(inout) :: h(0:2**n-1) ! Array of variable mnemonics type(correct_stuff) correct_data ! Indices for degenerate problems type(odd_tail_stuff) odd_tail_data(n) ! Indices for odd_tail type(scaled_odd_tail_stuff) scaled_odd_tail_data(n-2) ! Indices for scaled_odd_tail integer, intent(in) :: iunit ! Write code to this unit type(var_stuff) h0(0:2**(n-1)-1) ! even freq. mnemonics type(var_stuff) h1(0:2**(n-1)-1) ! odd freq. mnemonics integer max_depth ! Max depth of previous integer M ! Blocksize integer i ! Current index if(n < 3) then call real_fft_correct(h,n,iunit) 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) call WX(' '//write_var(h0(i))//' = '// & write_input(h(i))//'+'//write_input(h(i+M)),iunit) end do ! h0 = h(0:2**(n-1)-1)+h(2**(n-1):) call real_fft(h0,n-1,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) do i = 0, 2**(n-1)-1 h1(i) = var_stuff(max_depth+1, h(i+M)%index, h(i+M)%label) call WX(' '//write_var(h1(i))//' = '// & write_input(h(i))//'-'//write_input(h(i+M)),iunit) 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) call WX(' '//write_var(h0(i))//' = '// & write_var(h(i))//'+'//write_var(h(i+M)),iunit) end do ! h0 = h(0:2**(n-1)-1)+h(2**(n-1):) call real_fft(h0,n-1,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) do i = 0, 2**(n-1)-1 h1(i) = var_stuff(max_depth+1, h(i+M)%index, h(i+M)%label) call WX(' '//write_var(h1(i))//' = '// & write_var(h(i))//'-'//write_var(h(i+M)),iunit) end do ! h1 = h(0:2**(n-1)-1)-h(2**(n-1):) end if call odd_tail(h1,n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) 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(n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) ! Computes full complex fft of order 2**n and immortalizes results in code ! subroutine complex_fft(h,n) !................................................................ subroutine complex_fft(n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) integer, intent(in) :: n ! Order of the problem type(correct_stuff) correct_data ! Indices for degenerate problems type(odd_tail_stuff) odd_tail_data(n) ! Indices for odd tail type(scaled_odd_tail_stuff) scaled_odd_tail_data(n-2) ! Indices for ! scaled form of odd tail integer, intent(in) :: iunit ! Write code to this unit type(var_stuff) hreal(0:2**n-1) ! Array of real variable mnemonics type(var_stuff) himag(0:2**n-1) ! Array of imaginary variable mnemonics integer M ! Blocksize integer i ! Row number 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,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) ! himag = aimag(h) do i = 0, 2**n-1 himag(i) = var_stuff(0,i,'xi') end do call real_fft(himag,n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) ! h(0) = cmplx(hreal(0),himag(0),dp) call WX(' h(0) = cmplx('// & write_var(hreal(0))//','//write_var(himag(0))//',wp)',iunit) do i = 1, (M-1)/2 call WX(' '//write_output(i)//' = cmplx('// & write_var(hreal(i))//'+'//write_var(himag(M/2+i))//','// & write_var(himag(i))//'-'//write_var(hreal(M/2+i))//',wp)',iunit) ! h(i) = cmplx(hreal(i)-himag(M/2+i),hreal(M/2+i)+himag(i),dp) call WX(' '//write_output(M-i)//' = cmplx('// & write_var(hreal(i))//'-'//write_var(himag(M/2+i))//','// & write_var(hreal(M/2+i))//'+'//write_var(himag(i))//',wp)',iunit) ! h(M-i) = cmplx(hreal(i)+himag(M/2+i),-hreal(M/2+i)+himag(i),dp) end do call WX(' '//write_output(M/2)//' = cmplx('// & write_var(hreal(M/2))//','//write_var(himag(M/2))//',wp)',iunit) ! h(M/2) = cmplx(hreal(M/2),himag(M/2),dp) end subroutine complex_fft !................................................................ !................................................................ ! function write_var(var) ! Writes out variable name corresponding to mnemonic !................................................................ function write_var(var) type(var_stuff), intent(in) :: var ! Variable mnemonic character(write_var_len(var)) write_var ! Printed variable name write(write_var,'(a,i0,a,i0)') trim(adjustl(var%label)), & var%depth, '_', var%index end function write_var !................................................................ !................................................................ ! function write_var_len(var) ! Specification function to compute len(write_var) !................................................................ pure function write_var_len(var) integer write_var_len ! Len(write_var(var)) type(var_stuff), intent(in) :: var ! Variable mnemonic character(80) test_string ! Holds name temporarily 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) ! Write initial variable name in h array !................................................................ function write_input(var) type(var_stuff), intent(in) :: var ! Variable mnemonic character(write_input_len(var)) write_input ! Variable name 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 !................................................................ !................................................................ ! function write_input_len(var) ! Specification function to compute len(write_input(var)) !................................................................ pure function write_input_len(var) type(var_stuff), intent(in) :: var ! Variable mnemonic integer write_input_len ! Len(write_input(var)) character(80) test_string ! Holds name temporarily 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 !................................................................ !................................................................ ! function write_const(var) ! Write constant name !................................................................ function write_const(var) integer, intent(in) :: var ! Constant index character(write_const_len(var)) write_const ! Constant name write(write_const,'(a,i0)') 'd',var end function write_const !................................................................ !................................................................ ! function write_const_len(var) ! Specification function to compute len(write_const(var)) !................................................................ pure function write_const_len(var) integer, intent(in) :: var ! Constant index integer write_const_len ! Len(write_const(var)) character(80) test_string ! Holds name temporarily write(test_string,'(a,i0)') 'd',var write_const_len = len_trim(test_string) end function write_const_len !................................................................ !................................................................ ! function write_output(var) ! Write array element name !................................................................ function write_output(var) integer, intent(in) :: var ! Array element index ! The following line causes CVF 6.6C to crash at runtime. The ! workaround, commented out below, was found by James Giles. character(write_output_len(var)) write_output ! Array element name ! character(write_output_len((var))) write_output ! Array element name write(write_output,'(a,i0,a)') 'h(',var,')' end function write_output !................................................................ !................................................................ ! function write_output_len(var) ! Specification function to compute len(write_output(var)) !................................................................ pure function write_output_len(var) integer, intent(in) :: var ! Array element index integer write_output_len ! Len(write_output(var)) character(80) test_string ! Holds name temporarily write(test_string,'(a,i0,a)') 'h(',var,')' write_output_len = len_trim(test_string) end function write_output_len !................................................................ !................................................................ ! subroutine WX(string,iunit) ! Writes string to iunit !................................................................ subroutine WX(string,iunit) integer, intent(in) :: iunit character(*), intent(in) :: string write(iunit,'(a)') string end subroutine WX !................................................................ end module code_gen ! *************************************************************** ! *************************************************************** ! program test ! Prompts user for size of the problem and creates a dft subroutine ! of order 2**n along with a test program. ! *************************************************************** program test use mykinds use code_gen use avl_mod implicit none integer n ! Order of the problem integer i ! Subproblem order real(dp) t0, t1 ! Timer start and stop times real(qp), allocatable :: qconst_array(:) ! Array of required constants character(80) filename1 ! Name of output *.f90 file character(80) filename2 ! Name of output *.i90 file character(8) date ! Date for DATE_AND_TIME character(10) time ! Time for DATE_AND_TIME integer An, Mn ! Theoretical opcounts type(get_factor_stuff), allocatable :: get_factor_data(:) ! Scale ! factors from scaled_odd_tail type(correct_stuff) correct_data ! Indices for degenerate problems type(odd_tail_stuff), allocatable :: odd_tail_data(:) ! Indices for ! unscaled odd_tail problem type(scaled_odd_tail_stuff), allocatable :: scaled_odd_tail_data(:) ! Indices ! for scaled odd tail problem type(avl_bag) body ! Holds actual data integer iunit ! write code to this unit integer junit ! test code unit character(80) temp_out ! Staging area for numeric output write(*,'(a)',advance='no') ' Enter the exponent of the problem:> ' read(*,*) n if(n <= 0) then write(*,'(a)') ' Sorry, the exponent must be positive.' stop else if(n > 24) then write(*,'(a)') ' Sorry, the maximum value is 2**24.' write(*,'(a,i0,a)') ' Your input corresponds to 2**',n,'.' stop end if write(*,*) ' Creating DFT code:' call cpu_time(t0) ! Write test program write(filename1,'(a,i0,a)') 'fft', 2**n, '.f90' junit = 11 open(junit,file=filename1,status='replace') write(junit,'(a,a)') '! File: ', trim(filename1) call date_and_time(date, time) write(junit,'(12a)') '! Generated by scaled2h.f90 ', date(5:6), & '/', date(7:8), '/', date(1:4), ' ', time(1:2), & ':', time(3:4), ':', time(5:10) write(junit,'(a)') 'module mykinds' write(junit,'(a)') ' implicit none' write(junit,'(a)') ' integer, parameter :: sp = selected_real_kind(6,30)' write(junit,'(a)') ' integer, parameter :: dp = selected_real_kind(15,300)' write(junit,'(a)') 'end module mykinds' write(junit,'(a)') '' write(junit,'(a)') 'module dp_mod' write(junit,'(a)') ' use mykinds, only: wp => dp' write(junit,'(a)') ' implicit none' write(junit,'(a)') ' private wp' write(junit,'(a)') ' contains' write(junit,'(a,i0,a)') 'include ''fft',2**n,'.i90''' write(junit,'(a)') 'end module dp_mod' write(junit,'(a)') '' write(junit,'(a)') 'program test' write(junit,'(a)') ' use mykinds' write(junit,'(a)') ' use dp_mod' write(junit,'(a)') ' implicit none' write(junit,'(a,i0)') ' integer, parameter :: N = ',2**n write(junit,'(a)') ' complex(dp) h1(0:N-1)' write(junit,'(a)') ' complex(dp) h2(0:N-1)' write(junit,'(a)') ' complex(dp), parameter :: ri(2) = (/(1,0),(0,1)/)' write(junit,'(a)') ' integer i, j, k, L' write(junit,'(a)') ' real(dp) pi' write(junit,'(a)') ' real(dp) err, maxerr' write(junit,'(a)') '' write(junit,'(a)') ' pi = 4*atan(1.0_dp)' write(junit,'(a)') ' maxerr = 0' write(junit,'(a)') ' do i = 0, N-1' write(junit,'(a)') ' do j = 1, 2' write(junit,'(a)') ' h2 = 0' write(junit,'(a)') ' h2(i) = ri(j)' write(junit,'(a)') ' h1 = (/(h2(i)*exp(-2*pi*(0,1)*mod(i*k,N)/N),k=0,N-1)/)' write(junit,'(a,i0,a)') ' call fft',2**n,'(h2)' write(junit,'(a)') ' err = max(maxval(abs(real(h1-h2))),maxval(abs(aimag(h1-h2))))' write(junit,'(a)') ' write(*,*) i, j, err' write(junit,'(a)') ' maxerr = max(err, maxerr)' write(junit,'(a)') ' end do' write(junit,'(a)') ' end do' write(junit,'(a)') ' write(*,*)' write(junit,'(a)') " write(*,*) 'maxerr = ',maxerr" write(junit,'(a)') 'end program test' close(junit) ! Write DFT subroutine file write(filename2,'(a,i0,a)') 'fft', 2**n, '.i90' iunit = 10 open(iunit,file=filename2,status='replace') write(iunit,'(a,a)') '! File: ', trim(filename2) call date_and_time(date, time) write(iunit,'(12a)') '! Generated by scaled2h.f90 ', date(5:6), & '/', date(7:8), '/', date(1:4), ' ', time(1:2), & ':', time(3:4), ':', time(5:10) write(iunit,'(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(iunit,'(a,i0,a,i0,a)') '! ',An,' adds, ',Mn,' muls' write(iunit,'(a)') ' implicit real(wp) (x,y)' write(iunit,'(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, body, correct_data) allocate(odd_tail_data(n)) do i = 1, n call odd_tail_gen(odd_tail_data(i),i,body,get_factor_data) 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,body,get_factor_data) end do allocate(qconst_array(body%num_nodes)) call bag2array(body, qconst_array) write(*,'()') call cpu_time(t1) do i = 1, size(qconst_array) write(temp_out,'(g41.34)') qconst_array(i) write(iunit,'(a,i0,a,a,a)') ' real(wp), parameter :: d', & i, ' = ', trim(adjustl(temp_out)), '_wp' end do deallocate(qconst_array) write(*,'(a,f0.3)') ' Prep time = ', t1-t0 call cpu_time(t0) call complex_fft(n,correct_data,odd_tail_data,scaled_odd_tail_data,iunit) call cpu_time(t1) write(*,'(a,f0.3)') ' Execution time = ', t1-t0 write(iunit,'(a,i0)') 'end subroutine fft', 2**n write(*,'(a)') ' Files created: '//trim(filename1)//', '//trim(filename2) end program test ! ***************************************************************