Add routines for calculating the determinant of a 4d matrix.

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

3 files changed, 81 insertions, 0 deletions

diff --git a/src/lapack.F90 b/src/lapack.F90
index cbe41b6..58e64e1 100644
--- a/src/lapack.F90
+++ b/src/lapack.F90
@@ -129,4 +129,33 @@ module lapack
      end subroutine dsysv
   end interface
   
+  
+  
+  ! trf: LU factorisation
+  ! gb (general banded), ge (general), gt (general tridiagonal),
+  ! pb (symmetric positive banded), po (symmetric positive),
+  ! pp (symmetric positive packed), pt (symmetric positive tridiagonal),
+  ! sp (symmetric packed), sy (symmetric)
+  
+  interface sytrf
+     subroutine ssytrf (m, n, a, lda, ipiv, info)
+       implicit none
+       integer m
+       integer n
+       integer lda
+       real    a(lda,n)
+       integer ipiv(n)
+       integer info
+     end subroutine ssytrf
+     subroutine dsytrf (m, n, a, lda, ipiv, info)
+       implicit none
+       integer          m
+       integer          n
+       integer          lda
+       double precision a(lda,n)
+       integer          ipiv(n)
+       integer          info
+     end subroutine dsytrf
+  end interface
+  
 end module lapack
diff --git a/src/make.code.defn b/src/make.code.defn
index b589337..69e3078 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -10,6 +10,7 @@ SRCS =	cactus.F90				\
 	derivs.F90				\
 	derivs2.F90				\
 	lapack.F90				\
+	matdet.F90				\
 	matexp.F90				\
 	matinv.F90				\
 	rotation.F90				\
diff --git a/src/matdet.F90 b/src/matdet.F90
new file mode 100644
index 0000000..2aaa754
--- /dev/null
+++ b/src/matdet.F90
@@ -0,0 +1,51 @@
+! $Header$
+
+#include "cctk.h"
+
+module matdet
+  use lapack
+  implicit none
+  private
+  
+  public calc_symdet4
+  
+contains
+  
+  subroutine calc_symdet4 (g4, dtg4)
+    CCTK_REAL, intent(in)  :: g4(4,4)
+    CCTK_REAL, intent(out) :: dtg4
+    CCTK_REAL :: tmp(4,4)
+    integer   :: ipiv(4)
+    integer   :: info
+    integer   :: nperms
+    integer   :: i
+    character :: msg*1000
+    
+    tmp = g4
+    call sytrf (4, 4, tmp, 4, ipiv, info)
+    
+    if (info < 0) then
+       write (msg, '("Error in call to SYTRF, info=",i4)') info
+       call CCTK_WARN (0, trim(msg))
+    end if
+    
+    if (info > 0) then
+       dtg4 = 0
+       return
+    end if
+    
+    nperms = 0
+    do i=1,4
+       if (mod(ipiv(i),2) /= mod(i,2)) nperms = nperms+1
+    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
+    do i=1,4
+       dtg4 = dtg4 * tmp(i,i)
+    end do
+  end subroutine calc_symdet4
+  
+end module matdet