FrontISTR  5.2.0
Large-scale structural analysis program with finit element method
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 
7 contains
8  !
9  !C***
10  !C*** hecmw_solve_nn
11  !C***
12  !
13  subroutine hecmw_solve_iterative (hecMESH, hecMAT)
14 
15  use hecmw_util
16  use hecmw_solver_cg
21  use m_hecmw_comm_f
23  use hecmw_precond
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
179  use m_hecmw_comm_f
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
208  use m_hecmw_comm_f
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
214  logical :: hecmw_solve_check_zerorhs
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
221  hecmw_solve_check_zerorhs = .false.
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 
305 end module hecmw_solver_iterative
subroutine, public hecmw_mat_dump(hecMAT, hecMESH)
subroutine, public hecmw_mat_dump_solution(hecMAT)
real(kind=kreal) function, public hecmw_mat_get_sigma_diag(hecMAT)
integer(kind=kint) function, public hecmw_mat_get_iterpremax(hecMAT)
integer(kind=kint) function, public hecmw_mat_get_nrest(hecMAT)
subroutine, public hecmw_mat_set_flag_diverged(hecMAT, flag_diverged)
real(kind=kreal) function, public hecmw_mat_get_resid(hecMAT)
subroutine, public hecmw_mat_set_flag_converged(hecMAT, flag_converged)
subroutine, public hecmw_mat_recycle_precond_setting(hecMAT)
subroutine, public hecmw_mat_set_sigma_diag(hecMAT, sigma_diag)
integer(kind=kint) function, public hecmw_mat_get_iterlog(hecMAT)
integer(kind=kint) function, public hecmw_mat_get_timelog(hecMAT)
integer(kind=kint) function, public hecmw_mat_get_method2(hecMAT)
real(kind=kreal) function, public hecmw_mat_get_thresh(hecMAT)
subroutine, public hecmw_mat_set_flag_mpcmatvec(hecMAT, flag_mpcmatvec)
integer(kind=kint) function, public hecmw_mat_get_method(hecMAT)
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_mpc_method(hecMAT)
real(kind=kreal) function, public hecmw_mat_get_sigma(hecMAT)
integer(kind=kint) function, public hecmw_mat_get_iter(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_iterative_printmsg(hecMESH, hecMAT, METHOD)
logical function hecmw_solve_check_zerorhs(hecMESH, hecMAT)
subroutine hecmw_solve_check_zerodiag(hecMESH, hecMAT)
subroutine hecmw_solve_iterative(hecMESH, hecMAT)
subroutine, public hecmw_matvec_clear_timer
subroutine, public hecmw_matvec_set_async(hecMAT)
real(kind=kreal) function, public hecmw_rel_resid_l2(hecMESH, hecMAT, COMMtime)
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_i1(hecMESH, s, ntag)
subroutine hecmw_allreduce_r(hecMESH, val, n, 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
integer(kind=kint), parameter hecmw_solver_error_zero_rhs
subroutine hecmw_solve_error(hecMESH, IFLAG)
integer(kind=kint), parameter hecmw_solver_error_diverge_mat