PROGRAM TEST_BUG
  integer :: i
  real :: d0(2), H(2,2), param(5), x(2), y(2), d(2), x0(2) = [3, 2]
  !
  print *, 'Test of xROTMG; y2 expected to be 0.0 for any inputs d1 and d2'
  !
  d = [1.0, 1.0]
  x = x0
  print *, 'd1 = ', d(1), ' d2 = ', d(2)
  call srotmg(d(1), d(2), x(1), x(2), param)
  call par2h(param, H)
  y = matmul(H, x0)
  print *, 'y2 = ', y(2)
  !
  d = [1e-8, 1e-8]
  x = x0
  print *, 'd1 = ', d(1), ' d2 = ', d(2)
  call srotmg(d(1), d(2), x(1), x(2), param)
  call par2h(param, H)
  y = matmul(H, x0)
  print *, 'y2 = ', y(2)
    
CONTAINS
  subroutine par2h(p, H)
    ! Construct H matrix according to documentation of xROTMG
    real :: p(5), H(2,2)
    select case(nint(p(1)))
    case (1)
      H(1,1) = p(2)
      H(1,2) = 1.0
      H(2,1) = -1.0
      H(2,2) = p(5)
    case(0)
      H(1,1) = 1.0
      H(1,2) = p(4)
      H(2,1) = p(3)
      H(2,2) = 1.0
    case(-1)
      H(1,1) = p(2)
      H(1,2) = p(4)
      H(2,1) = p(3)
      H(2,2) = p(5)
    case(-2)
      H(1,1) = 1.0
      H(1,2) = 0.0
      H(2,1) = 0.0
      H(2,2) = 1.0
    end select
  end subroutine par2h

END PROGRAM TEST_BUG