Add more Lapack routines

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

2 files changed, 147 insertions, 1 deletions

diff --git a/src/lapack.F90 b/src/lapack.F90
index 5bcd70b..cbe41b6 100644
--- a/src/lapack.F90
+++ b/src/lapack.F90
@@ -6,6 +6,10 @@ module lapack
   implicit none
   public
   
+  ! ev: eigenvalues
+  ! ge (general), gg (generalised eigenvalues), sb (symmetric banded),
+  ! sp (symmetric packed), st (symmetric tridiagonal), sy (symmetric)
+  
   interface geev
      subroutine sgeev (jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, &
           &            ldvr, work, lwork, info)
@@ -22,7 +26,28 @@ module lapack
             &           wi(n), work(lwork), wr(n)
      end subroutine dgeev
   end interface
-
+  
+  interface syev
+     subroutine ssyev (jobz, uplo, n, a, lda, w, work, lwork, info)
+       character jobz, uplo
+       integer   info, lda, lwork, n
+       real      a(lda,n), w(n), work(lwork)
+     end subroutine ssyev
+     subroutine dsyev (jobz, uplo, n, a, lda, w, work, lwork, info)
+       character        jobz, uplo
+       integer          info, lda, lwork, n
+       double precision a(lda,n), w(n), work(lwork)
+     end subroutine dsyev
+  end interface
+  
+  
+  
+  ! sv: solve
+  ! 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 gesv
      subroutine sgesv (n, nrhs, a, lda, ipiv, b, ldb, info)
        implicit none
@@ -47,5 +72,61 @@ module lapack
        integer          info
      end subroutine dgesv
   end interface
+
+  interface posv
+     subroutine sposv (uplo, n, nrhs, a, lda, b, ldb, info)
+       implicit none
+       character uplo
+       integer   n
+       integer   nrhs
+       integer   lda
+       real      a(lda,n)
+       integer   ldb
+       real      b(ldb,nrhs)
+       integer   info
+     end subroutine sposv
+     subroutine dposv (uplo, n, nrhs, a, lda, b, ldb, info)
+       implicit none
+       character        uplo
+       integer          n
+       integer          nrhs
+       integer          lda
+       double precision a(lda,n)
+       integer          ldb
+       double precision b(ldb,nrhs)
+       integer          info
+     end subroutine dposv
+  end interface
+
+  interface sysv
+     subroutine ssysv (uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info)
+       implicit none
+       character uplo
+       integer   n
+       integer   nrhs
+       integer   lda
+       real      a(lda,n)
+       integer   ipiv(n)
+       integer   ldb
+       real      b(ldb,nrhs)
+       integer   lwork
+       real      work(lwork)
+       integer   info
+     end subroutine ssysv
+     subroutine dsysv (uplo, n, nrhs, a, lda, ipiv, b, ldb, work, lwork, info)
+       implicit none
+       character        uplo
+       integer          n
+       integer          nrhs
+       integer          lda
+       double precision a(lda,n)
+       integer          ipiv(n)
+       integer          ldb
+       double precision b(ldb,nrhs)
+       integer          lwork
+       double precision work(lwork)
+       integer          info
+     end subroutine dsysv
+  end interface
   
 end module lapack
diff --git a/src/matinv.F90 b/src/matinv.F90
index 4b6d8ec..e6fc6a3 100644
--- a/src/matinv.F90
+++ b/src/matinv.F90
@@ -8,11 +8,53 @@ module matinv
   implicit none
   private
   
+  public calc_posinv3
+  public calc_syminv3
   public calc_inv3
+  
+  public calc_syminv4
   public calc_inv4
   
 contains
   
+  subroutine calc_posinv3 (g3, gu3)
+    CCTK_REAL, intent(in)  :: g3(3,3)
+    CCTK_REAL, intent(out) :: gu3(3,3)
+    CCTK_REAL :: tmp(3,3)
+    integer   :: info
+    character :: msg*1000
+    
+    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 (0, trim(msg))
+    end if
+  end subroutine calc_posinv3
+  
+  subroutine calc_syminv3 (g3, gu3)
+    CCTK_REAL, intent(in)  :: g3(3,3)
+    CCTK_REAL, intent(out) :: gu3(3,3)
+    CCTK_REAL :: tmp(3,3)
+    integer   :: ipiv(3)
+    integer   :: info
+    character :: msg*1000
+    
+    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 (0, trim(msg))
+    end if
+  end subroutine calc_syminv3
+  
   subroutine calc_inv3 (g3, gu3)
     CCTK_REAL, intent(in)  :: g3(3,3)
     CCTK_REAL, intent(out) :: gu3(3,3)
@@ -31,6 +73,29 @@ contains
     end if
   end subroutine calc_inv3
   
+  
+  
+  subroutine calc_syminv4 (g4, gu4)
+    CCTK_REAL, intent(in)  :: g4(0:3,0:3)
+    CCTK_REAL, intent(out) :: gu4(0:3,0:3)
+    CCTK_REAL :: tmp(0:3,0:3)
+    integer   :: ipiv(0:3)
+    integer   :: info
+    character :: msg*1000
+    
+    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 (0, trim(msg))
+    end if
+  end subroutine calc_syminv4
+  
   subroutine calc_inv4 (g4, gu4)
     CCTK_REAL, intent(in)  :: g4(0:3,0:3)
     CCTK_REAL, intent(out) :: gu4(0:3,0:3)