1!-------------------------------------------------------------------------------
2! Copyright (c) 2019 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
6module hecmw_solver_iterative
7contains
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
17    use hecmw_solver_BiCGSTAB
18    use hecmw_solver_GMRES
19    use hecmw_solver_GPBiCG
20    use m_hecmw_solve_error
21    use m_hecmw_comm_f
22    use hecmw_solver_las
23    use hecmw_precond
24    use hecmw_matrix_misc
25    use hecmw_matrix_dump
26    use hecmw_solver_misc
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
84    call hecmw_mat_recycle_precond_setting(hecMAT)
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
98      call hecmw_matvec_clear_timer()
99      call hecmw_precond_clear_timer()
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
120      if (error==HECMW_SOLVER_ERROR_DIVERGE_PC .or. error==HECMW_SOLVER_ERROR_DIVERGE_MAT &
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)
149    call hecmw_matvec_unset_async
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
177    use hecmw_matrix_misc
178    use m_hecmw_solve_error
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
192          error=HECMW_SOLVER_ERROR_ZERO_DIAG
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
206    use hecmw_matrix_misc
207    use m_hecmw_solve_error
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
237      error= HECMW_SOLVER_ERROR_ZERO_RHS
238      call hecmw_solve_error (hecMESH, error)
239      hecMAT%X(:)=0.d0
240      hecmw_solve_check_zerorhs = .true.
241    endif
242
243  end function hecmw_solve_check_zerorhs
244
245  subroutine hecmw_solve_iterative_printmsg (hecMESH, hecMAT, METHOD)
246    use hecmw_util
247    use hecmw_solver_misc
248    use hecmw_matrix_misc
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
306