9 #include "sll_memory.h" 
   10 #include "sll_working_precision.h" 
   28       sll_int32, 
dimension(:, :, :), 
pointer  :: index 
 
   50       sll_real64, 
dimension(:), 
intent(in) :: 
data  
   51       sll_real64, 
dimension(:), 
intent(in) :: eta 
 
   53       if (interpolator%boundary == 0) 
then 
   67       sll_real64, 
dimension(:), 
intent(inout) :: data_in 
 
   68       sll_real64, 
dimension(:), 
intent(out) :: data_out 
 
   69       sll_int32, 
dimension(:), 
intent(in) :: dorder 
 
   70       sll_real64, 
intent(in) ::displacement 
 
   71       logical, 
intent(in) :: hiera
 
   73       sll_int32 :: i1, i2, i3, k2, k3, counter, j
 
   74       sll_int32, 
dimension(3) :: l, no, ind_order
 
   75       sll_int32, 
dimension(:, :), 
allocatable :: ind
 
   77       sll_allocate(ind(interpolator%max_level + 1, 3), i1); 
 
   78       ind_order(dorder(1)) = 0
 
   80       do i3 = 0, interpolator%levels(dorder(3))
 
   81          ind_order(dorder(3)) = i3
 
   83          no(dorder(3)) = max(2**(i3 - 1), 1); 
 
   84          do i2 = 0, min(interpolator%max_level - i3, interpolator%levels(dorder(2)))
 
   85             ind_order(dorder(2)) = i2
 
   86             no(dorder(2)) = max(2**(i2 - 1), 1); 
 
   87             ind(1, dorder(1)) = 0; 
 
   88             do k2 = 0, no(dorder(2)) - 1
 
   89                ind(ind_order(dorder(2)) + 1, dorder(2)) = k2; 
 
   90                do k3 = 0, no(dorder(3)) - 1
 
   91                   ind(ind_order(dorder(3)) + 1, dorder(3)) = k3; 
 
   92                   counter = interpolator%index( &
 
   93                             ind_order(1), ind_order(2), ind_order(3)) + &
 
   94                             ind(ind_order(1) + 1, 1)*no(2)*no(3) &
 
   95                             + ind(ind_order(2) + 1, 2)*no(3) + &
 
   96                             ind(ind_order(3) + 1, 3)
 
   98                   call interpolator%interpolate_disp_1d_periodic &
 
   99                      (displacement, dorder(1), &
 
  100                       min(interpolator%levels(dorder(1)), &
 
  101                           interpolator%max_level - ind_order(dorder(2)) - &
 
  102                           ind_order(dorder(3))), counter, data_in, data_out, hiera)
 
  108       if (hiera .EQV. .false.) 
then 
  110          do j = interpolator%order, 2, -1
 
  111             call interpolator%dehierarchical_part_order &
 
  113                 interpolator%dim, 2, dorder, j)
 
  116          call interpolator%dehierarchical_part(data_out, &
 
  117                                                interpolator%dim, 2, dorder)
 
  126       sll_int32 :: j, l1, l2, l3, level
 
  128       sll_real64, 
dimension(:), 
intent(in) :: 
data, eta
 
  129       sll_real64, 
dimension(3) :: eta_norm
 
  130       sll_real64, 
dimension(3) :: phi
 
  131       sll_int32, 
dimension(3) :: no, l
 
  132       sll_int32, 
dimension(:, :), 
allocatable :: ind
 
  136       sll_allocate(ind(0:interpolator%max_level, 1:3), j)
 
  141       do j = 1, interpolator%dim
 
  142          eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
 
  143          eta_norm(j) = modulo(eta_norm(j), 1.0_f64)
 
  146          do level = 2, interpolator%max_level
 
  147             ind(level, j) = ind(level - 1, j)*2
 
  148             if (eta_norm(j) > scale*real(ind(level, j) + 1, f64)) 
then 
  149                ind(level, j) = ind(level, j) + 1
 
  151             scale = scale*0.5_f64
 
  155       do level = 0, interpolator%max_level
 
  156          do l1 = 0, min(level, interpolator%levels(1))
 
  158             no(1) = max(2**(l1 - 1), 1)
 
  159             do l2 = max(0, level - l1 - interpolator%levels(3)), min(level - l1, interpolator%levels(2))
 
  161                no(2) = max(2**(l2 - 1), 1)
 
  162                l(3) = level - l1 - l2
 
  164                no(3) = max(2**(l(3) - 1), 1)
 
  166                index = interpolator%index(l1, l2, l3) + ind(l1, 1)*no(2)*no(3) &
 
  167                        + ind(l2, 2)*no(3) + ind(l3, 3)
 
  169                   call interpolator%basis_function(real(2**(max(l(j), 1)), f64)*eta_norm(j) &
 
  170                                                    - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
 
  171                                                    interpolator%hierarchy(index)%function_type(j))
 
  173                val = val + 
data(index)*phi(1)*phi(2)*phi(3)
 
  183       sll_int32 :: j, l1, l2, l3, level
 
  185       sll_real64, 
dimension(:), 
intent(in) :: 
data 
  186       sll_real64, 
dimension(:), 
intent(in) :: eta
 
  187       sll_real64, 
dimension(3) :: eta_norm
 
  188       sll_real64, 
dimension(3) :: phi
 
  189       sll_int32, 
dimension(3) :: no, l
 
  190       sll_int32, 
dimension(:, :), 
allocatable :: ind
 
  191       sll_real64 :: scale, phi1a, phi2a, phi3a
 
  194       sll_allocate(ind(0:interpolator%max_level, 1:interpolator%dim), j)
 
  197       ind(0:1, 1:interpolator%dim) = 0
 
  199       do j = 1, interpolator%dim
 
  200          eta_norm(j) = (eta(j) - interpolator%eta_min(j))/interpolator%length(j)
 
  203          do level = 2, interpolator%levels(j)
 
  204             ind(level, j) = ind(level - 1, j)*2
 
  205             if (eta_norm(j) > scale*real(ind(level, j) + 1, f64)) 
then 
  206                ind(level, j) = ind(level, j) + 1
 
  208             scale = scale*0.5_f64
 
  212       do level = 0, interpolator%max_level
 
  213          do l1 = 0, min(level, interpolator%levels(1))
 
  220             do l2 = max(0, level - l1 - interpolator%levels(3)), min(level - l1, interpolator%levels(2))
 
  225                   no(2) = max(2**(l2 - 1), 1); 
 
  227                l3 = level - l1 - l2; 
 
  231                   no(3) = max(2**(l3 - 1), 1); 
 
  236                         index = interpolator%index(l1, l2, l3) + &
 
  237                                 ind(l1, 1)*no(2)*no(3) &
 
  238                                 + ind(l2, 2)*no(3) + ind(l3, 3)
 
  240                         do j = 1, interpolator%dim
 
  241                            call interpolator%basis_function(real(2**l(j), f64)*eta_norm(j) &
 
  242                                                             - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
 
  243                                                             interpolator%hierarchy(index)%function_type(j))
 
  245                         val = val + 
data(index) &
 
  246                               *phi(1)*phi(2)*phi(3)
 
  248                         index = interpolator%index(l1, l2, l3) + &
 
  249                                 ind(l1, 1)*no(2)*no(3) &
 
  250                                 + ind(l2, 2)*no(3) + ind(l3, 3)
 
  251                         call interpolator%basis_function(real(2**l(1), f64)*eta_norm(1) &
 
  252                                                          - real(2*ind(l(1), 1), f64) - 1.0_f64, phi(1), &
 
  253                                                          interpolator%hierarchy(index)%function_type(1))
 
  254                         call interpolator%basis_function(real(2**l(2), f64)*eta_norm(2) &
 
  255                                                          - real(2*ind(l(2), 2), f64) - 1.0_f64, phi(2), &
 
  256                                                          interpolator%hierarchy(index)%function_type(2))
 
  257                         call interpolator%basis_function(eta_norm(3), phi(3), -1)
 
  258                         val = val + 
data(index)*phi(1)*phi(2)*phi(3); 
 
  259                         call interpolator%basis_function(eta_norm(3) - 1.0_f64, phi(3), -1)
 
  260                         val = val + 
data(index + 1)*phi(1)*phi(2)*phi(3); 
 
  264                         index = interpolator%index(l1, l2, l3) + &
 
  265                                 ind(l1, 1)*no(2)*no(3) &
 
  266                                 + ind(l2, 2)*no(3) + ind(l3, 3)
 
  267                         call interpolator%basis_function(real(2**l(1), f64)*eta_norm(1) &
 
  268                                                          - real(2*ind(l(1), 1), f64) - 1.0_f64, phi(1), &
 
  269                                                          interpolator%hierarchy(index)%function_type(1))
 
  270                         call interpolator%basis_function(real(2**l(3), f64)*eta_norm(3) &
 
  271                                                          - real(2*ind(l(3), 3), f64) - 1.0_f64, phi(3), &
 
  272                                                          interpolator%hierarchy(index)%function_type(3))
 
  273                         call interpolator%basis_function(eta_norm(2), phi(2), -1)
 
  274                         val = val + 
data(index)*phi(1)*phi(2)*phi(3); 
 
  275                         call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi(2), -1)
 
  276                         val = val + 
data(index + no(3))*phi(1)*phi(2)*phi(3); 
 
  278                         index = interpolator%index(l1, l2, l3) + &
 
  279                                 ind(l1, 1)*no(2)*no(3) + ind(l2, 2)*no(3) + ind(l3, 3); 
 
  280                         call interpolator%basis_function(real(2**l(1), f64)*eta_norm(1) &
 
  281                                                          - real(2*ind(l(1), 1), f64) - 1.0_f64, phi(1), &
 
  282                                                          interpolator%hierarchy(index)%function_type(1))
 
  283                         call interpolator%basis_function(eta_norm(3), phi(3), -1)
 
  284                         call interpolator%basis_function(eta_norm(3) - 1.0_f64, phi3a, -1)
 
  285                         call interpolator%basis_function(eta_norm(2), phi(2), -1)
 
  286                         call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi2a, -1)
 
  287                         val = val + 
data(index)*phi(1)*phi(2)*phi(3) + &
 
  288                               data(index + 1)*phi(1)*phi(2)*phi3a + &
 
  289                               data(index + 2)*phi(1)*phi2a*phi(3) + &
 
  290                               data(index + 3)*phi(1)*phi2a*phi3a
 
  296                         index = interpolator%index(l1, l2, l3) + &
 
  297                                 ind(l1, 1)*no(2)*no(3) &
 
  298                                 + ind(l2, 2)*no(3) + ind(l3, 3)
 
  300                         do j = 2, interpolator%dim
 
  301                            call interpolator%basis_function(real(2**l(j), f64)*eta_norm(j) &
 
  302                                                             - real(2*ind(l(j), j), f64) - 1.0_f64, phi(j), &
 
  303                                                             interpolator%hierarchy(index)%function_type(j))
 
  305                         call interpolator%basis_function(eta_norm(1), phi(1), -1)
 
  306                         val = val + 
data(index)*phi(1)*phi(2)*phi(3); 
 
  307                         call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi(1), -1)
 
  308                         val = val + 
data(index + no(2)*no(3))*phi(1)*phi(2)*phi(3)
 
  310                         index = interpolator%index(l1, l2, l3) + &
 
  311                                 ind(l1, 1)*no(2)*no(3) &
 
  312                                 + ind(l2, 2)*no(3) + ind(l3, 3)
 
  313                         call interpolator%basis_function(eta_norm(1), phi(1), -1)
 
  314                         call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi1a, -1)
 
  315                         call interpolator%basis_function(real(2**l(2), f64)*eta_norm(2) &
 
  316                                                          - real(2*ind(l(2), 2), f64) - 1.0_f64, phi(2), &
 
  317                                                          interpolator%hierarchy(index)%function_type(2))
 
  318                         call interpolator%basis_function(eta_norm(3), phi(3), -1)
 
  319                         call interpolator%basis_function(eta_norm(3) - 1.0_f64, phi3a, -1)
 
  320                         val = val + 
data(index)*phi(1)*phi(2)*phi(3) + &
 
  321                               data(index + 1)*phi(1)*phi(2)*phi3a + &
 
  322                               data(index + no(2)*no(3))*phi1a*phi(2)*phi(3) + &
 
  323                               data(index + no(2)*no(3) + 1)*phi1a*phi(2)*phi3a
 
  327                         index = interpolator%index(l1, l2, l3) + &
 
  328                                 ind(l1, 1)*no(2)*no(3) &
 
  329                                 + ind(l2, 2)*no(3) + ind(l3, 3)
 
  330                         call interpolator%basis_function(eta_norm(1), phi(1), -1)
 
  331                         call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi1a, -1)
 
  332                         call interpolator%basis_function(real(2**l(3), f64)*eta_norm(3) &
 
  333                                                          - real(2*ind(l(3), 3), f64) - 1.0_f64, phi(3), &
 
  334                                                          interpolator%hierarchy(index)%function_type(3))
 
  335                         call interpolator%basis_function(eta_norm(2), phi(2), -1)
 
  336                         call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi2a, -1)
 
  337                         val = val + 
data(index)*phi(1)*phi(2)*phi(3) + &
 
  338                               data(index + no(3))*phi(1)*phi2a*phi(3) + &
 
  339                               data(index + no(2)*no(3))*phi1a*phi(2)*phi(3) + &
 
  340                               data(index + no(2)*no(3) + no(3))*phi1a*phi2a*phi(3); 
 
  342                         index = interpolator%index(l1, l2, l3) + &
 
  343                                 ind(l1, 1)*no(2)*no(3) + ind(l2, 2)*no(3) + ind(l3, 3); 
 
  344                         call interpolator%basis_function(eta_norm(1), phi(1), -1)
 
  345                         call interpolator%basis_function(eta_norm(1) - 1.0_f64, phi1a, -1)
 
  346                         call interpolator%basis_function(eta_norm(3), phi(3), -1)
 
  347                         call interpolator%basis_function(eta_norm(3) - 1.0_f64, phi3a, -1)
 
  348                         call interpolator%basis_function(eta_norm(2), phi(2), -1)
 
  349                         call interpolator%basis_function(eta_norm(2) - 1.0_f64, phi2a, -1)
 
  350                         val = val + 
data(index)*phi(1)*phi(2)*phi(3) + &
 
  351                               data(index + 1)*phi(1)*phi(2)*phi3a + &
 
  352                               data(index + 2)*phi(1)*phi2a*phi(3) + &
 
  353                               data(index + 3)*phi(1)*phi2a*phi3a + &
 
  354                               data(index + 4)*phi1a*phi(2)*phi(3) + &
 
  355                               data(index + 5)*phi1a*phi(2)*phi3a + &
 
  356                               data(index + 6)*phi1a*phi2a*phi(3) + &
 
  357                               data(index + 7)*phi1a*phi2a*phi3a
 
  375       sll_int32, 
intent(in) :: dim 
 
  376       sll_real64, 
intent(in) :: displacment_in 
 
  377       sll_comp64, 
dimension(:), 
intent(inout)   :: data_in 
 
  378       sll_real64, 
dimension(:), 
intent(out) :: data_out 
 
  379       sll_real64:: displacement
 
  381       displacement = displacment_in*2.0_f64*
sll_p_pi/interpolator%length(dim)
 
  382       call displace(interpolator, dim, displacement, data_in); 
 
  383       call ispfft(interpolator, data_in, data_out); 
 
  398       interpolation_type, &
 
  405       sll_real64, 
dimension(:), 
intent(in)              :: eta_min 
 
  406       sll_real64, 
dimension(:), 
intent(in)             :: eta_max 
 
  407       sll_int32, 
dimension(:), 
intent(in)                :: levels
 
  409       sll_int32, 
intent(in)                             :: order 
 
  410       sll_int32, 
intent(in)                             :: interpolation 
 
  411       sll_int32, 
intent(in)                             :: interpolation_type 
 
  412       sll_int32, 
intent(in)                             :: modified
 
  417       sll_int32, 
intent(in)                             :: boundary
 
  420       sll_int32                                     :: i, j, k1, k2, k3, l1, l2, l3, l, counter
 
  421       sll_int32                                     :: ierr, no1, no2, no3
 
  422       sll_int32, 
dimension(:), 
allocatable          :: novec, lvec, kvec
 
  424       interpolator%dim = 3; 
 
  425       interpolator%modified = modified; 
 
  426       interpolator%boundary = boundary; 
 
  427       sll_allocate(lvec(interpolator%dim), ierr); 
 
  428       sll_allocate(kvec(interpolator%dim), ierr); 
 
  429       sll_allocate(novec(interpolator%dim), ierr); 
 
  430       interpolator%max_level = levels(1); 
 
  431       do l = 2, interpolator%dim
 
  432          interpolator%max_level = max(levels(l), interpolator%max_level); 
 
  434       if (interpolator%modified == 1) 
then 
  435          interpolator%max_level = interpolator%max_level + 2; 
 
  438       interpolator%size_basis = 0; 
 
  439       if (interpolator%boundary == 0) 
then 
  440          do l = 0, interpolator%max_level
 
  441             do l1 = 0, min(l, levels(1))
 
  442                do l2 = max(0, l - l1 - levels(3)), min(l - l1, levels(2))
 
  444                   interpolator%size_basis = &
 
  445                      interpolator%size_basis + &
 
  446                      max(2**(l1 - 1), 1)*max(2**(l2 - 1), 1)*max(2**(l3 - 1), 1)
 
  451          do l = 0, interpolator%max_level
 
  452             do l1 = 0, min(l, levels(1))
 
  453                do l2 = max(0, l - l1 - levels(3)), min(l - l1, levels(2))
 
  470                   interpolator%size_basis = &
 
  471                      interpolator%size_basis + no1*no2*no3; 
 
  477       call interpolator%initialize_sg(levels, order, interpolation, &
 
  478                                       interpolation_type, eta_min, eta_max); 
 
  479       sll_allocate(interpolator%index(0:interpolator%levels(1), 0:interpolator%levels(2), 0:interpolator%levels(3)), ierr)
 
  483       do l = 0, interpolator%max_level
 
  484          interpolator%level_mapping(l) = counter; 
 
  485          do l1 = 0, min(l, interpolator%levels(1))
 
  486             novec(1) = max(2**(l1 - 1), 1)
 
  487             if (interpolator%boundary == 1 .AND. l1 == 0) 
then 
  491             do l2 = max(0, l - l1 - interpolator%levels(3)), min(l - l1, interpolator%levels(2))
 
  492                novec(2) = max(2**(l2 - 1), 1)
 
  493                if (interpolator%boundary == 1 .AND. l2 == 0) 
then 
  497                lvec(3) = l - l1 - l2; 
 
  499                novec(3) = max(2**(l3 - 1), 1)
 
  500                if (interpolator%boundary == 1 .AND. l3 == 0) 
then 
  504                interpolator%index(l1, l2, l3) = counter
 
  505                do k1 = 0, novec(1) - 1
 
  507                   do k2 = 0, novec(2) - 1
 
  509                      do k3 = 0, novec(3) - 1
 
  511                         do j = 1, interpolator%dim
 
  512                            if (interpolator%boundary == 0) 
then 
  517                                  (interpolator, counter, j, lvec, kvec, novec); 
 
  520                         counter = counter + 1; 
 
  527       interpolator%level_mapping(interpolator%max_level + 1) = counter; 
 
  529       do i = 1, interpolator%size_basis
 
  530          do j = 1, interpolator%dim
 
  531             interpolator%hierarchy(i)%coordinate(j) = &
 
  532                interpolator%hierarchy(i)%coordinate(j)* &
 
  533                interpolator%length(j) + &
 
  534                interpolator%eta_min(j)
 
  546       sll_int32, 
intent(in) :: cdim 
 
  547       sll_int32, 
intent(in) :: counter 
 
  548       sll_int32, 
dimension(:), 
intent(in) :: lvecin 
 
  549       sll_int32, 
dimension(:), 
intent(in) :: kvecin 
 
  550       sll_int32, 
dimension(:), 
intent(in) :: novecin 
 
  552       sll_int32, 
dimension(3) :: lvec, kvec, novec
 
  553       sll_int32 :: jj, stride
 
  555       do jj = 1, interpolator%dim
 
  556          lvec(jj) = lvecin(jj); 
 
  557          kvec(jj) = kvecin(jj); 
 
  558          novec(jj) = novecin(jj); 
 
  562       interpolator%hierarchy(counter)%level(cdim) = ld; 
 
  563       interpolator%hierarchy(counter)%index_on_level(cdim) = kd; 
 
  566          interpolator%hierarchy(counter)%coordinate(cdim) = 0.0_f64
 
  567          interpolator%hierarchy(counter)%parent(stride) = &
 
  569          interpolator%hierarchy(counter)%parent(stride + 1) = &
 
  571          interpolator%hierarchy(counter)%function_type(cdim) = 0
 
  573          interpolator%hierarchy(counter)%coordinate(cdim) = &
 
  574             1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
 
  576          lvec(cdim) = lvec(cdim) - 1; 
 
  577          novec(cdim) = max(novec(cdim)/2, 1); 
 
  579          kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1)); 
 
  580          interpolator%hierarchy(counter)%parent( &
 
  581             modulo(kd, 2) + stride) = &
 
  582             interpolator%hierarchy( &
 
  583             interpolator%index(lvec(1), lvec(2), lvec(3)) + &
 
  584             kvec(1)*novec(2)*novec(3) + &
 
  586             kvec(3))%parent(stride)
 
  589          interpolator%hierarchy(counter)%parent( &
 
  590             modulo(kd + 1, 2) + stride) = &
 
  591             interpolator%index(lvec(1), lvec(2), lvec(3)) + &
 
  592             kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + &
 
  595          interpolator%hierarchy( &
 
  596             interpolator%hierarchy(counter)%parent( &
 
  597             modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
 
  599             interpolator%hierarchy( &
 
  600                interpolator%hierarchy(counter)%parent( &
 
  601                modulo(kd + 1, 2) + stride))%children(modulo(kd + 1, 2) + stride) = counter
 
  603          if (interpolator%order == 1) 
then 
  604             interpolator%hierarchy(counter)% &
 
  605                function_type(cdim) = -1
 
  606          elseif (ld == 1 .OR. interpolator%order == 2) 
then 
  607             interpolator%hierarchy(counter)% &
 
  608                function_type(cdim) = 1 + modulo(kd, 2)
 
  609          elseif (ld == 2 .OR. interpolator%order == 3) 
then  
  610             interpolator%hierarchy(counter)% &
 
  611                function_type(cdim) = 3 + modulo(kd, 4)
 
  613             interpolator%hierarchy(counter)% &
 
  614                function_type(cdim) = 7 + modulo(kd, 8)
 
  625       sll_int32, 
intent(in) :: cdim 
 
  626       sll_int32, 
intent(in) :: counter 
 
  627       sll_int32, 
dimension(:), 
intent(in) :: lvecin 
 
  628       sll_int32, 
dimension(:), 
intent(in) :: kvecin 
 
  629       sll_int32, 
dimension(:), 
intent(in) :: novecin 
 
  631       sll_int32, 
dimension(:), 
allocatable :: lvec, kvec, novec
 
  632       sll_int32 :: jj, stride
 
  634       sll_allocate(lvec(interpolator%dim), jj); 
 
  635       sll_allocate(kvec(interpolator%dim), jj); 
 
  636       sll_allocate(novec(interpolator%dim), jj); 
 
  637       do jj = 1, interpolator%dim
 
  638          lvec(jj) = lvecin(jj); 
 
  639          kvec(jj) = kvecin(jj); 
 
  640          novec(jj) = novecin(jj); 
 
  644       interpolator%hierarchy(counter)%level(cdim) = ld; 
 
  647          interpolator%hierarchy(counter)%coordinate(cdim) = real(kd, f64); 
 
  648          interpolator%hierarchy(counter)%parent(stride) = &
 
  650          interpolator%hierarchy(counter)%parent(stride + 1) = &
 
  652          interpolator%hierarchy(counter)%function_type(cdim) = -1
 
  654          interpolator%hierarchy(counter)%coordinate(cdim) = &
 
  655             1.0_f64/(2.0_f64**ld) + real(kd, f64)*1.0_f64/(2.0_f64**(ld - 1))
 
  657          lvec(cdim) = lvec(cdim) - 1; 
 
  658          novec(cdim) = novec(cdim)/2; 
 
  662             interpolator%hierarchy(counter)%parent(stride) = &
 
  663                interpolator%index(lvec(1), lvec(2), lvec(3)) + &
 
  664                kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + kvec(3)
 
  666                interpolator%hierarchy(counter)%parent( &
 
  667                   stride + 1) = interpolator%hierarchy(counter)%parent( &
 
  668                              stride) + novec(2)*novec(3)
 
  669             elseif (cdim == 2) 
then 
  670                interpolator%hierarchy(counter)%parent( &
 
  671                   stride + 1) = interpolator%hierarchy(counter)%parent( &
 
  674                interpolator%hierarchy(counter)%parent( &
 
  675                   stride + 1) = interpolator%hierarchy(counter)%parent( &
 
  678             interpolator%hierarchy(interpolator%hierarchy(counter)%parent(stride))%children &
 
  684             kvec(cdim) = modulo((kd + 1)/2, max(2**(ld - 2), 1)); 
 
  685             kvec(cdim) = (kd + 1)/2; 
 
  686             if (kvec(cdim) < novec(cdim)) 
then 
  687                interpolator%hierarchy(counter)%parent( &
 
  688                   modulo(kd, 2) + stride) = &
 
  689                   interpolator%hierarchy( &
 
  690                   interpolator%index(lvec(1), lvec(2), lvec(3)) + &
 
  691                   kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + kvec(3))% &
 
  694                kvec(cdim) = kvec(cdim) - 1; 
 
  695                interpolator%hierarchy(counter)%parent( &
 
  696                   modulo(kd, 2) + stride) = &
 
  697                   interpolator%hierarchy( &
 
  698                   interpolator%index(lvec(1), lvec(2), lvec(3)) + &
 
  699                   kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + kvec(3))% &
 
  704             interpolator%hierarchy(counter)%parent( &
 
  705                modulo(kd + 1, 2) + stride) = &
 
  706                interpolator%index(lvec(1), lvec(2), lvec(3)) + &
 
  707                kvec(1)*novec(2)*novec(3) + kvec(2)*novec(3) + kvec(3)
 
  709             interpolator%hierarchy( &
 
  710                interpolator%hierarchy(counter)%parent( &
 
  711                modulo(kd + 1, 2) + stride))%children(modulo(kd, 2) + stride) = counter
 
  714          if (interpolator%order == 1) 
then 
  715             interpolator%hierarchy(counter)% &
 
  716                function_type(cdim) = -1
 
  717          elseif (ld == 1 .OR. interpolator%order == 2) 
then 
  718             interpolator%hierarchy(counter)% &
 
  719                function_type(cdim) = 1 + modulo(kd, 2)
 
  720          elseif (ld == 2 .OR. interpolator%order == 3) 
then  
  721             interpolator%hierarchy(counter)% &
 
  722                function_type(cdim) = 3 + modulo(kd, 4)
 
  724             interpolator%hierarchy(counter)% &
 
  725                function_type(cdim) = 7 + modulo(kd, 8)
 
  736    subroutine fg_to_sg(interpolator, fg_values, sg_values)
 
  737       sll_real64, 
dimension(:, :, :), 
intent(in) :: fg_values
 
  738       sll_real64, 
dimension(:), 
intent(out) :: sg_values
 
  741       sll_int32, 
dimension(3) :: fg_ind
 
  743       do j = 1, interpolator%size_basis
 
  745          sg_values(j) = fg_values(fg_ind(1), fg_ind(2), fg_ind(3)); 
 
  752       sll_int32, 
intent(in) :: sg_index
 
  757       do j = 1, interpolator%dim
 
  758          fg_index(j) = 2**(interpolator%levels(j) - &
 
  759                            interpolator%hierarchy(sg_index)%level(j))* &
 
  760                        (1 + 2*interpolator%hierarchy(sg_index)%index_on_level(j)) + 1; 
 
  773       sll_real64, 
dimension(:), 
intent(in) :: data_in
 
  774       sll_comp64, 
dimension(:), 
intent(out) :: data_out
 
  775       sll_int32 :: i1, i2, i3, j
 
  778       do i2 = 0, interpolator%levels(2)
 
  779          do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
 
  780             do j = interpolator%index(0, i2, i3), &
 
  781                interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  782                call interpolator%ToHierarchical1D(1, &
 
  783                                                   min(interpolator%levels(1), interpolator%max_level - i2 - i3), &
 
  784                                                   j, data_in, data_out)
 
  789       do i1 = 0, interpolator%levels(1)
 
  790          do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
 
  791             do j = interpolator%index(i1, 0, i3), &
 
  792                interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  793                call interpolator%ToHierarchical1D_comp &
 
  795                    min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
 
  801       do i1 = 0, interpolator%levels(1)
 
  802          do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
 
  803             do j = interpolator%index(i1, i2, 0), &
 
  804                interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
 
  805                call interpolator%ToHierarchical1D_comp &
 
  807                    min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
 
  815    subroutine todehi(interpolator, data_array)
 
  817       sll_comp64, 
dimension(:), 
intent(inout) :: data_array
 
  818       sll_int32 :: i1, i2, i3, j
 
  820       do i1 = 0, interpolator%levels(1)
 
  821          do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
 
  822             do j = interpolator%index(i1, i2, 0), &
 
  823                interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
 
  824                call interpolator%ToDehi1D &
 
  826                    min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
 
  832       do i1 = 0, interpolator%levels(1)
 
  833          do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
 
  834             do j = interpolator%index(i1, 0, i3), &
 
  835                interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  836                call interpolator%ToDehi1D &
 
  838                    min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
 
  844       do i2 = 0, interpolator%levels(2)
 
  845          do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
 
  846             do j = interpolator%index(0, i2, i3), &
 
  847                interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  848                call interpolator%ToDehi1D(1, &
 
  849                                           min(interpolator%levels(1), interpolator%max_level - i2 - i3), &
 
  857    subroutine tohira(interpolator, data_array)
 
  859       sll_comp64, 
dimension(:), 
intent(inout) :: data_array
 
  860       sll_int32 :: i1, i2, i3, j
 
  862       do i2 = 0, interpolator%levels(2)
 
  863          do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
 
  864             do j = interpolator%index(0, i2, i3), &
 
  865                interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  866                call interpolator%ToHira1D(1, &
 
  867                                           min(interpolator%levels(1), interpolator%max_level - i2 - i3), &
 
  873       do i1 = 0, interpolator%levels(1)
 
  874          do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
 
  875             do j = interpolator%index(i1, 0, i3), &
 
  876                interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  877                call interpolator%ToHira1D &
 
  879                    min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
 
  885       do i1 = 0, interpolator%levels(1)
 
  886          do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
 
  887             do j = interpolator%index(i1, i2, 0), &
 
  888                interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
 
  889                call interpolator%ToHira1D &
 
  891                    min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
 
  899    subroutine tonodal(interpolator, data_in, data_out)
 
  901       sll_comp64, 
dimension(:), 
intent(inout) :: data_in
 
  902       sll_real64, 
dimension(:), 
intent(out) :: data_out
 
  903       sll_int32 :: i1, i2, i3, j
 
  905       do i1 = 0, interpolator%levels(1)
 
  906          do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
 
  907             do j = interpolator%index(i1, i2, 0), &
 
  908                interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
 
  909                call interpolator%ToNodal1D_comp &
 
  911                    min(interpolator%levels(3), interpolator%max_level - i1 - i2), j, &
 
  917       do i1 = 0, interpolator%levels(1)
 
  918          do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
 
  919             do j = interpolator%index(i1, 0, i3), &
 
  920                interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  921                call interpolator%ToNodal1D_comp &
 
  923                    min(interpolator%levels(2), interpolator%max_level - i1 - i3), j, &
 
  929       do i2 = 0, interpolator%levels(2)
 
  930          do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
 
  931             do j = interpolator%index(0, i2, i3), &
 
  932                interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  933                call interpolator%ToNodal1D(1, &
 
  934                                            min(interpolator%levels(1), interpolator%max_level - i2 - i3), &
 
  935                                            j, data_in, data_out)
 
  942    subroutine displace(interpolator, dim, displacement, data)
 
  944       sll_comp64, 
dimension(:), 
intent(inout) :: 
data 
  945       sll_real64, 
intent(in) :: displacement
 
  946       sll_int32, 
intent(in) :: dim
 
  947       sll_int32 :: i1, i2, i3, j
 
  950          do i2 = 0, interpolator%levels(2)
 
  951             do i3 = 0, min(interpolator%max_level - i2, interpolator%levels(3))
 
  952                do j = interpolator%index(0, i2, i3), &
 
  953                   interpolator%index(0, i2, i3) + max(2**(i2 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  954                   call interpolator%Displace1D(1, interpolator%levels(1) - i2 - i3, &
 
  955                                                j, displacement, data)
 
  959       else if (dim == 2) 
then 
  960          do i1 = 0, interpolator%levels(1)
 
  961             do i3 = 0, min(interpolator%max_level - i1, interpolator%levels(3))
 
  962                do j = interpolator%index(i1, 0, i3), &
 
  963                   interpolator%index(i1, 0, i3) + max(2**(i1 - 1), 1)*max(2**(i3 - 1), 1) - 1
 
  964                   call interpolator%Displace1D(2, interpolator%levels(2) - i1 - i3, &
 
  965                                                j, displacement, data)
 
  969       else if (dim == 3) 
then 
  970          do i1 = 0, interpolator%levels(1)
 
  971             do i2 = 0, min(interpolator%max_level - i1, interpolator%levels(2))
 
  972                do j = interpolator%index(i1, i2, 0), &
 
  973                   interpolator%index(i1, i2, 0) + max(2**(i1 - 1), 1)*max(2**(i2 - 1), 1) - 1
 
  974                   call interpolator%Displace1D(3, interpolator%levels(3) - i1 - i2, &
 
  975                                                j, displacement, data)
 
  984    subroutine spfft(interpolator, data_in, data_out)
 
  986       sll_real64, 
dimension(:), 
intent(in) :: data_in
 
  987       sll_comp64, 
dimension(:), 
intent(out) :: data_out
 
  990       call todehi(interpolator, data_out)
 
  995    subroutine ispfft(interpolator, data_in, data_out)
 
  997       sll_comp64, 
dimension(:), 
intent(inout) :: data_in
 
  998       sll_real64, 
dimension(:), 
intent(out) :: data_out
 
 1000       call tohira(interpolator, data_in)
 
 1001       call tonodal(interpolator, data_in, data_out)
 
Fortran module where set some physical and mathematical constants.
real(kind=f64), parameter, public sll_p_pi
Implementation of a 3D sparse grid with interpolation routines.
subroutine set_hierarchy_info_boundary(interpolator, counter, cdim, lvecin, kvecin, novecin)
Helfer function for initialization. Setting all the information needed for node counter of the sparse...
real(kind=f64) function interpolate_from_hierarchical_surplus_boundary(interpolator, data, eta)
implements interpolation from hierarchical surplus (interpolate_from_interpolant_value) non-periodic
subroutine interpolate_array_disp_sgfft(interpolator, dim, displacment_in, data_in, data_out)
Compute value at displaced grid points using trigonometric interpolation (based on SG FFT)
real(kind=f64) function interpolate_from_hierarchical_surplus(interpolator, data, eta)
Implements interpolate_from_interpolant_value for periodic sparse grid.
subroutine ispfft(interpolator, data_in, data_out)
Sparse grid inverse FFT.
subroutine initialize_sg3d(interpolator, levels, order, interpolation, interpolation_type, eta_min, eta_max, boundary, modified)
Initialization function. Set up the hierarchy of the sparse grid.
real(kind=f64) function interpolate_from_interpolant_value(interpolator, data, eta)
Compute the value of the sparse grid interpolant at position eta (using standard sparse grid interpol...
subroutine interpolate_const_disp(interpolator, dorder, displacement, data_in, data_out, hiera)
Interpolation function for interpolation at (constantly) displaced grid points; displacement only in ...
subroutine tohierarchical(interpolator, data_in, data_out)
subroutine tohira(interpolator, data_array)
integer(kind=i32) function, dimension(3) fg_index(interpolator, sg_index)
Compute the index of a sparse grid node on level "level" with index "index_on_level" on full grid wit...
subroutine fg_to_sg(interpolator, fg_values, sg_values)
Functions to evaluate fg on sg and sg on fg.
subroutine spfft(interpolator, data_in, data_out)
Sparse grid FFT.
subroutine displace(interpolator, dim, displacement, data)
subroutine set_hierarchy_info(interpolator, counter, cdim, lvecin, kvecin, novecin)
Helfer function for initialization. Setting all the information needed for node counter of the sparse...
subroutine todehi(interpolator, data_array)
subroutine tonodal(interpolator, data_in, data_out)
Dimension-independent functions for sparse grid with polynomial basis functions.
Sparse grid object for 3d with interpolation routines.
Class defining the sparse grid data structure.