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