subroutine detq ( a, n, d, ifault )

!*****************************************************************************80
!
!! detq() computes the determinant of an orthogonal matrix.
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 August 2021
!
!  Author:
!
!    Original FORTRAN77 version by J C Gower.
!    FORTRAN90 version by John Burkardt
!
!  Reference:
!
!    J C Gower,
!    Algorithm AS 82:
!    The determinant of an orthogonal matrix,
!    Applied Statistics,
!    Volume 24, Number 1, 1975, page 150-153.
!
!  Input:
!
!    real ( kind = rk ) A(N,N), the orthogonal matrix stored by rows or columns.
!
!    integer N, the order of the matrix.
!
!  Output:
!
!    real ( kind = rk ) D, the determinant of A.
!
!    integer IFAULT, 
!    0, no error occurred.
!    1, an error was detected.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  real ( kind = rk ) a(n,n)
  real ( kind = rk ) a2(n*n)
  real ( kind = rk ) d
  integer i
  integer ifault
  integer j
  integer k
  integer p
  integer q
  integer r
  integer s
  real ( kind = rk ) tol
  real ( kind = rk ) x
  real ( kind = rk ) y

  ifault = 0
  tol = 0.0001D+00

  if ( n <= 0 ) then
    write ( *, '(a)' ) ''
    write ( *, '(a)' ) 'detq - Fatal error!'
    write ( *, '(a)' ) '  n <= 0'
    ifault = 1
    return
  end if

  k = 0
  do i = 1, n
    do j = 1, n
      k = k + 1
      a2(k) = a(i,j)
    end do
  end do

  d = 1.0D+00
  r = 1

  do k = 2, n + 1

    q = r
    x = a2(r)
    y = sign ( 1.0D+00, x )
    d = d * y
    y = - 1.0D+00 / ( x + y )
    x = abs ( x ) - 1.0D+00

    if ( tol < abs ( x ) ) then

      if ( 0.0D+00 < x ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) 'detq - Fatal error!'
        write ( *, '(a)' ) '  x < 0.0'
        write ( *, '(a,g14.6)' ) '  x = ', x
        ifault = 1
        return
      end if

      if ( k == n + 1 ) then
        write ( *, '(a)' ) ''
        write ( *, '(a)' ) 'detq - Fatal error!'
        write ( *, '(a)' ) '  k == n + 1'
        ifault = 1
        return
      end if

      do i = k, n
        q = q + n
        x = a2(q) * y
        p = r
        s = q
        do j = k, n
          p = p + 1
          s = s + 1
          a2(s) = a2(s) + x * a2(p)
        end do

      end do

    end if

    r = r + n + 1

  end do

  return
end
subroutine r8mat_print ( m, n, a, title )

!*****************************************************************************80
!
!! r8mat_print() prints an R8MAT.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 August 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer M, the number of rows in A.
!
!    integer N, the number of columns in A.
!
!    real ( kind = rk ) A(M,N), the matrix.
!
!    character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer m
  integer n

  real ( kind = rk ) a(m,n)
  character ( len = * ) title

  call r8mat_print_some ( m, n, a, 1, 1, m, n, title )

  return
end
subroutine r8mat_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )

!*****************************************************************************80
!
!! r8mat_print_some() prints some of an R8MAT.
!
!  Discussion:
!
!    An R8MAT is an MxN array of R8's, stored by (I,J) -> [I+J*M].
!
!  Licensing:
!
!    This code is distributed under the MIT license.
!
!  Modified:
!
!    27 August 2021
!
!  Author:
!
!    John Burkardt
!
!  Input:
!
!    integer M, N, the number of rows and columns.
!
!    real ( kind = rk ) A(M,N), an M by N matrix to be printed.
!
!    integer ILO, JLO, the first row and column to print.
!
!    integer IHI, JHI, the last row and column to print.
!
!    character ( len = * ) TITLE, a title.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer, parameter :: incx = 5
  integer m
  integer n

  real ( kind = rk ) a(m,n)
  character ( len = 14 ) ctemp(incx)
  integer i
  integer i2hi
  integer i2lo
  integer ihi
  integer ilo
  integer inc
  integer j
  integer j2
  integer j2hi
  integer j2lo
  integer jhi
  integer jlo
  character ( len = * ) title

  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) trim ( title )

  if ( m <= 0 .or. n <= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) '  (None)'
    return
  end if

  do j2lo = max ( jlo, 1 ), min ( jhi, n ), incx

    j2hi = j2lo + incx - 1
    j2hi = min ( j2hi, n )
    j2hi = min ( j2hi, jhi )

    inc = j2hi + 1 - j2lo

    write ( *, '(a)' ) ' '

    do j = j2lo, j2hi
      j2 = j + 1 - j2lo
      write ( ctemp(j2), '(i8,6x)' ) j
    end do

    write ( *, '(''  Col   '',5a14)' ) ctemp(1:inc)
    write ( *, '(a)' ) '  Row'
    write ( *, '(a)' ) ' '

    i2lo = max ( ilo, 1 )
    i2hi = min ( ihi, m )

    do i = i2lo, i2hi

      do j2 = 1, inc

        j = j2lo - 1 + j2

        if ( a(i,j) == real ( int ( a(i,j) ), kind = rk ) ) then
          write ( ctemp(j2), '(f8.0,6x)' ) a(i,j)
        else
          write ( ctemp(j2), '(g14.6)' ) a(i,j)
        end if

      end do

      write ( *, '(i5,a,5a14)' ) i, ':', ( ctemp(j), j = 1, inc )

    end do

  end do

  return
end