1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
|
! $Header$
#include "cctk.h"
module matdet
use lapack
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
logical, optional, intent(out) :: lerr
CCTK_REAL :: tmp(4,4)
integer :: ipiv(4)
integer :: info
logical :: odd
integer :: i
character :: msg*100
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 (1, msg)
end if
if (present(lerr)) lerr = info /= 0
if (info > 0) then
dtg4 = 0
return
end if
odd = .false.
do i=1,4
if (ipiv(i) /= i) odd = .not. odd
end do
dtg4 = 1
if (odd) dtg4 = -dtg4
do i=1,4
dtg4 = dtg4 * tmp(i,i)
end do
end subroutine calc_symdet4
end module matdet
|
