! $Header$

#include "cctk.h"

module matinv
  use constants
  use lapack
  implicit none
  private
  
  public calc_posinv3
  public calc_syminv3
  public calc_inv3
  
  public calc_syminv4
  public calc_inv4
  
contains
  
  subroutine calc_posinv3 (g3, gu3, lerr)
    CCTK_REAL,           intent(in)  :: g3(3,3)
    CCTK_REAL,           intent(out) :: gu3(3,3)
    logical,   optional, intent(out) :: lerr
    CCTK_REAL :: tmp(3,3)
    integer   :: info
    character :: msg*100
    
    tmp = g3
    gu3 = delta3
    call posv ('u', 3, 3, tmp, 3, gu3, 3, info)
    
    if (info /= 0) then
       write (msg, '("Error in call to POSV, info=",i4)') info
       call CCTK_WARN (1, msg)
    end if
    
    if (present(lerr)) lerr = info /= 0
  end subroutine calc_posinv3
  
  subroutine calc_syminv3 (g3, gu3, lerr)
    CCTK_REAL,           intent(in)  :: g3(3,3)
    CCTK_REAL,           intent(out) :: gu3(3,3)
    logical,   optional, intent(out) :: lerr
    CCTK_REAL :: tmp(3,3)
    integer   :: ipiv(3)
    integer   :: info
    character :: msg*100
    
    integer, parameter :: lwork = 1000
    CCTK_REAL :: work(lwork)
    
    tmp = g3
    gu3 = delta3
    call sysv ('u', 3, 3, tmp, 3, ipiv, gu3, 3, work, lwork, info)
    
    if (info /= 0) then
       write (msg, '("Error in call to SYSV, info=",i4)') info
       call CCTK_WARN (1, msg)
    end if
    
    if (present(lerr)) lerr = info /= 0
  end subroutine calc_syminv3
  
  subroutine calc_inv3 (g3, gu3, lerr)
    CCTK_REAL,           intent(in)  :: g3(3,3)
    CCTK_REAL,           intent(out) :: gu3(3,3)
    logical,   optional, intent(out) :: lerr
    CCTK_REAL :: tmp(3,3)
    integer   :: ipiv(3)
    integer   :: info
    character :: msg*100
    
    tmp = g3
    gu3 = delta3
    call gesv (3, 3, tmp, 3, ipiv, gu3, 3, info)
    
    if (info /= 0) then
       write (msg, '("Error in call to GESV, info=",i2)') info
       call CCTK_WARN (1, msg)
    end if
    
    if (present(lerr)) lerr = info /= 0
  end subroutine calc_inv3
  
  
  
  subroutine calc_syminv4 (g4, gu4, lerr)
    CCTK_REAL,           intent(in)  :: g4(0:3,0:3)
    CCTK_REAL,           intent(out) :: gu4(0:3,0:3)
    logical,   optional, intent(out) :: lerr
    CCTK_REAL :: tmp(0:3,0:3)
    integer   :: ipiv(0:3)
    integer   :: info
    character :: msg*100
    
    integer, parameter :: lwork = 1000
    CCTK_REAL :: work(lwork)
    
    tmp = g4
    gu4 = delta4
    call sysv ('u', 4, 4, tmp, 4, ipiv, gu4, 4, work, lwork, info)
    
    if (info /= 0) then
       write (msg, '("Error in call to SYSV, info=",i4)') info
       call CCTK_WARN (1, msg)
    end if
    
    if (present(lerr)) lerr = info /= 0
  end subroutine calc_syminv4
  
  subroutine calc_inv4 (g4, gu4, lerr)
    CCTK_REAL,           intent(in)  :: g4(0:3,0:3)
    CCTK_REAL,           intent(out) :: gu4(0:3,0:3)
    logical,   optional, intent(out) :: lerr
    CCTK_REAL :: tmp(0:3,0:3)
    integer   :: ipiv(0:3)
    integer   :: info
    character :: msg*100
    
    tmp = g4
    gu4 = delta4
    call gesv (4, 4, tmp, 4, ipiv, gu4, 4, info)
    
    if (info /= 0) then
       write (msg, '("Error in call to GESV, info=",i2)') info
       call CCTK_WARN (1, msg)
    end if
    
    if (present(lerr)) lerr = info /= 0
  end subroutine calc_inv4
  
end module matinv