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
|
! $Header$
#include "cctk.h"
module gram_schmidt
implicit none
private
public gram_schmidt_project4
public gram_schmidt_normalise4
contains
subroutine gram_schmidt_project4 (g, y, dy, ddy, y2, x, dx, ddx)
CCTK_REAL, intent(in) :: g(0:3,0:3)
CCTK_REAL, intent(in) :: y(0:3), dy(0:3,0:3), ddy(0:3,0:3,0:3), y2
CCTK_REAL, intent(inout) :: x(0:3), dx(0:3,0:3), ddx(0:3,0:3,0:3)
CCTK_REAL :: z(0:3), dz(0:3,0:3), ddz(0:3,0:3,0:3)
integer :: a, b, c, d, e
do a=0,3
z(a) = x(a)
do d=0,3
do e=0,3
z(a) = z(a) - g(d,e) * x(d) * y(e) * y(a) / y2
end do
end do
end do
do a=0,3
do b=0,3
dz(a,b) = dx(a,b)
do d=0,3
do e=0,3
dz(a,b) = dz(a,b) - g(d,e) * ((dx(d,b) * y(e) + x(d) * dy(e,b)) * y(a) + x(d) * y(e) * dy(a,b)) / y2
end do
end do
end do
end do
do a=0,3
do b=0,3
do c=0,3
ddz(a,b,c) = ddx(a,b,c)
do d=0,3
do e=0,3
ddz(a,b,c) = ddz(a,b,c) - g(d,e) * ( (ddx(d,b,c) * y(e) + dx(d,b) * dy(e,c) + dx(d,c) * dy(e,b) + x(d) * ddy(e,b,c)) * y(a) + (dx(d,b) * y(e) + x(d) * dy(e,b)) * dy(a,c) &
& + dx(d,c) * y(e) * dy(a,b) + x(d) * dy(e,c) * dy(a,b) + x(d) * y(e) * ddy(a,b,c)) / y2
end do
end do
end do
end do
end do
x = z
dx = dz
ddx = ddz
end subroutine gram_schmidt_project4
subroutine gram_schmidt_normalise4 (g, x, dx, ddx, x2)
CCTK_REAL, intent(in) :: g(0:3,0:3)
CCTK_REAL, intent(inout) :: x(0:3), dx(0:3,0:3), ddx(0:3,0:3,0:3)
CCTK_REAL, intent(in) :: x2
CCTK_REAL, parameter :: two=2, half=1/two
CCTK_REAL :: z2, dz2(0:3), ddz2(0:3,0:3)
CCTK_REAL :: z(0:3), dz(0:3,0:3), ddz(0:3,0:3,0:3)
integer :: a, b, c, d, e
z2 = 0
do d=0,3
do e=0,3
z2 = z2 + g(d,e) * x(d) * x(e)
end do
end do
do a=0,3
dz2(a) = 0
do d=0,3
do e=0,3
dz2(a) = dz2(a) + 2 * g(d,e) * dx(d,a) * x(e)
end do
end do
end do
do a=0,3
do b=0,3
ddz2(a,b) = 0
do d=0,3
do e=0,3
ddz2(a,b) = ddz2(a,b) + 2 * g(d,e) * ddx(d,a,b) * x(e) + 2 * g(d,e) * dx(d,a) * dx(e,b)
end do
end do
end do
end do
do a=0,3
z(a) = x(a) / sqrt(z2 / x2)
end do
do a=0,3
do b=0,3
dz(a,b) = (dx(a,b) * z2 - half * x(a) * dz2(b) / x2) / sqrt(z2 / x2)**3
end do
end do
do a=0,3
do b=0,3
do c=0,3
ddz(a,b,c) = ((ddx(a,b,c) * z2 + dx(a,b) * dz2(c) / x2 - half * (dx(a,c) * dz2(b) / x2 + x(a) * ddz2(b,c) / x2**2)) * z2 &
& - 3*half * (dx(a,b) * z2 - half * x(a) * dz2(b) / x2) * dz2(c) / x2) / sqrt(z2 / x2)**5
end do
end do
end do
x = z
dx = dz
ddx = ddz
end subroutine gram_schmidt_normalise4
end module gram_schmidt
|
