gcc/testsuite/gfortran.dg/graphite/pr40982.f90

   1 ! { dg-options "-O3 -fgraphite-identity -floop-interchange " }
   2
   3 module mqc_m
   4
   5
   6 implicit none
   7
   8 private
   9 public :: mutual_ind_quad_cir_coil
  10
  11 integer, parameter, private :: longreal = selected_real_kind(15,90)
  12 real (kind = longreal), parameter, private :: pi = 3.141592653589793_longreal
  13 real (kind = longreal), parameter, private :: small = 1.0e-10_longreal
  14
  15 contains
  16
  17       subroutine mutual_ind_quad_cir_coil (r_coil, x_coil, y_coil, z_coil, h_coil, n_coil,  &
  18                                                       rotate_coil, m, mu, l12)
  19       real (kind = longreal), intent(in) :: r_coil, x_coil, y_coil, z_coil, h_coil, n_coil, &
  20                                             mu
  21       real (kind = longreal), dimension(:,:), intent(in) :: rotate_coil
  22       integer, intent(in) :: m
  23       real (kind = longreal), intent(out) :: l12
  24       real (kind = longreal), dimension(3,3) :: rotate_quad
  25       real (kind = longreal), dimension(9), save :: x2gauss, y2gauss, w2gauss, z1gauss,     &
  26                                                     w1gauss
  27       real (kind = longreal) :: xxvec, xyvec, xzvec, yxvec, yyvec, yzvec, zxvec, zyvec,     &
  28                                 zzvec, magnitude, l12_lower, l12_upper, dx, dy, dz, theta,  &
  29                                 a, b1, b2, numerator, denominator, coefficient, angle
  30       real (kind = longreal), dimension(3) :: c_vector, q_vector, rot_c_vector,             &
  31                                               rot_q_vector, current_vector,                 &
  32                                               coil_current_vec, coil_tmp_vector
  33       integer :: i, j, k
  34       logical, save :: first = .true.
  35
  36       do i = 1, 2*m
  37           theta = pi*real(i,longreal)/real(m,longreal)
  38           c_vector(1) = r_coil * cos(theta)
  39           c_vector(2) = r_coil * sin(theta)
  40           coil_tmp_vector(1) = -sin(theta)
  41           coil_tmp_vector(2) = cos(theta)
  42           coil_tmp_vector(3) = 0.0_longreal
  43           coil_current_vec(1) = dot_product(rotate_coil(1,:),coil_tmp_vector(:))
  44           coil_current_vec(2) = dot_product(rotate_coil(2,:),coil_tmp_vector(:))
  45           coil_current_vec(3) = dot_product(rotate_coil(3,:),coil_tmp_vector(:))
  46           do j = 1, 9
  47               c_vector(3) = 0.5 * h_coil * z1gauss(j)
  48               rot_c_vector(1) = dot_product(rotate_coil(1,:),c_vector(:)) + dx
  49               rot_c_vector(2) = dot_product(rotate_coil(2,:),c_vector(:)) + dy
  50               rot_c_vector(3) = dot_product(rotate_coil(3,:),c_vector(:)) + dz
  51               do k = 1, 9
  52                   q_vector(1) = 0.5_longreal * a * (x2gauss(k) + 1.0_longreal)
  53                   q_vector(2) = 0.5_longreal * b1 * (y2gauss(k) - 1.0_longreal)
  54                   q_vector(3) = 0.0_longreal
  55                   rot_q_vector(1) = dot_product(rotate_quad(1,:),q_vector(:))
  56                   rot_q_vector(2) = dot_product(rotate_quad(2,:),q_vector(:))
  57                   rot_q_vector(3) = dot_product(rotate_quad(3,:),q_vector(:))
  58                   numerator = w1gauss(j) * w2gauss(k) *                                     &
  59                                                  dot_product(coil_current_vec,current_vector)
  60                   denominator = sqrt(dot_product(rot_c_vector-rot_q_vector,                 &
  61                                                                   rot_c_vector-rot_q_vector))
  62                   l12_lower = l12_lower + numerator/denominator
  63               end do
  64           end do
  65       end do
  66       l12 = coefficient * (b1 * l12_lower + b2 * l12_upper)
  67       end subroutine mutual_ind_quad_cir_coil
  68
  69 end module mqc_m