FrontISTR 5.2.0
Large-scale structural analysis program with finit element method
Loading...
Searching...
No Matches
hecmw_solver_Iterative.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
7contains
8 !
9 !C***
10 !C*** hecmw_solve_nn
11 !C***
12 !
13 subroutine hecmw_solve_iterative (hecMESH, hecMAT)
14
15 use hecmw_util
27
28 implicit none
29
30 type (hecmwST_matrix), target :: hecMAT
31 type (hecmwST_local_mesh) :: hecMESH
32
33 integer(kind=kint) :: error
34 integer(kind=kint) :: ITER, METHOD, PRECOND, NSET, METHOD2
35 integer(kind=kint) :: iterPREmax
36 integer(kind=kint) :: ITERlog, TIMElog
37 real(kind=kreal) :: resid, sigma_diag, thresh, filter, resid2
38 real(kind=kreal) :: time_setup, time_comm, time_sol, tr
39 real(kind=kreal) :: time_ax, time_precond
40
41 integer(kind=kint) :: NREST
42 real(kind=kreal) :: sigma
43
44 integer(kind=kint) :: totalmpc, MPC_METHOD
45 integer(kind=kint) :: auto_sigma_diag
46
47 !C PARAMETERs
48 iter = hecmw_mat_get_iter(hecmat)
49 method = hecmw_mat_get_method(hecmat)
50 method2 = hecmw_mat_get_method2(hecmat)
51 precond = hecmw_mat_get_precond(hecmat)
52 nset = hecmw_mat_get_nset(hecmat)
53 iterpremax= hecmw_mat_get_iterpremax(hecmat)
54 nrest = hecmw_mat_get_nrest(hecmat)
55 iterlog = hecmw_mat_get_iterlog(hecmat)
56 timelog = hecmw_mat_get_timelog(hecmat)
57 time_setup= 0.d0
58 time_comm = 0.d0
59 time_sol = 0.d0
60 resid = hecmw_mat_get_resid(hecmat)
61 sigma_diag= hecmw_mat_get_sigma_diag(hecmat)
62 sigma = hecmw_mat_get_sigma(hecmat)
63 thresh = hecmw_mat_get_thresh(hecmat)
64 filter = hecmw_mat_get_filter(hecmat)
65 if (sigma_diag.lt.0.d0) then
66 auto_sigma_diag= 1
67 sigma_diag= 1.d0
68 else
69 auto_sigma_diag= 0
70 endif
71
72 !C ERROR CHECK
73 call hecmw_solve_check_zerodiag(hecmesh, hecmat) !C-- ZERO DIAGONAL component
74
75 !C-- IN CASE OF MPC-CG
76 totalmpc = hecmesh%mpc%n_mpc
77 call hecmw_allreduce_i1 (hecmesh, totalmpc, hecmw_sum)
78 mpc_method = hecmw_mat_get_mpc_method(hecmat)
79 if (totalmpc > 0 .and. mpc_method == 2) then
80 call hecmw_mat_set_flag_mpcmatvec(hecmat, 1)
81 endif
82
83 !C-- RECYCLE SETTING OF PRECONDITIONER
85
86 ! exchange diagonal elements of overlap region
87 call hecmw_mat_dump(hecmat, hecmesh)
88 call hecmw_matvec_set_async(hecmat)
89
90 !C ITERATIVE solver
91 error=0
92 !! Auto Sigma_diag loop
93 do
94 call hecmw_mat_set_flag_converged(hecmat, 0)
95 call hecmw_mat_set_flag_diverged(hecmat, 0)
96 if (auto_sigma_diag.eq.1) call hecmw_mat_set_sigma_diag(hecmat, sigma_diag)
97
100 call hecmw_solve_iterative_printmsg(hecmesh,hecmat,method)
101
102 select case(method)
103 case (1) !--CG
104 hecmat%symmetric = .true.
105 call hecmw_solve_cg( hecmesh, hecmat, iter, resid, error, time_setup, time_sol, time_comm )
106 case (2) !--BiCGSTAB
107 hecmat%symmetric = .false.
108 call hecmw_solve_bicgstab( hecmesh,hecmat, iter, resid, error,time_setup, time_sol, time_comm )
109 case (3) !--GMRES
110 hecmat%symmetric = .false.
111 call hecmw_solve_gmres( hecmesh,hecmat, iter, resid, error, time_setup, time_sol, time_comm )
112 case (4) !--GPBiCG
113 hecmat%symmetric = .false.
114 call hecmw_solve_gpbicg( hecmesh,hecmat, iter, resid, error, time_setup, time_sol, time_comm )
115 case default
116 error = hecmw_solver_error_incons_pc !!未定義なMETHOD!!
117 call hecmw_solve_error (hecmesh, error)
118 end select
119
121 .or. error==hecmw_solver_error_diverge_nan) then
122 call hecmw_mat_set_flag_diverged(hecmat, 1)
123 if ((precond>=10 .and. precond<20) .and. auto_sigma_diag==1 .and. sigma_diag<2.d0) then
124 sigma_diag = sigma_diag + 0.1
125 if (hecmesh%my_rank.eq.0) write(*,*) 'Increasing SIGMA_DIAG to', sigma_diag
126 cycle
127 elseif (method==1 .and. method2>1) then
128 if (auto_sigma_diag.eq.1) sigma_diag = 1.0
129 method = method2
130 cycle
131 endif
132 endif
133
134 if (auto_sigma_diag.eq.1) call hecmw_mat_set_sigma_diag(hecmat, -1.d0)
135 exit
136 enddo
137
138 if (error.ne.0) then
139 call hecmw_solve_error (hecmesh, error)
140 endif
141
142 resid2=hecmw_rel_resid_l2(hecmesh,hecmat)
143 if (hecmesh%my_rank.eq.0 .and. (iterlog.eq.1 .or. timelog.ge.1)) then
144 write(*,"(a,1pe12.5)")'### Relative residual =', resid2
145 endif
146 if (resid2 < hecmw_mat_get_resid(hecmat)) call hecmw_mat_set_flag_converged(hecmat, 1)
147
148 call hecmw_mat_dump_solution(hecmat)
150
151 !C-- IN CASE OF MPC-CG
152 if (totalmpc > 0 .and. mpc_method == 2) then
153 call hecmw_mat_set_flag_mpcmatvec(hecmat, 0)
154 endif
155
156 time_ax = hecmw_matvec_get_timer()
157 time_precond = hecmw_precond_get_timer()
158
159 if (hecmesh%my_rank.eq.0 .and. timelog.ge.1) then
160 tr= (time_sol-time_comm)/(time_sol+1.d-24)*100.d0
161 write (*,'(/a)') '### summary of linear solver'
162 write (*,'(i10,a, 1pe16.6)') iter, ' iterations ', resid
163 write (*,'(a, 1pe16.6 )') ' set-up time : ', time_setup
164 write (*,'(a, 1pe16.6 )') ' solver time : ', time_sol
165 write (*,'(a, 1pe16.6 )') ' solver/comm time : ', time_comm
166 write (*,'(a, 1pe16.6 )') ' solver/matvec : ', time_ax
167 write (*,'(a, 1pe16.6 )') ' solver/precond : ', time_precond
168 if (iter > 0) &
169 write (*,'(a, 1pe16.6 )') ' solver/1 iter : ', time_sol / iter
170 write (*,'(a, 1pe16.6/)') ' work ratio (%) : ', tr
171 endif
172
173 end subroutine hecmw_solve_iterative
174
175 subroutine hecmw_solve_check_zerodiag (hecMESH, hecMAT)
176 use hecmw_util
180 implicit none
181 type (hecmwST_local_mesh) :: hecMESH
182 type (hecmwST_matrix), target :: hecMAT
183 integer (kind=kint)::PRECOND,iterPREmax,i,j,error
184 precond = hecmw_mat_get_precond(hecmat)
185 iterpremax= hecmw_mat_get_iterpremax(hecmat)
186 !C
187 !C-- ZERO DIAGONAL component
188 error= 0
189 do i= 1, hecmat%N
190 do j = 1, hecmat%NDOF
191 if (dabs(hecmat%D(hecmat%NDOF*hecmat%NDOF*(i-1)+(j-1)*(hecmat%NDOF+1)+1)).eq.0.d0) then
193 end if
194 end do
195 enddo
196
197 call hecmw_allreduce_i1 (hecmesh, error, hecmw_max)
198 if (error.ne.0 .and. (precond.lt.10 .and. iterpremax.gt.0)) then
199 call hecmw_solve_error (hecmesh, error)
200 endif
201
202 end subroutine hecmw_solve_check_zerodiag
203
204 function hecmw_solve_check_zerorhs (hecMESH, hecMAT)
205 use hecmw_util
209 implicit none
210 type (hecmwst_local_mesh) :: hecmesh
211 type (hecmwst_matrix), target :: hecmat
212 real(kind=kreal), dimension(1) :: rhs
213 integer (kind=kint)::precond,iterpremax,i,j,error
215
216 precond = hecmw_mat_get_precond(hecmat)
217 iterpremax= hecmw_mat_get_iterpremax(hecmat)
218 !C
219 !C-- ZERO RHS norm
220 error= 0
222
223 rhs(1)= 0.d0
224 do i= 1, hecmat%N
225 do j = 1, hecmat%NDOF
226 rhs(1)=rhs(1) + hecmat%B(hecmat%NDOF*(i-1)+j)**2
227 end do
228 enddo
229 if (hecmesh%mpc%n_mpc > 0) then
230 do i= 1, hecmesh%mpc%n_mpc
231 rhs(1)= rhs(1) + hecmesh%mpc%mpc_const(i)**2
232 enddo
233 endif
234 call hecmw_allreduce_r (hecmesh, rhs, 1, hecmw_sum)
235
236 if (rhs(1).eq.0.d0) then
238 call hecmw_solve_error (hecmesh, error)
239 hecmat%X(:)=0.d0
241 endif
242
243 end function hecmw_solve_check_zerorhs
244
245 subroutine hecmw_solve_iterative_printmsg (hecMESH, hecMAT, METHOD)
246 use hecmw_util
249
250 implicit none
251 type (hecmwST_local_mesh) :: hecMESH
252 type (hecmwST_matrix), target :: hecMAT
253 integer(kind=kint) :: METHOD
254 integer(kind=kint) :: ITER, PRECOND, NSET, iterPREmax, NREST
255 integer(kind=kint) :: ITERlog, TIMElog
256
257 character(len=30) :: msg_precond
258 character(len=30) :: msg_method
259
260 iter = hecmw_mat_get_iter(hecmat)
261 ! METHOD = hecmw_mat_get_method(hecMAT)
262 precond = hecmw_mat_get_precond(hecmat)
263 nset = hecmw_mat_get_nset(hecmat)
264 iterpremax= hecmw_mat_get_iterpremax(hecmat)
265 nrest = hecmw_mat_get_nrest(hecmat)
266 iterlog= hecmw_mat_get_iterlog(hecmat)
267 timelog= hecmw_mat_get_timelog(hecmat)
268
269 select case(method)
270 case (1) !--CG
271 msg_method="CG"
272 case (2) !--BiCGSTAB
273 msg_method="BiCGSTAB"
274 case (3) !--GMRES
275 msg_method="GMRES"
276 case (4) !--GPBiCG
277 msg_method="GPBiCG"
278 case default
279 msg_method="Unlabeled"
280 end select
281 select case(precond)
282 case (1,2)
283 msg_precond="SSOR"
284 case (3)
285 msg_precond="DIAG"
286 case (5)
287 msg_precond="ML"
288 case (7)
289 msg_precond="DirectMUMPS"
290 case (10, 11, 12)
291 write(msg_precond,"(a,i0,a)") "ILU(",precond-10,")"
292 case (20)
293 msg_precond="SAINV"
294 case (21)
295 msg_precond="RIF"
296 case default
297 msg_precond="Unlabeled"
298 end select
299 if (hecmesh%my_rank.eq.0 .and. (iterlog.eq.1 .or. timelog.ge.1)) then
300 write (*,'(a,i0,a,i0,a,a,a,a,a,i0)') '### ',hecmat%NDOF,'x',hecmat%NDOF,' BLOCK ', &
301 & trim(msg_method),", ",trim(msg_precond),", ", iterpremax
302 end if
303 end subroutine hecmw_solve_iterative_printmsg
304
305end module hecmw_solver_iterative
subroutine, public hecmw_mat_dump(hecmat, hecmesh)
subroutine, public hecmw_mat_dump_solution(hecmat)
integer(kind=kint) function, public hecmw_mat_get_iterpremax(hecmat)
integer(kind=kint) function, public hecmw_mat_get_nrest(hecmat)
integer(kind=kint) function, public hecmw_mat_get_method2(hecmat)
integer(kind=kint) function, public hecmw_mat_get_method(hecmat)
subroutine, public hecmw_mat_set_sigma_diag(hecmat, sigma_diag)
subroutine, public hecmw_mat_recycle_precond_setting(hecmat)
integer(kind=kint) function, public hecmw_mat_get_timelog(hecmat)
real(kind=kreal) function, public hecmw_mat_get_resid(hecmat)
subroutine, public hecmw_mat_set_flag_converged(hecmat, flag_converged)
integer(kind=kint) function, public hecmw_mat_get_iter(hecmat)
subroutine, public hecmw_mat_set_flag_mpcmatvec(hecmat, flag_mpcmatvec)
real(kind=kreal) function, public hecmw_mat_get_sigma(hecmat)
integer(kind=kint) function, public hecmw_mat_get_mpc_method(hecmat)
real(kind=kreal) function, public hecmw_mat_get_thresh(hecmat)
subroutine, public hecmw_mat_set_flag_diverged(hecmat, flag_diverged)
integer(kind=kint) function, public hecmw_mat_get_precond(hecmat)
integer(kind=kint) function, public hecmw_mat_get_nset(hecmat)
real(kind=kreal) function, public hecmw_mat_get_filter(hecmat)
integer(kind=kint) function, public hecmw_mat_get_iterlog(hecmat)
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecmat)
subroutine, public hecmw_precond_clear_timer
real(kind=kreal) function, public hecmw_precond_get_timer()
subroutine hecmw_solve_bicgstab(hecmesh, hecmat, iter, resid, error, tset, tsol, tcomm)
subroutine, public hecmw_solve_cg(hecmesh, hecmat, iter, resid, error, tset, tsol, tcomm)
subroutine, public hecmw_solve_gmres(hecmesh, hecmat, iter, resid, error, tset, tsol, tcomm)
subroutine, public hecmw_solve_gpbicg(hecmesh, hecmat, iter, resid, error, tset, tsol, tcomm)
subroutine hecmw_solve_check_zerodiag(hecmesh, hecmat)
logical function hecmw_solve_check_zerorhs(hecmesh, hecmat)
subroutine hecmw_solve_iterative_printmsg(hecmesh, hecmat, method)
subroutine hecmw_solve_iterative(hecmesh, hecmat)
subroutine, public hecmw_matvec_clear_timer
real(kind=kreal) function, public hecmw_rel_resid_l2(hecmesh, hecmat, commtime)
subroutine, public hecmw_matvec_set_async(hecmat)
real(kind=kreal) function, public hecmw_matvec_get_timer()
subroutine, public hecmw_matvec_unset_async
I/O and Utility.
Definition: hecmw_util_f.F90:7
integer(kind=kint), parameter hecmw_sum
integer(kind=kint), parameter hecmw_max
integer(kind=4), parameter kreal
subroutine hecmw_allreduce_r(hecmesh, val, n, ntag)
subroutine hecmw_allreduce_i1(hecmesh, s, ntag)
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_incons_pc
integer(kind=kint), parameter hecmw_solver_error_zero_diag
subroutine hecmw_solve_error(hecmesh, iflag)
integer(kind=kint), parameter hecmw_solver_error_zero_rhs
integer(kind=kint), parameter hecmw_solver_error_diverge_mat