From 1496881ca2aaf851b69b4112cd3c1397aa5ee657 Mon Sep 17 00:00:00 2001
From: schnetter <schnetter@b716e942-a2de-43ad-8f52-f3dfe468e4e7>
Date: Sun, 27 Feb 2005 13:06:14 +0000
Subject: 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
---
 src/matdet.F90 | 93 ++++++++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file 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
-- 
cgit v1.2.3