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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
! $Header$
#include "cctk.h"
module matinv
use constants
use lapack
implicit none
private
public calc_inv2
public calc_posinv3
public calc_syminv3
public calc_inv3
public calc_syminv4
public calc_inv4
contains
subroutine calc_inv2 (g2, gu2, lerr)
CCTK_REAL, intent(in) :: g2(2,2)
CCTK_REAL, intent(out) :: gu2(2,2)
logical, optional, intent(out) :: lerr
CCTK_REAL :: tmp(2,2)
integer :: ipiv(2)
integer :: info
character :: msg*100
tmp = g2
gu2 = delta2
call gesv (2, 2, tmp, 2, ipiv, gu2, 2, info)
if (.not. present(lerr) .and. info /= 0) then
write (msg, '("Error in call to GESV, info=",i2)') info
call CCTK_WARN (1, msg)
end if
if (present(lerr)) lerr = info /= 0
end subroutine calc_inv2
subroutine calc_posinv3 (g3, gu3, lerr)
CCTK_REAL, intent(in) :: g3(3,3)
CCTK_REAL, intent(out) :: gu3(3,3)
logical, optional, intent(out) :: lerr
CCTK_REAL :: tmp(3,3)
integer :: info
character :: msg*100
tmp = g3
gu3 = delta3
call posv ('u', 3, 3, tmp, 3, gu3, 3, info)
if (.not. present(lerr) .and. info /= 0) then
write (msg, '("Error in call to POSV, info=",i4)') info
call CCTK_WARN (1, msg)
end if
if (present(lerr)) lerr = info /= 0
end subroutine calc_posinv3
subroutine calc_syminv3 (g3, gu3, lerr)
CCTK_REAL, intent(in) :: g3(3,3)
CCTK_REAL, intent(out) :: gu3(3,3)
logical, optional, intent(out) :: lerr
CCTK_REAL :: tmp(3,3)
integer :: ipiv(3)
integer :: info
character :: msg*100
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 (.not. present(lerr) .and. info /= 0) then
write (msg, '("Error in call to SYSV, info=",i4)') info
call CCTK_WARN (1, msg)
end if
if (present(lerr)) lerr = info /= 0
end subroutine calc_syminv3
subroutine calc_inv3 (g3, gu3, lerr)
CCTK_REAL, intent(in) :: g3(3,3)
CCTK_REAL, intent(out) :: gu3(3,3)
logical, optional, intent(out) :: lerr
CCTK_REAL :: tmp(3,3)
integer :: ipiv(3)
integer :: info
character :: msg*100
tmp = g3
gu3 = delta3
call gesv (3, 3, tmp, 3, ipiv, gu3, 3, info)
if (.not. present(lerr) .and. info /= 0) then
write (msg, '("Error in call to GESV, info=",i2)') info
call CCTK_WARN (1, msg)
end if
if (present(lerr)) lerr = info /= 0
end subroutine calc_inv3
subroutine calc_syminv4 (g4, gu4, lerr)
CCTK_REAL, intent(in) :: g4(0:3,0:3)
CCTK_REAL, intent(out) :: gu4(0:3,0:3)
logical, optional, intent(out) :: lerr
CCTK_REAL :: tmp(0:3,0:3)
integer :: ipiv(0:3)
integer :: info
character :: msg*100
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 (.not. present(lerr) .and. info /= 0) then
write (msg, '("Error in call to SYSV, info=",i4)') info
call CCTK_WARN (1, msg)
end if
if (present(lerr)) lerr = info /= 0
end subroutine calc_syminv4
subroutine calc_inv4 (g4, gu4, lerr)
CCTK_REAL, intent(in) :: g4(0:3,0:3)
CCTK_REAL, intent(out) :: gu4(0:3,0:3)
logical, optional, intent(out) :: lerr
CCTK_REAL :: tmp(0:3,0:3)
integer :: ipiv(0:3)
integer :: info
character :: msg*100
tmp = g4
gu4 = delta4
call gesv (4, 4, tmp, 4, ipiv, gu4, 4, info)
if (.not. present(lerr) .and. info /= 0) then
write (msg, '("Error in call to GESV, info=",i2)') info
call CCTK_WARN (1, msg)
end if
if (present(lerr)) lerr = info /= 0
end subroutine calc_inv4
end module matinv
|
