! File: conv2c.f90 ! Public domain 2004 James Van Buskirk ! PURPOSE: Generates fast convolution algorithms of order 2**n. ! Prompts the user for the exponent of the problem, and if the ! user wants to generate fast convolution code for problems of ! order 2**n, he enters n at the prompt. Three files are then ! created: ! 1) const###.i90 contains a subroutine const###(h,a) that ! takes a vector h of order ### = 2**n and prepares an array ! of constants of order Na = max(2**n,3*2**n/2-2) necessary ! to convolve another vector by h. ! 2) eval###.i90 contains a subroutine eval###(h,a) that ! converts a vector h of order ### = 2**n into its convolution ! with the vector from which the vector a of order Na = ! max(2**n,3*2**n-2) was created. ! 3) conv###.f90 contains a test program that shows how one has ! to INCLUDE the above two files in a module that defines ! the parameter wp to be a kind type parameter suitable for ! real types, and the sequence of calls, first one to const### ! to get an a vector, then one or more calls to eval### to ! convolve other vectors with the first. ! If everything works OK, the conv###.f90 file should compile ! to a program that prints out 2**n lines of output that are ! near EPSILON(1.0d0) and a final line of output that is the ! maxmimum of the previous 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 real convolution 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 factor_stuff ! Contains pointers to block factors !................................................................ type factor_stuff integer, pointer :: factor(:) ! Pointer to array of indices end type factor_stuff !................................................................ !................................................................ ! type rot_mat_entry ! Contains indices of rotation matrix entries !................................................................ type rot_mat_entry integer cos integer sin end type rot_mat_entry !................................................................ !................................................................ ! type rot_mat_stuff ! Contains pointer to rotation matrix entry indices !................................................................ type rot_mat_stuff type(rot_mat_entry), pointer :: rot_mats(:) ! Pointer to array of rotation matrix indices end type rot_mat_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(4). 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(2**n) !................................................................ recursive subroutine get_factor(factor, n, old_factor) integer, intent(in) :: n ! Order of the problem = 2**n real(qp), intent(out) :: factor(2**(n-3)) ! Scale factor real(qp), intent(in) :: old_factor(2**(n-5)) ! 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 <= 2) then factor = 1 else factor = (/(cos(2*pi*i/2**(n+1)),i=1,2**(n-2),2)/) if(n > 4) then k = 2**(n-5) 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 = 2**n type(correct_stuff) correct_data ! Holds indices of keys created real(qp) constant ! Current constant created real(qp) pi ! Pi if(n > 3) 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 prepare_factors(n,body,factor_data,get_factor_data) ! Prepares constant block factors !................................................................ subroutine prepare_factors(n, body, factor_data,get_factor_data) type(avl_bag), intent(inout) :: body ! Holds any keys created integer, intent(in) :: n ! Order of overall problem = 2**n type(factor_stuff), intent(out) :: factor_data(n-1) ! Holds indices of keys created type(get_factor_stuff), intent(in) :: get_factor_data(0:n-2) real(qp) constant ! Current constant created real(qp) four_2_n ! 4/2**n integer k ! Order of current negacyclic block = 2**k integer j ! Factor index four_2_n = 4/real(2,qp)**n k = 1 if(n >= 2) then allocate(factor_data(1)%factor(1)) constant = four_2_n/2 call key2index(constant, body, factor_data(1)%factor(1)) end if do k = 2, min(3,n-1) allocate(factor_data(k)%factor(1)) constant = four_2_n call key2index(constant, body, factor_data(k)%factor(1)) end do do k = 4, n-1 allocate(factor_data(k)%factor(2**(k-4))) do j = 1, 2**(k-4) constant = four_2_n*get_factor_data(k-1)%factor(j)**3 call key2index(constant, body, factor_data(k)%factor(j)) end do end do end subroutine prepare_factors !................................................................ !................................................................ ! Subroutine prepare_rot_mat(n,body,rot_mat_data) ! Prepares constants block factors !................................................................ Subroutine prepare_rot_mat(n,body,rot_mat_data) type(avl_bag), intent(inout) :: body ! Holds any keys created integer, intent(in) :: n ! Order of overall problem = 2**n type(rot_mat_stuff), intent(out) :: rot_mat_data(:) ! Holds indices of keys created integer k ! Order of current negacyclic block = 2**k integer j ! Current block index integer m ! Order of current negacyclic block real(qp) pi ! pi real(qp) constant ! Current constant craeted pi = 4*atan(1.0_qp) do k = 1, min(n-1,2) allocate(rot_mat_data(k)%rot_mats(0)) end do do k = 3, n-1 allocate(rot_mat_data(k)%rot_mats(2**(k-3))) m = 2**k do j = 1, m/8 constant = cos(pi*(2*j-1)/(m/2)) call key2index(constant, body, rot_mat_data(k)%rot_mats(j)%cos) constant = sin(pi*(2*j-1)/(m/2)) call key2index(constant, body, rot_mat_data(k)%rot_mats(j)%sin) end do end do end subroutine prepare_rot_mat !................................................................ !................................................................ ! Subroutine odd_tail_correct(h,n,correct_data,iunit) ! Computes odd_tail(2**n)*h using the theoretical formula for ! orders too degenerate for recursion and immortalizes the ! result in code !................................................................ subroutine odd_tail_correct(h,n,correct_data,iunit) integer, intent(in) :: n ! Order of the problem = 2**n type(var_stuff), intent(inout) :: h(2**n) ! 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), h2(2**n) ! Temporary variable mnemonics integer max_depth ! Previous maximum depth integer i ! General index select case(n) case(0) case(1) case(2) 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 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 h(i) = var_stuff(max_depth+2,h2(i)%index,h2(i)%label) end do case default stop ' bigger n than 2 in odd_tail_correct' end select end subroutine odd_tail_correct !................................................................ !................................................................ ! Subroutine odd_tail_correct_transpose(h,n,correct_data,iunit) ! Computes transpose of odd_tail(2**n)*h using the theoretical formula for ! orders too degenerate for recursion and immortalizes the result ! in code !................................................................ subroutine odd_tail_correct_transpose(h,n,correct_data,iunit) integer, intent(in) :: n ! Order of the problem = 2**n type(var_stuff), intent(inout) :: h(2**n) ! 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), h2(2**n) ! Temporary variable mnemonics integer max_depth ! Previous maximum depth integer i ! General index select case(n) case(0) case(1) case(2) 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(1))//'-'//write_var(h(2))//')',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(3))//'+'//write_var(h(4))//')',iunit) ! xs = const_array(correct_data%rsq2)*(h(2)+h(4)) do i = 1, 2**n 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(h(2)),iunit) call WX(' '//write_var(h2(2))//' = '// & write_var(h1(2))//'+'//write_var(h1(1)),iunit) call WX(' '//write_var(h2(3))//' = '// & write_var(h(3))//'-'//write_var(h(4)),iunit) call WX(' '//write_var(h2(4))//' = '// & write_var(h1(2))//'-'//write_var(h1(1)),iunit) ! h = (/h(1)+xc,h(1)-xc,xs+h(3),xs-h(3)/) do i = 1, 2**n h(i) = var_stuff(max_depth+2,h2(i)%index,h2(i)%label) end do case default stop ' bigger n than 2 in odd_tail_correct_transpose' end select end subroutine odd_tail_correct_transpose !................................................................ !................................................................ ! 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 = 2**n 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(0:n) ! Factors ! thrown out by scaled_odd_tail real(qp) factor(2**(n-3)) ! Factors passed through integer i ! Row number integer j ! Column select integer k ! Blocksize real(qp) pi ! Pi real(qp) matrix(2**(n-2)) ! Actual array of coefficients real(qp) matrix1(min(2,2**(n-3)),2**(n-3)) ! Actual matrix of coefficients real(qp) diag(2**(n-2)) ! Actual array of coefficients pi = 4*atan(1.0_qp) allocate(stuff%matrixc(2**(n-2))) allocate(stuff%matrixs(2**(n-2))) allocate(stuff%matrix1c(min(2,2**(n-3)),2**(n-3))) allocate(stuff%matrix1s(min(2,2**(n-3)),2**(n-3))) allocate(stuff%diag(2**(n-2))) if(n == 2) then ! Insert garbage values: should never be used stuff%matrixc = 0 stuff%matrixs = 0 stuff%diag = 0 else if(n > 2) then factor = get_factor_data(n)%factor k = 2**(n-3) do i = 1, k matrix(2*k+1-i) = tan(2*pi*(2*i-1)/2**(n+1)) matrix(i) = -tan(2*pi*(2*i-1)/2**(n+1)) end do do i = 1, 2**(n-2) call key2index(abs(matrix(i)), body, stuff%matrixc(i)) end do k = 2**(n-3) do i = 1, k matrix(i) = tan(2*pi*(2*i-1)/2**(n+1)) matrix(2*k+1-i) = -tan(2*pi*(2*i-1)/2**(n+1)) end do do i = 1, 2**(n-2) call key2index(abs(matrix(i)), body, stuff%matrixs(i)) end do k = 2**(n-4) ! Next line just in case n == 3 matrix1(1,1) = cos(pi/4) ! General case do i = 1, k matrix1(1,i) = cos(2*pi*(2*i-1)/2**n) matrix1(1,2*k+1-i) = sin(2*pi*(2*i-1)/2**n) matrix1(2,i) = -sin(2*pi*(2*i-1)/2**n) matrix1(2,2*k+1-i) = cos(2*pi*(2*i-1)/2**n) end do k = 2**(n-6) 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-3) do i = 1, k matrix1(:,0*k+i) = matrix1(:,0*k+i)/factor(i) end do do i = 1, 2**(n-3) do j = 1, min(2,2**(n-3)) call key2index(abs(matrix1(j,i)), body, stuff%matrix1c(j,i)) end do end do k = 2**(n-4) ! Next line just in case n == 3 matrix1(1,1) = cos(pi/4) ! General case do i = 1, k matrix1(1,i) = sin(2*pi*(2*i-1)/2**n) matrix1(1,2*k+1-i) = cos(2*pi*(2*i-1)/2**n) matrix1(2,i) = cos(2*pi*(2*i-1)/2**n) matrix1(2,2*k+1-i) = -sin(2*pi*(2*i-1)/2**n) end do k = 2**(n-6) 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-3) do i = 1, k matrix1(:,0*k+i) = matrix1(:,0*k+i)/factor(i) end do do i = 1, 2**(n-3) do j = 1, min(2,2**(n-3)) call key2index(abs(matrix1(j,i)), body, stuff%matrix1s(j,i)) end do end do k = 2**(n-5) 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-3) 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-2) 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 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 = 2**n type(var_stuff), intent(inout) :: h(2**n) ! Array of variable mnemonics type(scaled_odd_tail_stuff), intent(in) :: scaled_odd_tail_data(0: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-2)), hcx(2**(n-2)) ! Temporary variable mnemonics type(var_stuff) hs(2**(n-2)), hsx(2**(n-2)) ! Temporary variable mnemonics type(var_stuff) h2(2**(n-1)) ! Temporary variable mnemonics type(var_stuff) hc1(2**(n-3)), hc1x(2**(n-3)) ! Temporary variable mnemonics type(var_stuff) hs1(2**(n-3)), hs1x(2**(n-3)) ! Temporary variable mnemonics type(var_stuff) h21(2**(n-2)), h21x(2**(n-2)) ! Temporary variable mnemonics integer i ! Row number integer k ! Blocksize integer M ! Blocksize integer max_depth ! Previous maximum depth if(n <= 2) then call odd_tail_correct(h,n,correct_data,iunit) else max_depth = maxval(h%depth) M = 2**(n-2) 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:-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-3) 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-2) 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:-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-3) 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-3) 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-1):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-4) if(n == 3) 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 == 3) 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-3) 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-1):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-4) if(n == 3) 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 == 3) 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-2) 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-3) 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-2) 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 scaled_odd_tail_transpose(h,n,correct_data,scaled_odd_tail_data,iunit) ! Computes action of transpose of scaled odd tail matrix of order 2**n and ! immortalizes in code. !................................................................ recursive subroutine scaled_odd_tail_transpose(h,n,correct_data, & scaled_odd_tail_data,iunit) integer, intent(in) :: n ! Order of the problem = 2**n type(var_stuff), intent(inout) :: h(2**n) ! Array of variable mnemonics type(scaled_odd_tail_stuff), intent(in) :: scaled_odd_tail_data(0: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-2)), hcx(2**(n-2)) ! Temporary variable mnemonics type(var_stuff) hs(2**(n-2)), hsx(2**(n-2)) ! Temporary variable mnemonics type(var_stuff) h2(2**(n-1)) ! Temporary variable mnemonics type(var_stuff) hc1(2**(n-3)), hc1x(2**(n-3)) ! Temporary variable mnemonics type(var_stuff) hs1(2**(n-3)), hs1x(2**(n-3)) ! Temporary variable mnemonics type(var_stuff) h21(2**(n-2)), h21x(2**(n-2)) ! Temporary variable mnemonics integer i ! Row number integer k ! Blocksize integer M ! Blocksize integer max_depth ! Previous maximum depth if(n <= 2) then call odd_tail_correct_transpose(h,n,correct_data,iunit) else !********************** Last RF(2) ************************* k = 2**(n-2) max_depth = maxval(h%depth) do i = 1, k h2(i) = var_stuff(max_depth+1,h(2*i-1)%index,h(2*i-1)%label) hcx(i) = var_stuff(max_depth+1,h(2*i)%index,h(2*i)%label) hsx(i) = var_stuff(max_depth+1,h(2*(i+k))%index,h(2*i+k)%label) h2(i+k) = var_stuff(max_depth+1,h(2*(i+k)-1)%index,h(2*(i+k)-1)%label) end do do i = 1, k call WX(' '//write_var(h2(i))//' = '// & write_var(h(i))//'+'//write_var(h(2*k+1-i)),iunit) end do ! h(0*k+1:1*k) = h2(:k)+hc do i = 1, k call WX(' '//write_var(hcx(i))//' = '// & write_var(h(i))//'-'//write_var(h(2*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(hsx(i))//' = '// & write_var(h(2*k+i))//'+'//write_var(h(4*k+1-i)),iunit) end do ! h(2*k+1:3*k) = hs+h2(k+1:) do i = 1, k call WX(' '//write_var(h2(k+i))//' = '// & write_var(h(2*k+i))//'-'//write_var(h(4*k+1-i)),iunit) end do ! h(3*k+1:4*k) = hs(k:1:-1)-h2(2*k:k+1:-1) !********************** Last RF(2) ************************* !******************** Penultimate RF(2) ************************ k = 2**(n-3) max_depth = maxval(h2%depth) do i = 1, k h21x(i) = var_stuff(max_depth+1,h2(2*i-1)%index,h2(2*i-1)%label) hc1x(i) = var_stuff(max_depth+1,h2(2*i)%index,h2(2*i)%label) hs1x(i) = var_stuff(max_depth+1,h2(2*(i+k))%index,h2(2*(i+k))%label) h21x(i+k) = var_stuff(max_depth+1,h2(2*(i+k)-1)%index,h2(2*(i+k)-1)%label) end do do i = 1, k call WX(' '//write_var(h21x(i))//' = '// & write_var(h2(i))//'+'//write_var(h2(2*k+1-i)),iunit) end do ! h2(0*k+1:1*k) = h21(:k)+hc1 do i = 1, k call WX(' '//write_var(hc1x(i))//' = '// & write_var(h2(i))//'-'//write_var(h2(2*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(hs1x(i))//' = '// & write_var(h2(2*k+i))//'+'//write_var(h2(4*k+1-i)),iunit) end do ! h2(2*k+1:3*k) = hs1+h21(k+1:) do i = 1, k call WX(' '//write_var(h21x(k+i))//' = '// & write_var(h2(2*k+i))//'-'//write_var(h2(4*k+1-i)),iunit) end do ! h2(3*k+1:4*k) = hs1(k:1:-1)-h21(2*k:k+1:-1) !******************** Penultimate RF(2) ************************ !********************* Diagonal product *********************** max_depth = maxval(h21x%depth) M = 2**(n-2) do i = 1, M h21(i) = var_stuff(max_depth+1,h21x(i)%index,h21x(i)%label) end do do i = 1, M call WX(' '//write_var(h21(i))//' = '// & write_const(scaled_odd_tail_data(n)%diag(i))// & '*'//write_var(h21x(i)),iunit) end do ! h21 = scaled_odd_tail_data(n)%diag*h21 !********************* Diagonal product *********************** !********************* 2-deep scaled odd tail ***************** call scaled_odd_tail_transpose(h21,n-2,correct_data, & scaled_odd_tail_data,iunit) !********************* 2-deep scaled odd tail ***************** !******************** 1-deep Sin 2X2 ********************** max_depth = maxval(hs1x%depth) M = 2**(n-3) do i = 1, M hs1(i) = var_stuff(max_depth+1,hs1x(i)%index,hs1x(i)%label) end do M = 2**(n-4) if(n == 3) then call WX(' '//write_var(hs1(1))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1s(1,1))// & '*'//write_var(hs1x(1)),iunit) else do i = 1, M call WX(' '//write_var(hs1(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1s(1,i))// & '*'//write_var(hs1x(i))//'+'// & write_const(scaled_odd_tail_data(n)%matrix1s(1,2*M+1-i))// & '*'//write_var(hs1x(2*M+1-i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hs1(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1s(2,i-M))// & '*'//write_var(hs1x(i-M))//'-'// & write_const(scaled_odd_tail_data(n)%matrix1s(2,3*M+1-i))// & '*'//write_var(hs1x(3*M+1-i)),iunit) end do end if ! if(n == 3) 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) !******************** 1-deep Sin 2X2 ********************** !******************** 1-deep SOT Sin ********************** call scaled_odd_tail_transpose(hs1,n-3,correct_data, & scaled_odd_tail_data,iunit) !******************** 1-deep SOT Sin ********************** !******************** 1-deep Cos 2X2 ********************** max_depth = maxval(hc1x%depth) M = 2**(n-3) do i = 1, M hc1(i) = var_stuff(max_depth+1,hc1x(i)%index,hc1x(i)%label) end do M = 2**(n-4) if(n == 3) then call WX(' '//write_var(hc1(1))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1c(1,1))// & '*'//write_var(hc1x(1)),iunit) else do i = 1, M call WX(' '//write_var(hc1(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1c(1,i))// & '*'//write_var(hc1x(i))//'+'// & write_const(scaled_odd_tail_data(n)%matrix1c(1,2*M+1-i))// & '*'//write_var(hc1x(2*M+1-i)), iunit) end do do i = M+1,2*M call WX(' '//write_var(hc1(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrix1c(2,3*M+1-i))// & '*'//write_var(hc1x(3*M+1-i))//'-'// & write_const(scaled_odd_tail_data(n)%matrix1c(2,i-M))// & '*'//write_var(hc1x(i-M)), iunit) end do end if ! if(n == 3) 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 !******************** 1-deep Cos 2X2 ********************** !******************** 1-deep SOT Cos ********************** call scaled_odd_tail_transpose(hc1,n-3,correct_data, & scaled_odd_tail_data,iunit) !******************** 1-deep SOT Cos ********************** !******************** Second RF(2) ************************ max_depth = max(maxval(h21%depth),maxval(hc1%depth),maxval(hs1%depth)) M = 2**(n-3) do i = 1, M h2(2*i-1) = h21(i) h2(2*i) = var_stuff(max_depth+1,hc1(i)%index,hc1(i)%label) h2(2*(i+M)) = var_stuff(max_depth+1,hs1(i)%index,hs1(i)%label) h2(2*(i+M)-1) = h21(i+M) end do do i = 1, M call WX(' '//write_var(h2(4*i-2))//' = '// & write_var(hs1(i))//'+'//write_var(hc1(i)),iunit) end do ! hc1 = h2(2::4)-h2(2**(n-1):1:-4) do i = 1, M call WX(' '//write_var(h2(4*(M+1-i)))//' = '// & write_var(hs1(i))//'-'//write_var(hc1(i)),iunit) end do ! hs1 = h2(2::4)+h2(2**(n-1):1:-4) !******************** Second RF(2) ************************ !*********************** Sin 2X2 **************************** max_depth = maxval(hsx%depth) M = 2**(n-2) do i = 1, M hs(i) = var_stuff(max_depth+1,hsx(i)%index,hsx(i)%label) end do M = 2**(n-3) do i = 1, M call WX(' '//write_var(hs(i))//' = '// & write_const(scaled_odd_tail_data(n)%matrixs(i))// & '*'//write_var(hsx(i))//'+'//write_var(hsx(2*M+1-i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hs(i))//' = '// & write_var(hsx(i-M))//'-'// & write_const(scaled_odd_tail_data(n)%matrixs(3*M+1-i))// & '*'//write_var(hsx(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) !*********************** Sin 2X2 **************************** !*********************** SOT Sin **************************** call scaled_odd_tail_transpose(hs,n-2,correct_data, & scaled_odd_tail_data,iunit) !*********************** SOT Sin **************************** !*********************** Cos 2X2 **************************** max_depth = maxval(hcx%depth) M = 2**(n-2) do i = 1, M hc(i) = var_stuff(max_depth+1,hcx(i)%index,hcx(i)%label) end do M = 2**(n-3) do i = 1, M call WX(' '//write_var(hc(i))//' = '// & write_var(hcx(i))//'+'// & write_const(scaled_odd_tail_data(n)%matrixc(2*M+1-i))// & '*'//write_var(hcx(2*M+1-i)),iunit) end do do i = M+1,2*M call WX(' '//write_var(hc(i))//' = '// & write_var(hcx(3*M+1-i))//'-'// & write_const(scaled_odd_tail_data(n)%matrixc(i-M))// & '*'//write_var(hcx(i-M)),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)/)/) !*********************** Cos 2X2 **************************** !*********************** SOT Cos **************************** call scaled_odd_tail_transpose(hc,n-2,correct_data, & scaled_odd_tail_data,iunit) !*********************** SOT Cos **************************** !*********************** First RF(2) ************************ max_depth = max(maxval(h2%depth),maxval(hc%depth),maxval(hs%depth)) M = 2**(n-2) do i = 1, M h(2*i-1) = h2(i) h(2*i) = var_stuff(max_depth+1,hc(i)%index,hc(i)%label) h(2*(i+M)) = var_stuff(max_depth+1,hs(i)%index,hs(i)%label) h(2*(i+M)-1) = h2(i+M) end do do i = 1, M call WX(' '//write_var(h(4*i-2))//' = '// & write_var(hs(i))//'+'//write_var(hc(i)),iunit) end do ! hc = h(2::4)-h(2**n:1:-4) M = 2**(n-2) do i = 1, M call WX(' '//write_var(h(4*(M+1-i)))//' = '// & write_var(hs(i))//'-'//write_var(hc(i)),iunit) end do ! hs = h(2::4)+h(2**n:1:-4) !*********************** First RF(2) ************************ end if end subroutine scaled_odd_tail_transpose !................................................................ !................................................................ ! Subroutine conv_const_correct(h,n,factor_data,iunit) ! Computes constants for real convolution of order 2**n < 8 and immortalizes ! results in code !................................................................ subroutine conv_const_correct(h,n,factor_data,iunit) integer, intent(in) :: n ! Order of the problem = 2**n type(var_stuff), intent(in) :: h(0:2**n-1) ! Array of variable mnemonics type(factor_stuff) factor_data(n-1) ! Indices for block factors integer, intent(in) :: iunit ! Write code to this unit integer max_depth ! Depth of previous problem select case(n) case(0) call WX(' a(0) = h(0)',iunit) case(1) call WX(' a(0) = h(0)',iunit) call WX(' a(1) = h(1)',iunit) case(2) max_depth = maxval(h%depth) if(max_depth == 0) then call WX(' a(0) = '//write_const(factor_data(1)%factor(1))// & '*(h(0)+h(2))',iunit) call WX(' a(1) = '//write_const(factor_data(1)%factor(1))// & '*(h(1)+h(3))',iunit) call WX(' a(2) = '//write_const(factor_data(1)%factor(1))// & '*(h(0)-h(2))',iunit) call WX(' a(3) = '//write_const(factor_data(1)%factor(1))// & '*(h(1)-h(3))',iunit) else call WX(' a(0) = '//write_const(factor_data(1)%factor(1))// & '*('//write_var(h(0))//'+'//write_var(h(2))//')',iunit) call WX(' a(1) = '//write_const(factor_data(1)%factor(1))// & '*('//write_var(h(1))//'+'//write_var(h(3))//')',iunit) call WX(' a(2) = '//write_const(factor_data(1)%factor(1))// & '*('//write_var(h(0))//'-'//write_var(h(2))//')',iunit) call WX(' a(3) = '//write_const(factor_data(1)%factor(1))// & '*('//write_var(h(1))//'-'//write_var(h(3))//')',iunit) end if case default stop ' bigger n than 2 in conv_cons_correct' end select end subroutine conv_const_correct !................................................................ !................................................................ ! recursive Subroutine conv_const(h,n,correct_data,factor_data,rot_mat_data, & ! scaled_odd_tail_data,iunit) ! Computes constants for real convolution of order 2**n and immortalizes ! results in code !................................................................ recursive subroutine conv_const(h,n,correct_data,factor_data,rot_mat_data, & scaled_odd_tail_data,iunit) integer, intent(in) :: n ! Order of the problem = 2**n type(correct_stuff) correct_data ! Indices for degenerate problems type(factor_stuff) factor_data(n-1) ! Indices for block factors type(rot_mat_stuff) rot_mat_data(n-1) ! Indices for rotation matrices type(scaled_odd_tail_stuff) scaled_odd_tail_data(0:n-2) ! Indices for ! scaled form of odd tail integer, intent(in) :: iunit ! Write code to this unit type(var_stuff), intent(inout) :: h(0:2**n-1) ! Array of variable mnemonics type(var_stuff) h0(0:2**(n-1)-1) ! Array of variable mnemonics type(var_stuff) h1(0:2**(n-1)-1) ! Array of variable mnemonics integer i ! Current index integer M ! Blocksize integer j ! Index of current factor integer m0 ! Current output array index type(var_stuff) hs(4) ! Array of variable mnemonics integer max_depth ! Depth of previous problem if(n < 3) then call conv_const_correct(h,n,factor_data,iunit) else max_depth = maxval(h%depth) if(max_depth == 0) then M = 2**(n-1) do i = 0, M-1 h0(i) = var_stuff(max_depth+1, h(i)%index, h(i)%label) call WX(' '//write_var(h0(i))//' = h('// & write_number(i)//')+h('//write_number(i+M)//')',iunit) end do call conv_const(h0,n-1,correct_data,factor_data,rot_mat_data, & scaled_odd_tail_data,iunit) do i = 0, M-1 h1(i) = var_stuff(max_depth+1, h(i+M)%index, h(i+M)%label) call WX(' '//write_var(h1(i))//' = h('// & write_number(i)//')-h('//write_number(i+M)//')',iunit) end do else M = 2**(n-1) do i = 0, M-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 call conv_const(h0,n-1,correct_data,factor_data,rot_mat_data, & scaled_odd_tail_data,iunit) do i = 0, M-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 end if if(n == 3) then max_depth = maxval(h1%depth) hs(1) = var_stuff(max_depth+1, h1(0)%index, h1(0)%label) hs(2) = var_stuff(max_depth+1, h1(1)%index, h1(1)%label) hs(3) = var_stuff(max_depth+1, h1(2)%index, h1(2)%label) hs(4) = var_stuff(max_depth+1, h1(3)%index, h1(3)%label) j = 1 call WX(' '//write_var(hs(1))//' = '// & write_const(factor_data(n-1)%factor(j))//'*'// & write_var(h1(0)),iunit) call WX(' '//write_var(hs(2))//' = '// & write_const(factor_data(n-1)%factor(j))//'*'// & write_var(h1(1)),iunit) call WX(' '//write_var(hs(3))//' = '// & write_const(factor_data(n-1)%factor(j))//'*'// & write_var(h1(2)),iunit) call WX(' '//write_var(hs(4))//' = '// & write_const(factor_data(n-1)%factor(j))//'*'// & write_var(h1(3)),iunit) call WX(' a(4) = '//write_var(hs(1)),iunit) call WX(' a(5) = '//write_var(hs(2))//'-'//write_var(hs(1)),iunit) call WX(' a(6) = '//write_var(hs(4))//'+'//write_var(hs(1)),iunit) call WX(' a(7) = '//write_var(hs(3)),iunit) call WX(' a(8) = '//write_var(hs(4))//'-'//write_var(hs(3)),iunit) call WX(' a(9) = '//write_var(hs(2))//'-'//write_var(hs(3)),iunit) else call scaled_odd_tail(h1(0::2),n-2,correct_data,scaled_odd_tail_data,iunit) call scaled_odd_tail(h1(1::2),n-2,correct_data,scaled_odd_tail_data,iunit) max_depth = maxval(h1%depth) m0 = 3*M/2-2 do i = 1, M/4 j = modulo(i-1,M/8)+1 if(j > M/16) then j = M/8+1-j end if hs(1) = var_stuff(max_depth+1, h1(2*(i-1))%index, h1(2*(i-1))%label) hs(2) = var_stuff(max_depth+1, h1(2*i-1)%index, h1(2*i-1)%label) hs(3) = var_stuff(max_depth+1, h1(2*(i-1)+M/2)%index, h1(2*(i-1)+M/2)%label) hs(4) = var_stuff(max_depth+1, h1(2*i-1+M/2)%index, h1(2*i-1+M/2)%label) call WX(' '//write_var(hs(1))//' = '// & write_const(factor_data(n-1)%factor(j))//'*'// & write_var(h1(2*(i-1))),iunit) call WX(' '//write_var(hs(2))//' = '// & write_const(factor_data(n-1)%factor(j))//'*'// & write_var(h1(2*i-1)),iunit) call WX(' '//write_var(hs(3))//' = '// & write_const(factor_data(n-1)%factor(j))//'*'// & write_var(h1(2*(i-1)+M/2)),iunit) call WX(' '//write_var(hs(4))//' = '// & write_const(factor_data(n-1)%factor(j))//'*'// & write_var(h1(2*i-1+M/2)),iunit) if(i <= M/8) then call WX(' a('//write_number(m0)//') = '// & write_var(hs(1)),iunit) call WX(' a('//write_number(m0+1)//') = '// & write_var(hs(2))//'-'//write_var(hs(1)),iunit) call WX(' a('//write_number(m0+2)//') = '// & write_const(rot_mat_data(n-1)%rot_mats(i)%cos)//'*'// & write_var(hs(2))//'-'// & write_const(rot_mat_data(n-1)%rot_mats(i)%sin)//'*'// & write_var(hs(4))//'-'//write_var(hs(1)),iunit) call WX(' a('//write_number(m0+3)//') = '// & write_var(hs(3)),iunit) call WX(' a('//write_number(m0+4)//') = '// & write_var(hs(4))//'-'//write_var(hs(3)),iunit) call WX(' a('//write_number(m0+5)//') = '// & write_const(rot_mat_data(n-1)%rot_mats(i)%sin)//'*'// & write_var(hs(2))//'+'// & write_const(rot_mat_data(n-1)%rot_mats(i)%cos)//'*'// & write_var(hs(4))//'-'//write_var(hs(3)),iunit) else call WX(' a('//write_number(m0)//') = '// & write_var(hs(1)),iunit) call WX(' a('//write_number(m0+1)//') = '// & write_var(hs(2))//'-'//write_var(hs(1)),iunit) call WX(' a('//write_number(m0+2)//') = '// & write_const(rot_mat_data(n-1)%rot_mats(M/4+1-i)%cos)//'*'// & write_var(hs(2))//'+'// & write_const(rot_mat_data(n-1)%rot_mats(M/4+1-i)%sin)//'*'// & write_var(hs(4))//'+'//write_var(hs(1)),iunit) call WX(' a('//write_number(m0+3)//') = '// & write_var(hs(3)),iunit) call WX(' a('//write_number(m0+4)//') = '// & write_var(hs(4))//'-'//write_var(hs(3)),iunit) call WX(' a('//write_number(m0+5)//') = '// & write_const(rot_mat_data(n-1)%rot_mats(M/4+1-i)%sin)//'*'// & write_var(hs(2))//'-'// & write_const(rot_mat_data(n-1)%rot_mats(M/4+1-i)%cos)//'*'// & write_var(hs(4))//'-'//write_var(hs(3)),iunit) end if m0 = m0+6 end do end if end if end subroutine conv_const !................................................................ !................................................................ ! Subroutine conv_eval_correct(h,n,iunit) ! Evaluates real convolution of order 2**n < 8 and immortalizes ! results in code !................................................................ subroutine conv_eval_correct(h,n,iunit) integer, intent(in) :: n ! Order of the problem = 2**n type(var_stuff), intent(inout) :: h(0:2**n-1) ! Array of variable mnemonics integer, intent(in) :: iunit ! Write code to this unit integer max_depth ! Depth of previous problem type(var_stuff) hs(4) ! Array of variable mnemonics select case(n) case(0) call WX(' h(0) = a(0)*'//write_input(h(0)),iunit) case(1) hs(1) = var_stuff(1,h(0)%index,h(0)%label) hs(2) = var_stuff(1,h(1)%index,h(1)%label) call WX(' '//write_var(hs(1))//' = a(0)*h(0)+a(1)*h(1)',iunit) call WX(' '//write_var(hs(2))//' = a(1)*h(0)+a(0)*h(1)',iunit) call WX(' h(0) = '//write_var(hs(1)),iunit) call WX(' h(1) = '//write_var(hs(2)),iunit) case(2) max_depth = maxval(h%depth) if(max_depth == 0) then hs(1) = var_stuff(max_depth+1,h(0)%index,h(0)%label) hs(2) = var_stuff(max_depth+1,h(1)%index,h(1)%label) hs(3) = var_stuff(max_depth+1,h(2)%index,h(2)%label) hs(4) = var_stuff(max_depth+1,h(3)%index,h(3)%label) call WX(' '//write_var(hs(1))//' = h(0)+h(2)',iunit) call WX(' '//write_var(hs(2))//' = h(1)+h(3)',iunit) call WX(' '//write_var(hs(3))//' = h(0)-h(2)',iunit) call WX(' '//write_var(hs(4))//' = h(1)-h(3)',iunit) max_depth = maxval(hs%depth) h(0) = var_stuff(max_depth+1,hs(1)%index,hs(1)%label) h(1) = var_stuff(max_depth+1,hs(2)%index,hs(2)%label) h(2) = var_stuff(max_depth+1,hs(3)%index,hs(3)%label) h(3) = var_stuff(max_depth+1,hs(4)%index,hs(4)%label) call WX(' '//write_var(h(0))//' = a(0)*'// & write_var(hs(1))//'+a(1)*'//write_var(hs(2)),iunit) call WX(' '//write_var(h(1))//' = a(1)*'// & write_var(hs(1))//'+a(0)*'//write_var(hs(2)),iunit) call WX(' '//write_var(h(2))//' = a(2)*'// & write_var(hs(3))//'-a(3)*'//write_var(hs(4)),iunit) call WX(' '//write_var(h(3))//' = a(3)*'// & write_var(hs(3))//'+a(2)*'//write_var(hs(4)),iunit) call WX(' h(0) = '//write_var(h(0))//'+'// & write_var(h(2)),iunit) call WX(' h(1) = '//write_var(h(1))//'+'// & write_var(h(3)),iunit) call WX(' h(2) = '//write_var(h(0))//'-'// & write_var(h(2)),iunit) call WX(' h(3) = '//write_var(h(1))//'-'// & write_var(h(3)),iunit) else hs(1) = var_stuff(max_depth+1,h(0)%index,h(0)%label) hs(2) = var_stuff(max_depth+1,h(1)%index,h(1)%label) hs(3) = var_stuff(max_depth+1,h(2)%index,h(2)%label) hs(4) = var_stuff(max_depth+1,h(3)%index,h(3)%label) call WX(' '//write_var(hs(1))//' = '// & write_var(h(0))//'+'//write_var(h(2)),iunit) call WX(' '//write_var(hs(2))//' = '// & write_var(h(1))//'+'//write_var(h(3)),iunit) call WX(' '//write_var(hs(3))//' = '// & write_var(h(0))//'-'//write_var(h(2)),iunit) call WX(' '//write_var(hs(4))//' = '// & write_var(h(1))//'-'//write_var(h(3)),iunit) max_depth = maxval(hs%depth) h(0) = var_stuff(max_depth+1,hs(1)%index,hs(1)%label) h(1) = var_stuff(max_depth+1,hs(2)%index,hs(2)%label) h(2) = var_stuff(max_depth+1,hs(3)%index,hs(3)%label) h(3) = var_stuff(max_depth+1,hs(4)%index,hs(4)%label) call WX(' '//write_var(h(0))//' = a(0)*'// & write_var(hs(1))//'+a(1)*'//write_var(hs(2)),iunit) call WX(' '//write_var(h(1))//' = a(1)*'// & write_var(hs(1))//'+a(0)*'//write_var(hs(2)),iunit) call WX(' '//write_var(h(2))//' = a(2)*'// & write_var(hs(3))//'-a(3)*'//write_var(hs(4)),iunit) call WX(' '//write_var(h(3))//' = a(3)*'// & write_var(hs(3))//'+a(2)*'//write_var(hs(4)),iunit) max_depth = maxval(h%depth) hs(1) = var_stuff(max_depth+1,h(0)%index,h(0)%label) hs(2) = var_stuff(max_depth+1,h(1)%index,h(1)%label) hs(3) = var_stuff(max_depth+1,h(2)%index,h(2)%label) hs(4) = var_stuff(max_depth+1,h(3)%index,h(3)%label) call WX(' '//write_var(hs(1))//' = '// & write_var(h(0))//'+'//write_var(h(2)),iunit) call WX(' '//write_var(hs(2))//' = '// & write_var(h(1))//'+'//write_var(h(3)),iunit) call WX(' '//write_var(hs(3))//' = '// & write_var(h(0))//'-'//write_var(h(2)),iunit) call WX(' '//write_var(hs(4))//' = '// & write_var(h(1))//'-'//write_var(h(3)),iunit) h = hs end if case default stop ' bigger n than 2 in conv_cons_correct' end select end subroutine conv_eval_correct !................................................................ !................................................................ ! recursive Subroutine conv_eval(h,n,correct_data, & ! scaled_odd_tail_data,iunit) ! evaluates real convolution of order 2**n and immortalizes ! results in code !................................................................ recursive subroutine conv_eval(h,n,correct_data, & scaled_odd_tail_data,iunit) integer, intent(in) :: n ! Order of the problem = 2**n type(correct_stuff) correct_data ! Indices for degenerate problems type(scaled_odd_tail_stuff) scaled_odd_tail_data(0:n-2) ! Indices for ! scaled form of odd tail integer, intent(in) :: iunit ! Write code to this unit type(var_stuff), intent(inout) :: h(0:2**n-1) ! Array of variable mnemonics type(var_stuff) h0(0:2**(n-1)-1) ! Array of variable mnemonics type(var_stuff) h1(0:2**(n-1)-1) ! Array of variable mnemonics integer i ! Current index integer M ! Blocksize integer m0 ! Current constant array index type(var_stuff) hs(4) ! Array of variable mnemonics type(var_stuff) ht(4) ! Array of variable mnemonics integer max_depth ! Depth of previous problem logical outer ! Marks highest recursion level if(n < 3) then call conv_eval_correct(h,n,iunit) else max_depth = maxval(h%depth) outer = max_depth == 0 if(outer) then M = 2**(n-1) do i = 0, M-1 h0(i) = var_stuff(max_depth+1, h(i)%index, h(i)%label) call WX(' '//write_var(h0(i))//' = h('// & write_number(i)//')+h('//write_number(i+M)//')',iunit) end do call conv_eval(h0,n-1,correct_data, & scaled_odd_tail_data,iunit) do i = 0, M-1 h1(i) = var_stuff(max_depth+1, h(i+M)%index, h(i+M)%label) call WX(' '//write_var(h1(i))//' = h('// & write_number(i)//')-h('//write_number(i+M)//')',iunit) end do else M = 2**(n-1) do i = 0, M-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 call conv_eval(h0,n-1,correct_data, & scaled_odd_tail_data,iunit) do i = 0, M-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 end if if(n == 3) then max_depth = maxval(h1%depth) hs(1) = var_stuff(max_depth+1, h1(0)%index, h1(0)%label) hs(2) = var_stuff(max_depth+1, h1(1)%index, h1(1)%label) hs(3) = var_stuff(max_depth+2, h1(2)%index, h1(2)%label) hs(4) = var_stuff(max_depth+2, h1(3)%index, h1(3)%label) call WX(' '//write_var(hs(1))//' = '// & write_var(h1(0))//'+'//write_var(h1(1)),iunit) call WX(' '//write_var(hs(2))//' = '// & write_var(h1(2))//'+'//write_var(h1(3)),iunit) call WX(' '//write_var(hs(3))//' = a(4)*'// & write_var(hs(1))//'-a(7)*'//write_var(hs(2)),iunit) call WX(' '//write_var(hs(4))//' = a(7)*'// & write_var(hs(1))//'+a(4)*'//write_var(hs(2)),iunit) max_depth = maxval(hs%depth) ht(1) = var_stuff(max_depth+1, hs(1)%index, hs(1)%label) ht(2) = var_stuff(max_depth+1, hs(2)%index, hs(2)%label) ht(3) = var_stuff(max_depth+1, hs(3)%index, hs(3)%label) ht(4) = var_stuff(max_depth+1, hs(4)%index, hs(4)%label) call WX(' '//write_var(ht(1))//' = '//write_var(hs(3))// & '-a(6)*'//write_var(h1(1))//'-a(9)*'//write_var(h1(3)),iunit) call WX(' '//write_var(ht(2))//' = '//write_var(hs(3))// & '+a(5)*'//write_var(h1(0))//'-a(8)*'//write_var(h1(2)),iunit) call WX(' '//write_var(ht(3))//' = '//write_var(hs(4))// & '+a(9)*'//write_var(h1(1))//'-a(6)*'//write_var(h1(3)),iunit) call WX(' '//write_var(ht(4))//' = '//write_var(hs(4))// & '+a(8)*'//write_var(h1(0))//'+a(5)*'//write_var(h1(2)),iunit) h1 = ht else call scaled_odd_tail(h1(0::2),n-2,correct_data,scaled_odd_tail_data,iunit) call scaled_odd_tail(h1(1::2),n-2,correct_data,scaled_odd_tail_data,iunit) max_depth = maxval(h1%depth) m0 = 3*M/2-2 do i = 1, M/4 hs(1) = var_stuff(max_depth+1, h1(2*(i-1))%index, h1(2*(i-1))%label) hs(2) = var_stuff(max_depth+1, h1(2*i-1)%index, h1(2*i-1)%label) hs(3) = var_stuff(max_depth+2, h1(2*(i-1)+M/2)%index, h1(2*(i-1)+M/2)%label) hs(4) = var_stuff(max_depth+2, h1(2*i-1+M/2)%index, h1(2*i-1+M/2)%label) call WX(' '//write_var(hs(1))//' = '// & write_var(h1(2*(i-1)))//'+'//write_var(h1(2*i-1)),iunit) call WX(' '//write_var(hs(2))//' = '// & write_var(h1(2*(i-1)+M/2))//'+'//write_var(h1(2*i-1+M/2)),iunit) call WX(' '//write_var(hs(3))//' = a('//write_number(m0)// & ')*'//write_var(hs(1))//'-a('//write_number(m0+3)//')*'// & write_var(hs(2)),iunit) call WX(' '//write_var(hs(4))//' = a('//write_number(m0+3)// & ')*'//write_var(hs(1))//'+a('//write_number(m0)//')*'// & write_var(hs(2)),iunit) ht(1) = var_stuff(max_depth+3, hs(1)%index, hs(1)%label) ht(2) = var_stuff(max_depth+3, hs(2)%index, hs(2)%label) ht(3) = var_stuff(max_depth+3, hs(3)%index, hs(3)%label) ht(4) = var_stuff(max_depth+3, hs(4)%index, hs(4)%label) if(i <= M/8) then call WX(' '//write_var(ht(1))//' = '//write_var(hs(3))// & '+a('//write_number(m0+2)//')*'//write_var(h1(2*i-1))// & '-a('//write_number(m0+5)//')*'//write_var(h1(2*i-1+M/2)),iunit) call WX(' '//write_var(ht(2))//' = '//write_var(hs(3))// & '+a('//write_number(m0+1)//')*'//write_var(h1(2*(i-1)))// & '-a('//write_number(m0+4)//')*'//write_var(h1(2*(i-1)+M/2)),iunit) call WX(' '//write_var(ht(3))//' = '//write_var(hs(4))// & '+a('//write_number(m0+5)//')*'//write_var(h1(2*i-1))// & '+a('//write_number(m0+2)//')*'//write_var(h1(2*i-1+M/2)),iunit) call WX(' '//write_var(ht(4))//' = '//write_var(hs(4))// & '+a('//write_number(m0+4)//')*'//write_var(h1(2*(i-1)))// & '+a('//write_number(m0+1)//')*'//write_var(h1(2*(i-1)+M/2)),iunit) else call WX(' '//write_var(ht(1))//' = '//write_var(hs(3))// & '-a('//write_number(m0+2)//')*'//write_var(h1(2*i-1))// & '-a('//write_number(m0+5)//')*'//write_var(h1(2*i-1+M/2)),iunit) call WX(' '//write_var(ht(2))//' = '//write_var(hs(3))// & '+a('//write_number(m0+1)//')*'//write_var(h1(2*(i-1)))// & '-a('//write_number(m0+4)//')*'//write_var(h1(2*(i-1)+M/2)),iunit) call WX(' '//write_var(ht(3))//' = '//write_var(hs(4))// & '+a('//write_number(m0+5)//')*'//write_var(h1(2*i-1))// & '-a('//write_number(m0+2)//')*'//write_var(h1(2*i-1+M/2)),iunit) call WX(' '//write_var(ht(4))//' = '//write_var(hs(4))// & '+a('//write_number(m0+4)//')*'//write_var(h1(2*(i-1)))// & '+a('//write_number(m0+1)//')*'//write_var(h1(2*(i-1)+M/2)),iunit) end if h1(2*(i-1)) = ht(1) h1(2*i-1) = ht(2) h1(2*(i-1)+M/2) = ht(3) h1(2*i-1+M/2) = ht(4) m0 = m0+6 end do call scaled_odd_tail_transpose(h1(0::2),n-2,correct_data, & scaled_odd_tail_data,iunit) call scaled_odd_tail_transpose(h1(1::2),n-2,correct_data, & scaled_odd_tail_data,iunit) end if if(outer) then M = 2**(n-1) do i = 0, M-1 call WX(' h('//write_number(i)//') = '// & write_var(h0(i))//'+'//write_var(h1(i)),iunit) end do do i = 0, M-1 call WX(' h('//write_number(i+M)//') = '// & write_var(h0(i))//'-'//write_var(h1(i)),iunit) end do else M = 2**(n-1) max_depth = max(maxval(h0%depth),maxval(h1%depth)) do i = 0, M-1 h(i) = var_stuff(max_depth+1,h0(i)%index,h0(i)%label) call WX(' '//write_var(h(i))//' = '// & write_var(h0(i))//'+'//write_var(h1(i)),iunit) end do do i = 0, M-1 h(i+M) = var_stuff(max_depth+1,h1(i)%index,h1(i)%label) call WX(' '//write_var(h(i+M))//' = '// & write_var(h0(i))//'-'//write_var(h1(i)),iunit) end do end if end if end subroutine conv_eval !................................................................ !................................................................ ! 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 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 !................................................................ !................................................................ ! function write_number(var) ! Write number !................................................................ function write_number(var) integer, intent(in) :: var ! Number ! The following declaration triggers a compiler bug on CVF 6.6C. ! The workaround was found by James Giles 8/20/04 and is ! commented out below. character(write_number_len(var)) write_number ! Corresponding string ! character(write_number_len((var))) write_number ! Corresponding string write(write_number,'(i0)') var end function write_number !................................................................ !................................................................ ! function write_number_len(var) ! Specification function to compute len(write_number(var)) !................................................................ pure function write_number_len(var) integer, intent(in) :: var ! number integer write_number_len ! Len(write_number(var)) character(80) test_string ! Holds string temporarily write(test_string,'(i0)') var write_number_len = len_trim(test_string) end function write_number_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 eval*.f90 file character(80) filename2 ! Name of output const*.i90 file character(80) filename3 ! Name of output conv*.f90 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(factor_stuff), allocatable :: factor_data(:) ! Indices for block factors type(rot_mat_stuff), allocatable :: rot_mat_data(:) ! Indices for rotation matrices type(scaled_odd_tail_stuff), allocatable :: scaled_odd_tail_data(:) ! Indices ! for scaled odd tail problem type(avl_bag) body ! Holds actual numeric data integer iunit ! write code to this unit type(var_stuff), allocatable :: h(:) ! Initial input array mnemonics integer Na ! Order of constants array 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(*,'(a)') ' Creating constants code:' call cpu_time(t0) ! Write constant subroutine file write(filename2,'(a,i0,a)') 'const', 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 conv2c.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 const', 2**n, '(h,a)' if(n < 3) then Na = 2*n else Na = 3*2**(n-1)-2 end if write(iunit,'(a)') ' implicit real(wp) (x,y)' write(iunit,'(a,i0,a)') ' real(wp), intent(in) :: h(0:', 2**n-1, ')' write(iunit,'(a,i0,a)') ' real(wp), intent(out) :: a(0:', Na-1, ')' allocate(get_factor_data(0:n-2)) do i = 0, n-2 allocate(get_factor_data(i)%factor(2**(i-3))) if(i >= 3) then call get_factor(get_factor_data(i)%factor,i,get_factor_data(i-2)%factor) end if end do allocate(factor_data(n-1),stat=i) call prepare_factors(n, body, factor_data, get_factor_data) allocate(rot_mat_data(n-1)) call prepare_rot_mat(n, body, rot_mat_data) call prepare_correct(n, body, correct_data) allocate(scaled_odd_tail_data(0:n-2)) do i = 0, 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) allocate(h(2**n)) do i = 1, size(h) h(i) = var_stuff(0,i-1,'xr') end do call conv_const(h,n,correct_data,factor_data,rot_mat_data,scaled_odd_tail_data,iunit) call cpu_time(t1) write(*,'(a,f0.3)') ' Execution time = ', t1-t0 write(iunit,'(a,i0)') 'end subroutine const', 2**n close(iunit) write(*,'(a)') ' Creating evaluation code:' call cpu_time(t0) ! Write evaluation subroutine file write(filename1,'(a,i0,a)') 'eval', 2**n, '.i90' iunit = 11 open(iunit,file=filename1,status='replace') write(iunit,'(a,a)') '! File: ', trim(filename1) call date_and_time(date, time) write(iunit,'(12a)') '! Generated by conv2c.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 eval', 2**n, '(h,a)' An = (8*3*n*2**n-31*2**n+4*(-1)**n)/9+4 if(n == 1) then Mn = 4 else Mn = (3*10*n*2**n-43*2**n+4*3*n*(-1)**n-38*(-1)**n)/27+6 end if write(iunit,'(a,i0,a,i0,a)') '! ',An,' adds, ',Mn,' muls' write(iunit,'(a)') ' implicit real(wp) (x,y)' write(iunit,'(a,i0,a)') ' real(wp), intent(inout) :: h(0:', 2**n-1, ')' write(iunit,'(a,i0,a)') ' real(wp), intent(in) :: a(0:', Na-1, ')' call prepare_correct(n, body, correct_data) do i = 0, 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) do i = 1, size(h) h(i) = var_stuff(0,i-1,'xr') end do call conv_eval(h,n,correct_data,scaled_odd_tail_data,iunit) call cpu_time(t1) write(*,'(a,f0.3)') ' Execution time = ', t1-t0 write(iunit,'(a,i0)') 'end subroutine eval', 2**n close(iunit) ! Write test program write(filename3,'(a,i0,a)') 'conv', 2**n, '.f90' iunit = 12 open(iunit,file=filename3,status='replace') write(iunit,'(a,a)') '! File: ', trim(filename1) call date_and_time(date, time) write(iunit,'(12a)') '! Generated by conv2c.f90 ', date(5:6), & '/', date(7:8), '/', date(1:4), ' ', time(1:2), & ':', time(3:4), ':', time(5:10) write(iunit,'(a)') 'module mykinds' write(iunit,'(a)') ' implicit none' write(iunit,'(a)') ' integer, parameter :: dp = selected_real_kind(15,300)' write(iunit,'(a)') 'end module mykinds' write(iunit,'(a)') '' write(iunit,'(a)') 'module test' write(iunit,'(a)') ' use mykinds, only: wp => dp' write(iunit,'(a)') ' implicit none' write(iunit,'(a)') ' private wp' write(iunit,'(a)') ' contains' write(iunit,'(a,i0,a)') "include 'const",2**n,".i90'" write(iunit,'(a,i0,a)') "include 'eval",2**n,".i90'" write(iunit,'(a)') 'end module test' write(iunit,'(a)') '' write(iunit,'(a)') 'program conv_test' write(iunit,'(a)') ' use mykinds' write(iunit,'(a)') ' use test' write(iunit,'(a)') ' implicit none' write(iunit,'(a,i0)') ' integer, parameter :: k = ',n write(iunit,'(a)') ' integer, parameter :: n = 2**k' write(iunit,'(a)') ' integer, parameter :: nta = n+(1+sign(1,n-4))/2*(n/2-2)' write(iunit,'(a)') ' real(dp) ta(0:nta-1)' write(iunit,'(a)') ' integer i' write(iunit,'(a)') ' integer j' write(iunit,'(a)') ' real(dp) h0(0:n-1)' write(iunit,'(a)') ' real(dp) h1(0:n-1)' write(iunit,'(a)') ' real(dp) error' write(iunit,'(a)') ' real(dp) max_error' write(iunit,'(a)') '' write(iunit,'(a)') ' max_error = 0' write(iunit,'(a)') ' do i = 0, n-1' write(iunit,'(a)') ' h0 = 0' write(iunit,'(a)') ' h0(i) = 1' write(iunit,'(a,i0,a)') ' call const',2**n,'(h0,ta)' write(iunit,'(a)') ' error = 0' write(iunit,'(a)') ' do j = 0, n-1' write(iunit,'(a)') ' h0 = 0' write(iunit,'(a)') ' h0(j) = 1' write(iunit,'(a,i0,a)') ' call eval',2**n,'(h0,ta)' write(iunit,'(a)') ' h1 = 0' write(iunit,'(a)') ' h1(modulo(i+j,n)) = 1' write(iunit,'(a)') ' error = max(error,maxval(abs(h0-h1)))' write(iunit,'(a)') ' end do' write(iunit,'(a)') ' write(*,*) i,error' write(iunit,'(a)') ' max_error = max(error,max_error)' write(iunit,'(a)') ' end do' write(iunit,'(a)') ' write(*,*)' write(iunit,'(a)') " write(*,*) 'max_error = ',max_error" write(iunit,'(a)') 'end program conv_test' close(iunit) write(*,'(a)') ' Files created: '//trim(filename1)// & ', '//trim(filename2)//', '//trim(filename3) end program test ! ***************************************************************