FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_solver_CG.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_CG
9!C***
10!C
12
13 public :: hecmw_solve_cg
14
15contains
16 !C
17 !C*** CG_nn
18 !C
19 subroutine hecmw_solve_cg( hecMESH, hecMAT, ITER, RESID, error, &
20 & Tset, Tsol, Tcomm )
21
22 use hecmw_util
32
33 implicit none
34
35 type(hecmwst_local_mesh) :: hecmesh
36 type(hecmwst_matrix) :: hecmat
37 integer(kind=kint ), intent(inout):: iter, error
38 real (kind=kreal), intent(inout):: resid, tset, tsol, tcomm
39
40 integer(kind=kint ) :: n, np, ndof, nndof
41 integer(kind=kint ) :: my_rank
42 integer(kind=kint ) :: iterlog, timelog
43 real(kind=kreal), pointer :: b(:), x(:)
44
45 real(kind=kreal), dimension(:,:), allocatable :: ww
46
47 integer(kind=kint), parameter :: r= 1
48 integer(kind=kint), parameter :: z= 2
49 integer(kind=kint), parameter :: q= 2
50 integer(kind=kint), parameter :: p= 3
51 integer(kind=kint), parameter :: wk= 4
52
53 integer(kind=kint ) :: maxit
54
55 ! local variables
56 real (kind=kreal) :: tol
57 integer(kind=kint )::i
58 real (kind=kreal)::s_time,s1_time,e_time,e1_time, start_time, end_time
59 real (kind=kreal)::bnrm2
60 real (kind=kreal)::rho,rho1,beta,c1,alpha,dnrm2
61 real (kind=kreal)::t_max,t_min,t_avg,t_sd
62 integer(kind=kint) ::estcond
63 real (kind=kreal), allocatable::d(:),e(:)
64 real (kind=kreal)::alpha1
65 integer(kind=kint) :: n_indef_precond
66
67 integer(kind=kint), parameter :: n_iter_recompute_r= 50
68
69 call hecmw_barrier(hecmesh)
70 s_time= hecmw_wtime()
71
72 !C===
73 !C +-------+
74 !C | INIT. |
75 !C +-------+
76 !C===
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 estcond = hecmw_mat_get_estcond( hecmat )
90
91 error = 0
92 n_indef_precond = 0
93 rho1 = 0.0d0
94 alpha1 = 0.0d0
95 beta = 0.0d0
96
97 allocate (ww(ndof*np, 4))
98 ww = 0.d0
99
100 !C
101 !C-- SCALING
102 call hecmw_solver_scaling_fw(hecmesh, hecmat, tcomm)
103
104 if (hecmw_mat_get_usejad(hecmat).ne.0) then
105 call hecmw_jad_init(hecmat)
106 endif
107
108 if (estcond /= 0 .and. hecmesh%my_rank == 0) then
109 allocate(d(maxit),e(maxit-1))
110 endif
111 !C===
112 !C +----------------------+
113 !C | SETUP PRECONDITIONER |
114 !C +----------------------+
115 !C===
116 call hecmw_precond_setup(hecmat, hecmesh, 1)
117
118 !C===
119 !C +---------------------+
120 !C | {r0}= {b} - [A]{x0} |
121 !C +---------------------+
122 !C===
123 call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
124
125 !C-- compute ||{b}||
126 call hecmw_innerproduct_r(hecmesh, ndof, b, b, bnrm2, tcomm)
127 if (bnrm2.eq.0.d0) then
128 iter = 0
129 maxit = 0
130 resid = 0.d0
131 x = 0.d0
132 endif
133
134 e_time = hecmw_wtime()
135 if (timelog.eq.2) then
136 call hecmw_time_statistics(hecmesh, e_time - s_time, &
137 t_max, t_min, t_avg, t_sd)
138 if (hecmesh%my_rank.eq.0) then
139 write(*,*) 'Time solver setup'
140 write(*,*) ' Max :',t_max
141 write(*,*) ' Min :',t_min
142 write(*,*) ' Avg :',t_avg
143 write(*,*) ' Std Dev :',t_sd
144 endif
145 tset = t_max
146 else
147 tset = e_time - s_time
148 endif
149
150 tcomm = 0.d0
151 call hecmw_barrier(hecmesh)
152 s1_time = hecmw_wtime()
153 !C
154 !C************************************************* Conjugate Gradient Iteration start
155 !C
156 do iter = 1, maxit
157
158 !C===
159 !C +----------------+
160 !C | {z}= [Minv]{r} |
161 !C +----------------+
162 !C===
163 call hecmw_precond_apply(hecmesh, hecmat, ww(:,r), ww(:,z), ww(:,wk), tcomm)
164
165
166 !C===
167 !C +---------------+
168 !C | {RHO}= {r}{z} |
169 !C +---------------+
170 !C===
171 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,z), rho, tcomm)
172 ! if RHO is NaN or Inf then no converge
173 if (rho == 0.d0) then
174 ! converged due to RHO==0
175 exit
176 elseif (iter > 1 .and. rho*rho1 <= 0) then
177 n_indef_precond = n_indef_precond + 1
178 if (n_indef_precond >= 3) then
179 ! diverged due to indefinite preconditioner
181 exit
182 endif
183 elseif (rho /= rho) then ! RHO is NaN
185 exit
186 endif
187
188 !C===
189 !C +-----------------------------+
190 !C | {p} = {z} if ITER=1 |
191 !C | BETA= RHO / RHO1 otherwise |
192 !C +-----------------------------+
193 !C===
194 if ( iter.eq.1 ) then
195 call hecmw_copy_r(hecmesh, ndof, ww(:,z), ww(:,p))
196 else
197 beta = rho / rho1
198 call hecmw_xpay_r(hecmesh, ndof, beta, ww(:,z), ww(:,p))
199 endif
200
201 !C===
202 !C +--------------+
203 !C | {q}= [A] {p} |
204 !C +--------------+
205 !C===
206 call hecmw_matvec(hecmesh, hecmat, ww(:,p), ww(:,q), tcomm)
207
208 !C===
209 !C +---------------------+
210 !C | ALPHA= RHO / {p}{q} |
211 !C +---------------------+
212 !C===
213 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,p), ww(:,q), c1, tcomm)
214 ! if C1 is NaN or Inf, no converge
215 if (c1 <= 0) then
216 ! diverged due to indefinite or negative definite matrix
218 exit
219 elseif (c1 /= c1) then ! C1 is NaN
221 exit
222 endif
223
224 alpha= rho / c1
225
226 !C===
227 !C +----------------------+
228 !C | {x}= {x} + ALPHA*{p} |
229 !C | {r}= {r} - ALPHA*{q} |
230 !C +----------------------+
231 !C===
232 call hecmw_axpy_r(hecmesh, ndof, alpha, ww(:,p), x)
233
234 if ( mod(iter,n_iter_recompute_r)==0 ) then
235 call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
236 else
237 call hecmw_axpy_r(hecmesh, ndof, -alpha, ww(:,q), ww(:,r))
238 endif
239
240 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
241
242 resid= dsqrt(dnrm2/bnrm2)
243
244 !C##### ITERATION HISTORY
245 if (my_rank.eq.0.and.iterlog.eq.1) write (*,'(i7, 1pe16.6)') iter, resid
246 !C#####
247
248 if (estcond /= 0 .and. hecmesh%my_rank == 0) then
249 if (iter == 1) then
250 d(1) = 1.d0 / alpha
251 else
252 d(iter) = 1.d0 / alpha + beta / alpha1
253 e(iter-1) = sqrt(beta) / alpha1
254 endif
255 alpha1 = alpha
256 if (mod(iter,estcond) == 0) call hecmw_estimate_condition_cg(iter, d, e)
257 endif
258
259 if ( resid.le.tol ) then
260 if ( mod(iter,n_iter_recompute_r)==0 ) exit
261 !C----- recompute R to make sure it is really converged
262 call hecmw_matresid(hecmesh, hecmat, x, b, ww(:,r), tcomm)
263 call hecmw_innerproduct_r(hecmesh, ndof, ww(:,r), ww(:,r), dnrm2, tcomm)
264 resid= dsqrt(dnrm2/bnrm2)
265 if ( resid.le.tol ) exit
266 endif
267 if ( iter .eq.maxit ) error = hecmw_solver_error_noconv_maxit
268
269 rho1 = rho
270
271 enddo
272 !C
273 !C************************************************* Conjugate Gradient Iteration end
274 !C
275 call hecmw_solver_scaling_bk(hecmat)
276 !C
277 !C-- INTERFACE data EXCHANGE
278 !C
279 start_time= hecmw_wtime()
280 call hecmw_update_m_r (hecmesh, x, hecmat%NP, hecmat%NDOF)
281 end_time = hecmw_wtime()
282 tcomm = tcomm + end_time - start_time
283
284 deallocate (ww)
285 !call hecmw_precond_clear(hecMAT)
286
287 if (hecmw_mat_get_usejad(hecmat).ne.0) then
288 call hecmw_jad_finalize(hecmat)
289 endif
290
291 if (estcond /= 0 .and. error == 0 .and. hecmesh%my_rank == 0) then
292 call hecmw_estimate_condition_cg(iter, d, e)
293 deallocate(d, e)
294 endif
295
296 e1_time = hecmw_wtime()
297 if (timelog.eq.2) then
298 call hecmw_time_statistics(hecmesh, e1_time - s1_time, &
299 t_max, t_min, t_avg, t_sd)
300 if (hecmesh%my_rank.eq.0) then
301 write(*,*) 'Time solver iterations'
302 write(*,*) ' Max :',t_max
303 write(*,*) ' Min :',t_min
304 write(*,*) ' Avg :',t_avg
305 write(*,*) ' Std Dev :',t_sd
306 endif
307 tsol = t_max
308 else
309 tsol = e1_time - s1_time
310 endif
311
312 end subroutine hecmw_solve_cg
313
314end module hecmw_solver_cg
subroutine hecmw_estimate_condition_cg(iter, d, e)
Jagged Diagonal Matrix storage for vector processors. Original code was provided by JAMSTEC.
Definition: hecmw_jadm.f90:8
subroutine, public hecmw_jad_init(hecmat)
Definition: hecmw_jadm.f90:29
subroutine, public hecmw_jad_finalize(hecmat)
Definition: hecmw_jadm.f90:42
integer(kind=kint) function, public hecmw_mat_get_usejad(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_cg(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_axpy_r(hecmesh, ndof, alpha, x, y)
subroutine hecmw_xpay_r(hecmesh, ndof, alpha, x, y)
subroutine hecmw_copy_r(hecmesh, ndof, x, y)
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_diverge_pc
integer(kind=kint), parameter hecmw_solver_error_diverge_nan
integer(kind=kint), parameter hecmw_solver_error_noconv_maxit
integer(kind=kint), parameter hecmw_solver_error_diverge_mat