FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_solver_GMRES.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*** module hecmw_solver_GMRES
8!C***
9!
11
12 public :: hecmw_solve_gmres
13
14contains
15 !C
16 !C*** hecmw_solve_GMRES
17 !C
18 subroutine hecmw_solve_gmres( hecMESH, hecMAT, ITER, RESID, error, &
19 & Tset, Tsol, Tcomm )
20 use hecmw_util
29
30 implicit none
31
32 type(hecmwst_local_mesh) :: hecmesh
33 type(hecmwst_matrix) :: hecmat
34 integer(kind=kint ), intent(inout):: iter, error
35 real (kind=kreal), intent(inout):: resid, tset, tsol, tcomm
36
37 integer(kind=kint ) :: n, np, ndof, nndof
38 integer(kind=kint ) :: my_rank
39 integer(kind=kint ) :: iterlog, timelog
40 real(kind=kreal), pointer :: b(:), x(:)
41
42 real(kind=kreal), dimension(:,:), allocatable :: ww
43
44 integer(kind=kint ) :: maxit, nrest
45
46 real (kind=kreal) :: tol
47
48 real (kind=kreal), dimension(:), allocatable :: ss
49 real (kind=kreal), dimension(:,:), allocatable :: h
50
51 integer(kind=kint ) :: cs, sn
52
53 real (kind=kreal) zero, one
54 parameter( zero = 0.0d+0, one = 1.0d+0 )
55
56 integer(kind=kint ) :: nrk,i,k,kk,jj,info,ik
57 integer(kind=kint ) :: irow
58 real (kind=kreal) :: s_time,e_time,s1_time,e1_time
59 real (kind=kreal) :: ldh,ldw,bnrm2,dnrm2,rnorm
60 real (kind=kreal) :: commtime,comptime, coef,val,vcs,vsn,dtemp,aa,bb,r0,scale,rr
61 integer(kind=kint ) :: estcond
62 real (kind=kreal) :: t_max,t_min,t_avg,t_sd
63
64 integer(kind=kint), parameter :: r = 1
65 integer(kind=kint), parameter :: zp = r + 1
66 integer(kind=kint), parameter :: zq = r + 2
67 integer(kind=kint), parameter :: s = r + 3
68 integer(kind=kint), parameter :: w = s + 1
69 integer(kind=kint), parameter :: y = w
70 integer(kind=kint), parameter :: av = y + 1
71 integer(kind=kint), parameter :: v = av + 1
72
73 call hecmw_barrier(hecmesh)
74 s_time= hecmw_wtime()
75 !C
76 !C-- INIT.
77 n = hecmat%N
78 np = hecmat%NP
79 ndof = hecmat%NDOF
80 nndof = n * ndof
81 my_rank = hecmesh%my_rank
82 x => hecmat%X
83 b => hecmat%B
84
85 iterlog = hecmw_mat_get_iterlog( hecmat )
86 timelog = hecmw_mat_get_timelog( hecmat )
87 maxit = hecmw_mat_get_iter( hecmat )
88 tol = hecmw_mat_get_resid( hecmat )
89 nrest = hecmw_mat_get_nrest( hecmat )
90 estcond = hecmw_mat_get_estcond( hecmat )
91
92 if (nrest >= ndof*np-1) nrest = ndof*np-2
93
94 error= 0
95 nrk= nrest + 7
96
97 allocate (h(nrk,nrk))
98 allocate (ww(ndof*np,nrk))
99 allocate (ss(nrk))
100
101 commtime= 0.d0
102 comptime= 0.d0
103
104 ldh= nrest + 2
105 ldw= n
106
107 !C
108 !C-- Store the Givens parameters in matrix H.
109 cs= nrest + 1
110 sn= cs + 1
111
112 !C
113 !C-- SCALING
114 call hecmw_solver_scaling_fw(hecmesh, hecmat, tcomm)
115
116 !C===
117 !C +----------------------+
118 !C | SETUP PRECONDITIONER |
119 !C +----------------------+
120 !C===
121 call hecmw_precond_setup(hecmat, hecmesh, 0)
122
123 !C
124 !C
125 !C +--------------------+
126 !C | {r}= {b} - [A]{x0} |
127 !C +--------------------+
128 !C===
129 call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
130 !C===
131
132 call hecmw_innerproduct_r(hecmesh, ndof, b, b, bnrm2, tcomm)
133 if (bnrm2.eq.0.d0) then
134 iter = 0
135 maxit = 0
136 resid = 0.d0
137 x = 0.d0
138 endif
139
140 e_time= hecmw_wtime()
141 if (timelog.eq.2) then
142 call hecmw_time_statistics(hecmesh, e_time - s_time, &
143 t_max, t_min, t_avg, t_sd)
144 if (hecmesh%my_rank.eq.0) then
145 write(*,*) 'Time solver setup'
146 write(*,*) ' Max :',t_max
147 write(*,*) ' Min :',t_min
148 write(*,*) ' Avg :',t_avg
149 write(*,*) ' Std Dev :',t_sd
150 endif
151 tset = t_max
152 else
153 tset = e_time - s_time
154 endif
155 !C===
156
157 call hecmw_barrier(hecmesh)
158 s1_time= hecmw_wtime()
159 iter= 0
160
161 outer: do
162
163 !C
164 !C************************************************ GMRES Iteration
165 !C
166 i= 0
167 !C
168 !C +---------------+
169 !C | {v1}= {r}/|r| |
170 !C +---------------+
171 !C===
172 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
173 if (dnrm2 == 0.d0) exit ! converged
174
175 rnorm= dsqrt(dnrm2)
176 coef= one/rnorm
177 do ik= 1, nndof
178 ww(ik,v)= ww(ik,r) * coef
179 enddo
180 !C===
181
182 !C
183 !C +--------------+
184 !C | {s}= |r|{e1} |
185 !C +--------------+
186 !C===
187 ww(1 ,s) = rnorm
188 do k = 2, nndof
189 ww(k,s) = zero
190 enddo
191 !C===
192
193 !C************************************************ GMRES(m) restart
194 do i = 1, nrest
195 iter= iter + 1
196
197 !C
198 !C +-------------------+
199 !C | {w}= [A][Minv]{v} |
200 !C +-------------------+
201 !C===
202 call hecmw_precond_apply(hecmesh, hecmat, ww(:,v+i-1), ww(:,zq), ww(:,zp), tcomm)
203
204 call hecmw_matvec(hecmesh, hecmat, ww(:,zq), ww(:,w), tcomm)
205 !C===
206
207 !C
208 !C +------------------------------+
209 !C | ORTH. BASIS by GRAMM-SCHMIDT |
210 !C +------------------------------+
211 !C Construct the I-th column of the upper Hessenberg matrix H
212 !C using the Gram-Schmidt process on V and W.
213 !C===
214 do k= 1, i
215 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,w), ww(:,v+k-1), val, tcomm)
216
217 do ik= 1, nndof
218 ww(ik,w)= ww(ik,w) - val * ww(ik,v+k-1)
219 enddo
220 h(k,i)= val
221 enddo
222
223 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,w), ww(:,w), val, tcomm)
224 if (val == 0.d0) exit ! converged
225
226 h(i+1,i)= dsqrt(val)
227 coef= one / h(i+1,i)
228 do ik= 1, nndof
229 ww(ik,v+i+1-1)= ww(ik,w) * coef
230 enddo
231 !C===
232
233 !C
234 !C +-----------------+
235 !C | GIVENS ROTARION |
236 !C +-----------------+
237 !C===
238
239 !C
240 !C-- Plane Rotation
241 do k = 1, i-1
242 vcs= h(k,cs)
243 vsn= h(k,sn)
244 dtemp = vcs*h(k ,i) + vsn*h(k+1,i)
245 h(k+1,i)= vcs*h(k+1,i) - vsn*h(k ,i)
246 h(k ,i)= dtemp
247 enddo
248
249 !C
250 !C-- Construct Givens Plane Rotation
251 aa = h(i ,i)
252 bb = h(i+1,i)
253 r0= bb
254 if (dabs(aa).gt.dabs(bb)) r0= aa
255 scale= dabs(aa) + dabs(bb)
256
257 if (scale.ne.0.d0) then
258 rr= scale * dsqrt((aa/scale)**2+(bb/scale)**2)
259 rr= dsign(1.d0,r0)*rr
260 h(i,cs)= aa/rr
261 h(i,sn)= bb/rr
262 else
263 h(i,cs)= 1.d0
264 h(i,sn)= 0.d0
265 rr = 0.d0
266 endif
267
268 !C
269 !C-- Plane Rotation
270 vcs= h(i,cs)
271 vsn= h(i,sn)
272 dtemp = vcs*h(i ,i) + vsn*h(i+1,i)
273 h(i+1,i)= vcs*h(i+1,i) - vsn*h(i ,i)
274 h(i ,i)= dtemp
275
276 dtemp = vcs*ww(i ,s) + vsn*ww(i+1,s)
277 ww(i+1,s)= vcs*ww(i+1,s) - vsn*ww(i ,s)
278 ww(i ,s)= dtemp
279
280 resid = dabs( ww(i+1,s))/dsqrt(bnrm2)
281
282 if (my_rank.eq.0 .and. iterlog.eq.1) &
283 & write (*, '(2i8, 1pe16.6)') iter,i+1, resid
284
285 if (estcond /= 0 .and. hecmesh%my_rank == 0) then
286 if (mod(iter,estcond) == 0) call hecmw_estimate_condition_gmres(i, h)
287 endif
288
289 if ( resid.le.tol ) then
290 !C-- [H]{y}= {s_tld}
291 do ik= 1, i
292 ss(ik)= ww(ik,s)
293 enddo
294 irow= i
295 ww(irow,y)= ss(irow) / h(irow,irow)
296
297 do kk= irow-1, 1, -1
298 do jj= irow, kk+1, -1
299 ss(kk)= ss(kk) - h(kk,jj)*ww(jj,y)
300 enddo
301 ww(kk,y)= ss(kk) / h(kk,kk)
302 enddo
303
304 !C-- {x}= {x} + {y}{V}
305 do kk= 1, nndof
306 ww(kk, av)= 0.d0
307 enddo
308
309 jj= irow
310 do jj= 1, irow
311 do kk= 1, nndof
312 ww(kk,av)= ww(kk,av) + ww(jj,y)*ww(kk,v+jj-1)
313 enddo
314 enddo
315
316 call hecmw_precond_apply(hecmesh, hecmat, ww(:,av), ww(:,zq), ww(:,zp), tcomm)
317
318 do kk= 1, nndof
319 x(kk)= x(kk) + ww(kk,zq)
320 enddo
321
322 exit outer
323 endif
324
325 if ( iter.gt.maxit ) then
327 exit outer
328 end if
329 end do
330 !C===
331
332 !C
333 !C +------------------+
334 !C | CURRENT SOLUTION |
335 !C +------------------+
336 !C===
337
338 !C-- [H]{y}= {s_tld}
339 do ik= 1, nrest
340 ss(ik)= ww(ik,s)
341 enddo
342 irow= nrest
343 ww(irow,y)= ss(irow) / h(irow,irow)
344
345 do kk= irow-1, 1, -1
346 do jj= irow, kk+1, -1
347 ss(kk)= ss(kk) - h(kk,jj)*ww(jj,y)
348 enddo
349 ww(kk,y)= ss(kk) / h(kk,kk)
350 enddo
351
352 !C-- {x}= {x} + {y}{V}
353 do kk= 1, nndof
354 ww(kk, av)= 0.d0
355 enddo
356
357 jj= irow
358 do jj= 1, irow
359 do kk= 1, nndof
360 ww(kk,av)= ww(kk,av) + ww(jj,y)*ww(kk,v+jj-1)
361 enddo
362 enddo
363
364 call hecmw_precond_apply(hecmesh, hecmat, ww(:,av), ww(:,zq), ww(:,zp), tcomm)
365
366 do kk= 1, nndof
367 x(kk)= x(kk) + ww(kk,zq)
368 enddo
369
370 !C
371 !C-- Compute residual vector R, find norm, then check for tolerance.
372 call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
373
374 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
375
376 ww(i+1,s)= dsqrt(dnrm2/bnrm2)
377 resid = ww( i+1,s )
378
379 if ( resid.le.tol ) exit outer
380 if ( iter .gt.maxit ) then
382 exit outer
383 end if
384 !C
385 !C-- RESTART
386 end do outer
387
388 !C
389 !C-- iteration FAILED
390
391 if (error == hecmw_solver_error_noconv_maxit) then
392 info = iter
393
394 !C-- [H]{y}= {s_tld}
395 do ik= 1, i
396 ss(ik)= ww(ik,s)
397 enddo
398 irow= i
399 ww(irow,y)= ss(irow) / h(irow,irow)
400
401 do kk= irow-1, 1, -1
402 do jj= irow, kk+1, -1
403 ss(kk)= ss(kk) - h(kk,jj)*ww(jj,y)
404 enddo
405 ww(kk,y)= ss(kk) / h(kk,kk)
406 enddo
407
408 !C-- {x}= {x} + {y}{V}
409 do kk= 1, nndof
410 ww(kk, av)= 0.d0
411 enddo
412
413 jj= irow
414 do jj= 1, irow
415 do kk= 1, nndof
416 ww(kk,av)= ww(kk ,av) + ww(jj,y)*ww(kk ,v+jj-1)
417 enddo
418 enddo
419
420 call hecmw_precond_apply(hecmesh, hecmat, ww(:,av), ww(:,zq), ww(:,zp), tcomm)
421
422 do kk= 1, nndof
423 x(kk)= x(kk) + ww(kk,zq)
424 enddo
425 end if
426
427 call hecmw_solver_scaling_bk(hecmat)
428
429 if (estcond /= 0 .and. hecmesh%my_rank == 0) then
431 endif
432 !C
433 !C-- INTERFACE data EXCHANGE
434 s_time = hecmw_wtime()
435 call hecmw_update_m_r (hecmesh, x, hecmat%NP, hecmat%NDOF)
436 e_time = hecmw_wtime()
437 tcomm = tcomm + e_time - s_time
438
439 deallocate (h, ww, ss)
440 !call hecmw_precond_clear(hecMAT)
441
442 e1_time= hecmw_wtime()
443 if (timelog.eq.2) then
444 call hecmw_time_statistics(hecmesh, e1_time - s1_time, &
445 t_max, t_min, t_avg, t_sd)
446 if (hecmesh%my_rank.eq.0) then
447 write(*,*) 'Time solver iterations'
448 write(*,*) ' Max :',t_max
449 write(*,*) ' Min :',t_min
450 write(*,*) ' Avg :',t_avg
451 write(*,*) ' Std Dev :',t_sd
452 endif
453 tsol = t_max
454 else
455 tsol = e1_time - s1_time
456 endif
457
458 end subroutine hecmw_solve_gmres
459
460end module hecmw_solver_gmres
subroutine hecmw_estimate_condition_gmres(i, h)
integer(kind=kint) function, public hecmw_mat_get_nrest(hecmat)
integer(kind=kint) function, public hecmw_mat_get_timelog(hecmat)
integer(kind=kint) function, public hecmw_mat_get_estcond(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, public hecmw_solve_gmres(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_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=4), parameter kreal
real(kind=kreal) function hecmw_wtime()
subroutine hecmw_barrier(hecmesh)
subroutine hecmw_update_m_r(hecmesh, val, n, m)
integer(kind=kint), parameter hecmw_solver_error_noconv_maxit