Correct routines that calculate the determinant.

Add routine for dim=3.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinUtils/TGRtensor/trunk@19 b716e942-a2de-43ad-8f52-f3dfe468e4e7

1 files changed, 87 insertions, 6 deletions

diff --git a/src/matdet.F90 b/src/matdet.F90
index 3ac9226..da3d673 100644
--- a/src/matdet.F90
+++ b/src/matdet.F90
@@ -7,10 +7,93 @@ module matdet
   implicit none
   private
   
+  public calc_symdet3
+  public calc_det3
+  
   public calc_symdet4
   
 contains
   
+  subroutine calc_symdet3 (gg, dtg, lerr)
+    CCTK_REAL,           intent(in)  :: gg(3,3)
+    CCTK_REAL,           intent(out) :: dtg
+    logical,   optional, intent(out) :: lerr
+    CCTK_REAL :: tmp(3,3)
+    integer   :: ipiv(3)
+    integer   :: info
+    logical   :: odd
+    integer   :: i
+    character :: msg*100
+    
+    tmp = gg
+    call sytrf (3, 3, tmp, 3, ipiv, info)
+    
+    if (info < 0) then
+       write (msg, '("Error in call to SYTRF, info=",i4)') info
+       call CCTK_WARN (1, msg)
+    end if
+    
+    if (present(lerr)) lerr = info /= 0
+    
+    if (info > 0) then
+       dtg = 0
+       return
+    end if
+    
+    odd = .false.
+    do i=1,3
+       if (ipiv(i) /= i) odd = .not. odd
+    end do
+    
+    dtg = 1
+    if (odd) dtg = -dtg
+    do i=1,3
+       dtg = dtg * tmp(i,i)
+    end do
+  end subroutine calc_symdet3
+  
+  
+  
+  subroutine calc_det3 (gg, dtg, lerr)
+    CCTK_REAL,           intent(in)  :: gg(3,3)
+    CCTK_REAL,           intent(out) :: dtg
+    logical,   optional, intent(out) :: lerr
+    CCTK_REAL :: tmp(3,3)
+    integer   :: ipiv(3)
+    integer   :: info
+    logical   :: odd
+    integer   :: i
+    character :: msg*100
+    
+    tmp = gg
+    call getrf (3, 3, tmp, 3, ipiv, info)
+    
+    if (info < 0) then
+       write (msg, '("Error in call to GETRF, info=",i4)') info
+       call CCTK_WARN (1, msg)
+    end if
+    
+    if (present(lerr)) lerr = info /= 0
+    
+    if (info > 0) then
+       dtg = 0
+       return
+    end if
+    
+    odd = .false.
+    do i=1,3
+       if (ipiv(i) /= i) odd = .not. odd
+    end do
+    
+    dtg = 1
+    if (odd) dtg = -dtg
+    do i=1,3
+       dtg = dtg * tmp(i,i)
+    end do
+  end subroutine calc_det3
+  
+  
+  
   subroutine calc_symdet4 (g4, dtg4, lerr)
     CCTK_REAL,           intent(in)  :: g4(4,4)
     CCTK_REAL,           intent(out) :: dtg4
@@ -18,7 +101,7 @@ contains
     CCTK_REAL :: tmp(4,4)
     integer   :: ipiv(4)
     integer   :: info
-    integer   :: nperms
+    logical   :: odd
     integer   :: i
     character :: msg*100
     
@@ -37,15 +120,13 @@ contains
        return
     end if
     
-    nperms = 0
+    odd = .false.
     do i=1,4
-       if (mod(ipiv(i),2) /= mod(i,2)) nperms = nperms+1
+       if (ipiv(i) /= i) odd = .not. odd
     end do
-    if (mod(nperms,2) /=0) call CCTK_WARN (0, "internal error")
-    nperms = nperms / 2
     
     dtg4 = 1
-    if (mod(nperms,2)/=0) dtg4 = -1
+    if (odd) dtg4 = -dtg4
     do i=1,4
        dtg4 = dtg4 * tmp(i,i)
     end do