FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_solver_BiCGSTAB.f90
Go to the documentation of this file.
1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
6!C
7!C***
8!C*** module hecmw_solver_BiCGSTAB
9!C***
10!C
12contains
13 !C
14 !C*** BiCGSTAB_3
15 !C
16 subroutine hecmw_solve_bicgstab( hecMESH, hecMAT, ITER, RESID, error, &
17 & Tset, Tsol, Tcomm )
18
19 use hecmw_util
27
28 implicit none
29
30 type(hecmwst_local_mesh) :: hecMESH
31 type(hecmwst_matrix) :: hecMAT
32 integer(kind=kint ), intent(inout):: ITER, error
33 real (kind=kreal), intent(inout):: resid, tset, tsol, tcomm
34
35 integer(kind=kint ) :: N, NP, NDOF, NNDOF
36 integer(kind=kint ) :: my_rank
37 integer(kind=kint ) :: ITERlog, TIMElog
38 real(kind=kreal), pointer :: b(:), x(:)
39
40 real(kind=kreal), dimension(:,:), allocatable :: ww
41 real(kind=kreal), dimension(2) :: cg
42
43 integer(kind=kint ) :: MAXIT
44
45 ! local variables
46 real (kind=kreal):: tol
47 integer(kind=kint )::i
48 real (kind=kreal)::s_time,s1_time,e_time,e1_time, start_time, end_time
49 real (kind=kreal)::bnrm2,c2
50 real (kind=kreal)::rho,rho1,beta,alpha,dnrm2
51 real (kind=kreal)::omega
52 real (kind=kreal)::t_max,t_min,t_avg,t_sd
53
54 integer(kind=kint), parameter :: R = 1
55 integer(kind=kint), parameter :: RT= 2
56 integer(kind=kint), parameter :: P = 3
57 integer(kind=kint), parameter :: PT= 4
58 integer(kind=kint), parameter :: S = 5
59 integer(kind=kint), parameter :: ST= 1
60 integer(kind=kint), parameter :: T = 6
61 integer(kind=kint), parameter :: V = 7
62 integer(kind=kint), parameter :: WK= 8
63
64 integer(kind=kint), parameter :: N_ITER_RECOMPUTE_R= 100
65
66 call hecmw_barrier(hecmesh)
67 s_time= hecmw_wtime()
68
69 !C===
70 !C +-------+
71 !C | INIT. |
72 !C +-------+
73 !C===
74 n = hecmat%N
75 np = hecmat%NP
76 ndof = hecmat%NDOF
77 nndof = n * ndof
78 my_rank = hecmesh%my_rank
79 x => hecmat%X
80 b => hecmat%B
81
82 iterlog = hecmw_mat_get_iterlog( hecmat )
83 timelog = hecmw_mat_get_timelog( hecmat )
84 maxit = hecmw_mat_get_iter( hecmat )
85 tol = hecmw_mat_get_resid( hecmat )
86
87 error = 0
88 rho1 = 0.0d0
89 alpha = 0.0d0
90 beta = 0.0d0
91 omega = 0.0d0
92
93 allocate (ww(ndof*np, 8))
94 ww = 0.d0
95
96 !C
97 !C-- SCALING
98 call hecmw_solver_scaling_fw(hecmesh, hecmat, tcomm)
99
100 !C===
101 !C +----------------------+
102 !C | SETUP PRECONDITIONER |
103 !C +----------------------+
104 !C===
105 call hecmw_precond_setup(hecmat, hecmesh, 0)
106
107 !C===
108 !C +---------------------+
109 !C | {r0}= {b} - [A]{x0} |
110 !C +---------------------+
111 !C===
112 call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
113
114 !C-- set arbitrary {r_tld}
115 do i=1, nndof
116 ww(i,rt) = ww(i,r)
117 enddo
118
119 !C-- compute ||{b}||
120 call hecmw_innerproduct_r(hecmesh, ndof, b, b, bnrm2, tcomm)
121 if (bnrm2.eq.0.d0) then
122 iter = 0
123 maxit = 0
124 resid = 0.d0
125 x = 0.d0
126 endif
127
128 e_time = hecmw_wtime()
129 if (timelog.eq.2) then
130 call hecmw_time_statistics(hecmesh, e_time - s_time, &
131 t_max, t_min, t_avg, t_sd)
132 if (hecmesh%my_rank.eq.0) then
133 write(*,*) 'Time solver setup'
134 write(*,*) ' Max :',t_max
135 write(*,*) ' Min :',t_min
136 write(*,*) ' Avg :',t_avg
137 write(*,*) ' Std Dev :',t_sd
138 endif
139 tset = t_max
140 else
141 tset = e_time - s_time
142 endif
143
144 tcomm = 0.d0
145 call hecmw_barrier(hecmesh)
146 s1_time = hecmw_wtime()
147 !C
148 !C*************************************************************** iterative procedures start
149 !C
150 do iter = 1, maxit
151
152 !C===
153 !C +-----------------+
154 !C | RHO= {r}{r_tld} |
155 !C +-----------------+
156 !C===
157 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,rt), rho, tcomm)
158
159 !C===
160 !C +----------------------------------------+
161 !C | BETA= (RHO/RHO1) * (ALPHA/OMEGA) |
162 !C | {p} = {r} + BETA * ( {p} - OMEGA*{v} ) |
163 !C +----------------------------------------+
164 !C===
165 if ( iter.gt.1 ) then
166 beta = (rho/rho1) * (alpha/omega)
167 do i = 1, nndof
168 ww(i,p) = ww(i,r) + beta * (ww(i,p) - omega * ww(i,v))
169 enddo
170 else
171 do i = 1, nndof
172 ww(i,p) = ww(i,r)
173 enddo
174 endif
175
176 !C===
177 !C +--------------------+
178 !C | {p_tld}= [Minv]{p} |
179 !C +--------------------+
180 !C===
181 call hecmw_precond_apply(hecmesh, hecmat, ww(:, p), ww(:, pt), ww(:, wk), tcomm)
182
183 !C===
184 !C +-------------------------+
185 !C | {v}= [A] {p_tld} |
186 !C +-------------------------+
187 !C===
188 call hecmw_matvec(hecmesh, hecmat, ww(:,pt), ww(:,v), tcomm)
189
190 !C
191 !C-- calc. ALPHA
192 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,rt), ww(:,v), c2, tcomm)
193
194 alpha = rho / c2
195
196 !C
197 !C-- {s}= {r} - ALPHA*{V}
198 do i = 1, nndof
199 ww(i,s) = ww(i,r) - alpha * ww(i,v)
200 enddo
201
202 !C===
203 !C +--------------------+
204 !C | {s_tld}= [Minv]{s} |
205 !C +--------------------+
206 !C===
207 call hecmw_precond_apply(hecmesh, hecmat, ww(:, s), ww(:, st), ww(:, wk), tcomm)
208
209 !C===
210 !C +-------------------------+
211 !C | {t}= [A] {s_tld} |
212 !C +-------------------------+
213 !C===
214 call hecmw_matvec(hecmesh, hecmat, ww(:,st), ww(:,t), tcomm)
215
216 !C===
217 !C +----------------------------+
218 !C | OMEGA= ({t}{s}) / ({t}{t}) |
219 !C +----------------------------+
220 !C===
221 call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,t), ww(:,s), cg(1))
222 call hecmw_innerproduct_r_nocomm(hecmesh, ndof, ww(:,t), ww(:,t), cg(2))
223 s_time= hecmw_wtime()
224 call hecmw_allreduce_r(hecmesh, cg, 2, hecmw_sum)
225 e_time= hecmw_wtime()
226 tcomm = tcomm + e_time - s_time
227
228 omega = cg(1) / cg(2)
229
230 !C===
231 !C +----------------+
232 !C | update {x},{r} |
233 !C +----------------+
234 !C===
235 do i = 1, nndof
236 x(i) = x(i) + alpha * ww(i,pt) + omega * ww(i,st)
237 enddo
238 !C
239 !C--- recompute R sometimes
240 if ( mod(iter,n_iter_recompute_r)==0 ) then
241 call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
242 else
243 do i = 1, nndof
244 ww(i,r) = ww(i,s) - omega * ww(i,t)
245 enddo
246 endif
247
248 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
249
250 resid= dsqrt(dnrm2/bnrm2)
251
252 !C##### ITERATION HISTORY
253 if (my_rank.eq.0.and.iterlog.eq.1) write (*,'(i7, 1pe16.6)') iter, resid
254 !C#####
255
256 if ( resid.le.tol ) then
257 if ( mod(iter,n_iter_recompute_r)==0 ) exit
258 !C----- recompute R to make sure it is really converged
259 call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
260 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
261 resid= dsqrt(dnrm2/bnrm2)
262 if ( resid.le.tol ) exit
263 endif
264 if ( iter .eq.maxit ) error = hecmw_solver_error_noconv_maxit
265
266 rho1 = rho
267
268 enddo
269 !C
270 !C*************************************************************** iterative procedures end
271 !C
272
273 call hecmw_solver_scaling_bk(hecmat)
274 !C
275 !C-- INTERFACE data EXCHANGE
276 !C
277 start_time = hecmw_wtime()
278 call hecmw_update_m_r (hecmesh, x, hecmat%NP, hecmat%NDOF)
279 end_time = hecmw_wtime()
280 tcomm = tcomm + end_time - start_time
281
282 deallocate (ww)
283 !call hecmw_precond_clear(hecMAT)
284
285 e1_time = hecmw_wtime()
286 if (timelog.eq.2) then
287 call hecmw_time_statistics(hecmesh, e1_time - s1_time, &
288 t_max, t_min, t_avg, t_sd)
289 if (hecmesh%my_rank.eq.0) then
290 write(*,*) 'Time solver iterations'
291 write(*,*) ' Max :',t_max
292 write(*,*) ' Min :',t_min
293 write(*,*) ' Avg :',t_avg
294 write(*,*) ' Std Dev :',t_sd
295 endif
296 tsol = t_max
297 else
298 tsol = e1_time - s1_time
299 endif
300
301 end subroutine hecmw_solve_bicgstab
302end module hecmw_solver_bicgstab
integer(kind=kint) function, public hecmw_mat_get_timelog(hecmat)
real(kind=kreal) function, public hecmw_mat_get_resid(hecmat)
integer(kind=kint) function, public hecmw_mat_get_iter(hecmat)
integer(kind=kint) function, public hecmw_mat_get_iterlog(hecmat)
subroutine, public hecmw_precond_setup(hecmat, hecmesh, sym)
subroutine, public hecmw_precond_apply(hecmesh, hecmat, r, z, zp, commtime)
subroutine hecmw_solve_bicgstab(hecmesh, hecmat, iter, resid, error, tset, tsol, tcomm)
subroutine, public hecmw_matresid(hecmesh, hecmat, x, b, r, commtime)
subroutine, public hecmw_matvec(hecmesh, hecmat, x, y, commtime)
subroutine hecmw_innerproduct_r_nocomm(hecmesh, ndof, x, y, sum)
subroutine hecmw_time_statistics(hecmesh, time, t_max, t_min, t_avg, t_sd)
subroutine hecmw_innerproduct_r(hecmesh, ndof, x, y, sum, commtime)
subroutine, public hecmw_solver_scaling_fw(hecmesh, hecmat, commtime)
subroutine, public hecmw_solver_scaling_bk(hecmat)
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=kint), parameter hecmw_sum
integer(kind=4), parameter kreal
real(kind=kreal) function hecmw_wtime()
subroutine hecmw_barrier(hecmesh)
subroutine hecmw_update_m_r(hecmesh, val, n, m)
subroutine hecmw_allreduce_r(hecmesh, val, n, ntag)
integer(kind=kint), parameter hecmw_solver_error_noconv_maxit