1 !-------------------------------------------------------------------------------
2 ! Copyright (c) 2016 The University of Tokyo
3 ! This software is released under the MIT License, see LICENSE.txt
4 !-------------------------------------------------------------------------------
5
6 !C
7 !C***
8 !C*** module hecmw_solver_CG_33
9 !C***
10 !C
11 module hecmw_solver_CG_33
12
13 public :: hecmw_solve_CG_33
14
15 contains
16 !C
17 !C*** CG_33
18 !C
19 subroutine hecmw_solve_CG_33( hecMESH, hecMAT, ITER, RESID, ERROR, &
20 & Tset, Tsol, Tcomm )
21
22 use hecmw_util
23 use m_hecmw_solve_error
24 use m_hecmw_comm_f
25 use hecmw_matrix_misc
26 use hecmw_solver_misc
27 use hecmw_solver_las_33
28 use hecmw_solver_scaling_33
29 use hecmw_precond_33
30 use hecmw_jad_type
31 use hecmw_estimate_condition
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
94 allocate (WW(NDOF*NP, 4))
95 WW = 0.d0
96
97 !C
98 !C-- SCALING
99 call hecmw_solver_scaling_fw_33(hecMESH, hecMAT, Tcomm)
100
101 IF (hecmw_mat_get_usejad(hecMAT).ne.0) THEN
102 call hecmw_JAD_INIT(hecMAT)
103 ENDIF
104
105 if (ESTCOND /= 0 .and. hecMESH%my_rank == 0) then
106 allocate(D(MAXIT),E(MAXIT-1))
107 endif
108 !C===
109 !C +----------------------+
110 !C | SETUP PRECONDITIONER |
111 !C +----------------------+
112 !C===
113 call hecmw_precond_33_setup(hecMAT, hecMESH, 1)
114
115 !C===
116 !C +---------------------+
117 !C | {r0}= {b} - [A]{x0} |
118 !C +---------------------+
119 !C===
120 call hecmw_matresid_33(hecMESH, hecMAT, X, B, WW(:,R), Tcomm)
121
122 !C-- compute ||{b}||
123 call hecmw_InnerProduct_R(hecMESH, NDOF, B, B, BNRM2, Tcomm)
124 if (BNRM2.eq.0.d0) then
125 iter = 0
126 MAXIT = 0
127 RESID = 0.d0
128 X = 0.d0
129 endif
130
131 E_TIME = HECMW_WTIME()
132 if (TIMElog.eq.2) then
133 call hecmw_time_statistics(hecMESH, E_TIME - S_TIME, &
134 t_max, t_min, t_avg, t_sd)
135 if (hecMESH%my_rank.eq.0) then
136 write(*,*) 'Time solver setup'
137 write(*,*) ' Max :',t_max
138 write(*,*) ' Min :',t_min
139 write(*,*) ' Avg :',t_avg
140 write(*,*) ' Std Dev :',t_sd
141 endif
142 Tset = t_max
143 else
144 Tset = E_TIME - S_TIME
145 endif
146
147 Tcomm = 0.d0
148 call hecmw_barrier(hecMESH)
149 S1_TIME = HECMW_WTIME()
150 !C
151 !C************************************************* Conjugate Gradient Iteration start
152 !C
153 do iter = 1, MAXIT
154
155 !C===
156 !C +----------------+
157 !C | {z}= [Minv]{r} |
158 !C +----------------+
159 !C===
160 call hecmw_precond_33_apply(hecMESH, hecMAT, WW(:,R), WW(:,Z), WW(:,WK), Tcomm)
161
162 !C===
163 !C +---------------+
164 !C | {RHO}= {r}{z} |
165 !C +---------------+
166 !C===
167 call hecmw_InnerProduct_R(hecMESH, NDOF, WW(:,R), WW(:,Z), RHO, Tcomm)
168 ! if RHO is NaN or Inf then no converge
169 if (RHO == 0.d0) then
170 ! converged due to RHO==0
171 exit
172 elseif (iter > 1 .and. RHO*RHO1 <= 0) then
173 n_indef_precond = n_indef_precond + 1
174 if (n_indef_precond >= 3) then
175 ! diverged due to indefinite preconditioner
176 ERROR = HECMW_SOLVER_ERROR_DIVERGE_PC
177 exit
178 endif
179 endif
180
181 !C===
182 !C +-----------------------------+
183 !C | {p} = {z} if ITER=1 |
184 !C | BETA= RHO / RHO1 otherwise |
185 !C +-----------------------------+
186 !C===
187 if ( ITER.eq.1 ) then
188 do i = 1, NNDOF
189 WW(i,P) = WW(i,Z)
190 enddo
191 else
192 BETA = RHO / RHO1
193 do i = 1, NNDOF
194 WW(i,P) = WW(i,Z) + BETA*WW(i,P)
195 enddo
196 endif
197
198 !C===
199 !C +--------------+
200 !C | {q}= [A] {p} |
201 !C +--------------+
202 !C===
203 call hecmw_matvec_33(hecMESH, hecMAT, WW(:,P), WW(:,Q), Tcomm)
204
205 !C===
206 !C +---------------------+
207 !C | ALPHA= RHO / {p}{q} |
208 !C +---------------------+
209 !C===
210 call hecmw_InnerProduct_R(hecMESH, NDOF, WW(:,P), WW(:,Q), C1, Tcomm)
211 ! if C1 is NaN or Inf, no converge
212 if (C1 <= 0) then
213 ! diverged due to indefinite or negative definite matrix
214 ERROR = HECMW_SOLVER_ERROR_DIVERGE_MAT
215 exit
216 endif
217
218 ALPHA= RHO / C1
219
220 !C===
221 !C +----------------------+
222 !C | {x}= {x} + ALPHA*{p} |
223 !C | {r}= {r} - ALPHA*{q} |
224 !C +----------------------+
225 !C===
226 do i = 1, NNDOF
227 X(i) = X(i) + ALPHA * WW(i,P)
228 ! WW(i,R)= WW(i,R) - ALPHA * WW(i,Q)
229 enddo
230
231 if ( mod(ITER,N_ITER_RECOMPUTE_R)==0 ) then
232 call hecmw_matresid_33(hecMESH, hecMAT, X, B, WW(:,R), Tcomm)
233 else
234 do i = 1, NNDOF
235 WW(i,R)= WW(i,R) - ALPHA * WW(i,Q)
236 enddo
237 endif
238
239 call hecmw_InnerProduct_R(hecMESH, NDOF, WW(:,R), WW(:,R), DNRM2, Tcomm)
240
241 RESID= dsqrt(DNRM2/BNRM2)
242
243 !C##### ITERATION HISTORY
244 if (my_rank.eq.0.and.ITERLog.eq.1) write (*,'(i7, 1pe16.6)') ITER, RESID
245 !C#####
246
247 if (ESTCOND /= 0 .and. hecMESH%my_rank == 0) then
248 if (ITER == 1) then
249 D(1) = 1.d0 / ALPHA
250 else
251 D(ITER) = 1.d0 / ALPHA + BETA / ALPHA1
252 E(ITER-1) = sqrt(BETA) / ALPHA1
253 endif
254 ALPHA1 = ALPHA
255 if (mod(ITER,ESTCOND) == 0) call hecmw_estimate_condition_CG(ITER, D, E)
256 endif
257
258 if ( RESID.le.TOL ) then
259 if ( mod(ITER,N_ITER_RECOMPUTE_R)==0 ) exit
260 !C----- recompute R to make sure it is really converged
261 call hecmw_matresid_33(hecMESH, hecMAT, X, B, WW(:,R), Tcomm)
262 call hecmw_InnerProduct_R(hecMESH, NDOF, WW(:,R), WW(:,R), DNRM2, Tcomm)
263 RESID= dsqrt(DNRM2/BNRM2)
264 if ( RESID.le.TOL ) exit
265 endif
266 if ( ITER .eq.MAXIT ) ERROR = HECMW_SOLVER_ERROR_NOCONV_MAXIT
267
268 RHO1 = RHO
269
270 enddo
271 !C
272 !C************************************************* Conjugate Gradient Iteration end
273 !C
274 call hecmw_solver_scaling_bk_33(hecMAT)
275 !C
276 !C-- INTERFACE data EXCHANGE
277 !C
278 START_TIME= HECMW_WTIME()
279 call hecmw_update_3_R (hecMESH, X, hecMAT%NP)
280 END_TIME = HECMW_WTIME()
281 Tcomm = Tcomm + END_TIME - START_TIME
282
283 deallocate (WW)
284 !call hecmw_precond_33_clear(hecMAT)
285
286 IF (hecmw_mat_get_usejad(hecMAT).ne.0) THEN
287 call hecmw_JAD_FINALIZE()
288 ENDIF
289
290 if (ESTCOND /= 0 .and. ERROR == 0 .and. hecMESH%my_rank == 0) then
291 call hecmw_estimate_condition_CG(ITER, D, E)
292 deallocate(D, E)
293 endif
294
295 E1_TIME = HECMW_WTIME()
296 if (TIMElog.eq.2) then
297 call hecmw_time_statistics(hecMESH, E1_TIME - S1_TIME, &
298 t_max, t_min, t_avg, t_sd)
299 if (hecMESH%my_rank.eq.0) then
300 write(*,*) 'Time solver iterations'
301 write(*,*) ' Max :',t_max
302 write(*,*) ' Min :',t_min
303 write(*,*) ' Avg :',t_avg
304 write(*,*) ' Std Dev :',t_sd
305 endif
306 Tsol = t_max
307 else
308 Tsol = E1_TIME - S1_TIME
309 endif
310
311 end subroutine hecmw_solve_CG_33
312
313 end module hecmw_solver_CG_33