forked from UK-MAC/TeaLeaf_Trilinos
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathtea_leaf_ppcg.f90
134 lines (105 loc) · 3.91 KB
/
tea_leaf_ppcg.f90
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
MODULE tea_leaf_kernel_ppcg_module
USE tea_leaf_kernel_module
USE tea_leaf_kernel_cheby_module
IMPLICIT NONE
CONTAINS
SUBROUTINE tea_leaf_kernel_ppcg_init_sd(x_min, &
x_max, &
y_min, &
y_max, &
r, &
sd, &
theta )
IMPLICIT NONE
INTEGER(KIND=4):: x_min,x_max,y_min,y_max
REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: r, sd
REAL(KIND=8) :: theta, theta_r
INTEGER :: j,k
theta_r = 1.0_8/theta
!$OMP PARALLEL
!$OMP DO
DO k=y_min,y_max
DO j=x_min,x_max
sd(j, k) = r(j, k)*theta_r
ENDDO
ENDDO
!$OMP END DO
!$OMP END PARALLEL
END SUBROUTINE tea_leaf_kernel_ppcg_init_sd
SUBROUTINE tea_leaf_kernel_ppcg_inner(x_min, &
x_max, &
y_min, &
y_max, &
step, &
alpha, &
beta, &
rx, ry, &
u, &
r, &
Kx, &
Ky, &
sd )
IMPLICIT NONE
INTEGER(KIND=4):: x_min,x_max,y_min,y_max
REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: u, r, Kx, Ky, sd
INTEGER(KIND=4) :: j,k, step
REAL(KIND=8), DIMENSION(:) :: alpha, beta
REAL(KIND=8) :: smvp, rx, ry
!$OMP PARALLEL
!$OMP DO
DO k=y_min,y_max
DO j=x_min,x_max
smvp = (1.0_8 &
+ ry*(Ky(j, k+1) + Ky(j, k)) &
+ rx*(Kx(j+1, k) + Kx(j, k)))*sd(j, k) &
- ry*(Ky(j, k+1)*sd(j, k+1) + Ky(j, k)*sd(j, k-1)) &
- rx*(Kx(j+1, k)*sd(j+1, k) + Kx(j, k)*sd(j-1, k))
r(j, k) = r(j, k) - smvp
u(j, k) = u(j, k) + sd(j, k)
ENDDO
ENDDO
!$OMP END DO
!$OMP DO
DO k=y_min,y_max
DO j=x_min,x_max
sd(j, k) = alpha(step)*sd(j, k) + beta(step)*r(j, k)
ENDDO
ENDDO
!$OMP END DO
!$OMP END PARALLEL
END SUBROUTINE
SUBROUTINE tea_leaf_kernel_ppcg_init_p(x_min, &
x_max, &
y_min, &
y_max, &
p, &
r, &
rro )
IMPLICIT NONE
INTEGER(KIND=4):: x_min,x_max,y_min,y_max
REAL(KIND=8), DIMENSION(x_min-2:x_max+2,y_min-2:y_max+2) :: p, r
INTEGER :: j,k
REAL(KIND=8) :: rro
rro = 0.0_8
!$OMP PARALLEL
!$OMP DO REDUCTION(+:rro)
DO k=y_min,y_max
DO j=x_min,x_max
p(j, k) = r(j, k)
rro = rro + p(j, k)*r(j, k)
ENDDO
ENDDO
!$OMP END DO
!$OMP END PARALLEL
END SUBROUTINE
SUBROUTINE tea_calc_ls_coefs(ch_alphas, ch_betas, eigmin, eigmax, &
theta, ppcg_inner_steps)
INTEGER :: n, ppcg_inner_steps
REAL(KIND=8), DIMENSION(ppcg_inner_steps) :: ch_alphas, ch_betas
REAL(KIND=8) :: eigmin, eigmax
REAL(KIND=8) :: theta, delta, sigma, rho_old, rho_new, cur_alpha, cur_beta
! TODO
CALL tea_calc_ch_coefs(ch_alphas, ch_betas, eigmin, eigmax, &
theta, ppcg_inner_steps)
END SUBROUTINE
END MODULE