subroutine test_invert

  use vector_tools
  use unittest_tools
  implicit none

  real, dimension(10, 10) :: mat, inv, tmp, id
  logical :: fail, warn
  integer :: i
  character(len=20) :: buf

  id = get_matrix_identity(10)


!!$  mat=reshape((/0.0002083333333333334,   -7.9576252343771758e-05, &
!!$        -1.0842021724855044e-19, -7.9576252343771758e-05, &
!!$        -0.00033709041432289469, -1.0842021724855044e-19, &
!!$        -7.9576252343771758e-05, -0.00033709041432289469, &
!!$        -0.00033709041432289469, -1.0842021724855044e-19, &
!!$        -7.9576252343771758e-05, 0.0011516383427084211, &
!!$        -7.9576252343771758e-05, 0.00083333333333333339, &
!!$        0.00083333333333333328 , -0.00033709041432289469, &
!!$        0.00083333333333333339 , 0.00083333333333333328, &
!!$        0.0005150283239582457  , -0.00033709041432289469, &
!!$        -1.0842021724855044e-19, -7.9576252343771758e-05, &
!!$        0.0002083333333333334  , -0.00033709041432289469, &
!!$        -7.9576252343771758e-05, -1.0842021724855044e-19, &
!!$        -0.00033709041432289469, -7.9576252343771758e-05, &
!!$        -0.00033709041432289469, -1.0842021724855044e-19, &
!!$        -7.9576252343771758e-05, 0.00083333333333333339, &
!!$        -0.00033709041432289469, 0.0011516383427084211, &
!!$        0.00083333333333333328 , -7.9576252343771758e-05, &
!!$        0.00083333333333333328 , 0.0005150283239582457, &
!!$        0.00083333333333333328 , -0.00033709041432289469, &
!!$        -0.00033709041432289469, 0.00083333333333333328, &
!!$        -7.9576252343771758e-05, 0.00083333333333333328, &
!!$        0.0011516383427084211  , -7.9576252343771758e-05, &
!!$        0.0005150283239582457  , 0.00083333333333333339, &
!!$        0.00083333333333333328 , -0.00033709041432289469, &
!!$        -1.0842021724855044e-19, -0.00033709041432289469, &
!!$        -1.0842021724855044e-19, -7.9576252343771758e-05, &
!!$        -7.9576252343771758e-05, 0.00020833333333333337, &
!!$        -0.00033709041432289469, -0.00033709041432289469, &
!!$        -7.9576252343771758e-05, -1.0842021724855044e-19, &
!!$        -7.9576252343771758e-05, 0.00083333333333333339, &
!!$        -0.00033709041432289469, 0.00083333333333333328, &
!!$        0.0005150283239582457  , -0.00033709041432289469, &
!!$        0.0011516383427084211  , 0.00083333333333333328, &
!!$        0.00083333333333333328 , -7.9576252343771758e-05, &
!!$        -0.00033709041432289469, 0.00083333333333333328, &
!!$        -7.9576252343771758e-05, 0.0005150283239582457, &
!!$        0.00083333333333333339 , -0.00033709041432289469, &
!!$        0.00083333333333333328 , 0.0011516383427084211, &
!!$        0.00083333333333333328 , -7.9576252343771758e-05, &
!!$        -0.00033709041432289469, 0.0005150283239582457, &
!!$        -0.00033709041432289469, 0.00083333333333333328, &
!!$        0.00083333333333333328 , -7.9576252343771758e-05, &
!!$        0.00083333333333333328 , 0.00083333333333333328, &
!!$        .0011516383427084211   ,-7.9576252343771758e-05, &
!!$        -1.0842021724855044e-19, -0.00033709041432289469, &
!!$        -1.0842021724855044e-19, -0.00033709041432289469, &
!!$        -0.00033709041432289469, -1.0842021724855044e-19, &
!!$        -7.9576252343771758e-05, -7.9576252343771758e-05, &
!!$        -7.9576252343771758e-05, 0.00020833333333333337/),(/10,10/))
!!$
!!$  buf="mm"
!!$  
!!$  inv = mat; call invert(inv)
!!$  tmp = matmul(mat, inv)
!!$  tmp = tmp - id
!!$
!!$  if (.not. mat_zero(tmp)) fail = .true.
!!$  call report_test("[invert " // trim(buf) // "]", fail, warn, "Invert should produce inverses.")

  do i=1,5
    write(buf,'(i0)') i
    fail = .false.; warn = .false.
    mat = random_symmetric_matrix(10)
    inv = mat; call invert(inv)
    tmp = matmul(mat, inv)
    tmp = tmp - id

    if (.not. mat_zero(tmp)) fail = .true.
    call report_test("[invert " // trim(buf) // "]", fail, warn, "Invert should produce inverses.")
  end do

!  mat(1, :) = (/1.0, 5.0, 2.0/)
!  mat(2, :) = (/1.0, 1.0, 7.0/)
!  mat(3, :) = (/0.0,-3.0, 4.0/)
!  inv = mat; call invert(inv)
!  do i=1,3
!    write(0,*) "inv(", i, ", :) == ", inv(i, :)
!  end do

end subroutine test_invert
