| 1 | !======================================================================!
|
| 2 | ! !
|
| 3 | ! Software Name : HEC-MW Library for PC-cluster !
|
| 4 | ! Version : 2.7 !
|
| 5 | ! !
|
| 6 | ! Last Update : 2014/12/14 !
|
| 7 | ! Category : Linear Solver !
|
| 8 | ! !
|
| 9 | ! Written by Kengo Nakajima (Univ. of Tokyo) !
|
| 10 | ! Kazuya Goto (PExProCS LLC) !
|
| 11 | ! !
|
| 12 | ! Contact address : IIS,The University of Tokyo RSS21 project !
|
| 13 | ! !
|
| 14 | ! "Structural Analysis System for General-purpose Coupling !
|
| 15 | ! Simulations Using High End Computing Middleware (HEC-MW)" !
|
| 16 | ! !
|
| 17 | !======================================================================!
|
| 18 |
|
| 19 | !C
|
| 20 | !C***
|
| 21 | !C*** module hecmw_solver_CG_33
|
| 22 | !C***
|
| 23 | !C
|
| 24 | module hecmw_solver_CG_33
|
| 25 |
|
| 26 | public :: hecmw_solve_CG_33
|
| 27 |
|
| 28 | contains
|
| 29 | !C
|
| 30 | !C*** CG_33
|
| 31 | !C
|
| 32 | subroutine hecmw_solve_CG_33( hecMESH, hecMAT, ITER, RESID, ERROR, &
|
| 33 | & Tset, Tsol, Tcomm )
|
| 34 |
|
| 35 | use hecmw_util
|
| 36 | use m_hecmw_solve_error
|
| 37 | use m_hecmw_comm_f
|
| 38 | use hecmw_matrix_misc
|
| 39 | use hecmw_solver_misc
|
| 40 | use hecmw_solver_las_33
|
| 41 | use hecmw_solver_scaling_33
|
| 42 | use hecmw_precond_33
|
| 43 | use hecmw_jad_type
|
| 44 | use hecmw_estimate_condition
|
| 45 |
|
| 46 | implicit none
|
| 47 |
|
| 48 | type(hecmwST_local_mesh) :: hecMESH
|
| 49 | type(hecmwST_matrix) :: hecMAT
|
| 50 | integer(kind=kint ), intent(inout):: ITER, ERROR
|
| 51 | real (kind=kreal), intent(inout):: RESID, Tset, Tsol, Tcomm
|
| 52 |
|
| 53 | integer(kind=kint ) :: N, NP, NDOF, NNDOF
|
| 54 | integer(kind=kint ) :: my_rank
|
| 55 | integer(kind=kint ) :: ITERlog, TIMElog
|
| 56 | real(kind=kreal), pointer :: B(:), X(:)
|
| 57 |
|
| 58 | real(kind=kreal), dimension(:,:), allocatable :: WW
|
| 59 |
|
| 60 | integer(kind=kint), parameter :: R= 1
|
| 61 | integer(kind=kint), parameter :: Z= 2
|
| 62 | integer(kind=kint), parameter :: Q= 2
|
| 63 | integer(kind=kint), parameter :: P= 3
|
| 64 | integer(kind=kint), parameter :: WK= 4
|
| 65 |
|
| 66 | integer(kind=kint ) :: MAXIT
|
| 67 |
|
| 68 | ! local variables
|
| 69 | real (kind=kreal) :: TOL
|
| 70 | integer(kind=kint )::i
|
| 71 | real (kind=kreal)::S_TIME,S1_TIME,E_TIME,E1_TIME, START_TIME, END_TIME
|
| 72 | real (kind=kreal)::BNRM2
|
| 73 | real (kind=kreal)::RHO,RHO1,BETA,C1,ALPHA,DNRM2
|
| 74 | real (kind=kreal)::t_max,t_min,t_avg,t_sd
|
| 75 | integer(kind=kint) ::ESTCOND
|
| 76 | real (kind=kreal), allocatable::D(:),E(:)
|
| 77 | real (kind=kreal)::ALPHA1
|
| 78 | integer(kind=kint) :: n_indef_precond
|
| 79 |
|
| 80 | integer(kind=kint), parameter :: N_ITER_RECOMPUTE_R= 50
|
| 81 |
|
| 82 | call hecmw_barrier(hecMESH)
|
| 83 | S_TIME= HECMW_WTIME()
|
| 84 |
|
| 85 | !C===
|
| 86 | !C +-------+
|
| 87 | !C | INIT. |
|
| 88 | !C +-------+
|
| 89 | !C===
|
| 90 | N = hecMAT%N
|
| 91 | NP = hecMAT%NP
|
| 92 | NDOF = hecMAT%NDOF
|
| 93 | NNDOF = N * NDOF
|
| 94 | my_rank = hecMESH%my_rank
|
| 95 | X => hecMAT%X
|
| 96 | B => hecMAT%B
|
| 97 |
|
| 98 | ITERlog = hecmw_mat_get_iterlog( hecMAT )
|
| 99 | TIMElog = hecmw_mat_get_timelog( hecMAT )
|
| 100 | MAXIT = hecmw_mat_get_iter( hecMAT )
|
| 101 | TOL = hecmw_mat_get_resid( hecMAT )
|
| 102 | ESTCOND = hecmw_mat_get_estcond( hecMAT )
|
| 103 |
|
| 104 | ERROR = 0
|
| 105 | n_indef_precond = 0
|
| 106 |
|
| 107 | allocate (WW(NDOF*NP, 4))
|
| 108 | WW = 0.d0
|
| 109 |
|
| 110 | !C
|
| 111 | !C-- SCALING
|
| 112 | call hecmw_solver_scaling_fw_33(hecMESH, hecMAT, Tcomm)
|
| 113 |
|
| 114 | IF (hecmw_mat_get_usejad(hecMAT).ne.0) THEN
|
| 115 | call hecmw_JAD_INIT(hecMAT)
|
| 116 | ENDIF
|
| 117 |
|
| 118 | if (ESTCOND /= 0 .and. hecMESH%my_rank == 0) then
|
| 119 | allocate(D(MAXIT),E(MAXIT-1))
|
| 120 | endif
|
| 121 | !C===
|
| 122 | !C +----------------------+
|
| 123 | !C | SETUP PRECONDITIONER |
|
| 124 | !C +----------------------+
|
| 125 | !C===
|
| 126 | call hecmw_precond_33_setup(hecMAT, hecMESH, 1)
|
| 127 |
|
| 128 | !C===
|
| 129 | !C +---------------------+
|
| 130 | !C | {r0}= {b} - [A]{x0} |
|
| 131 | !C +---------------------+
|
| 132 | !C===
|
| 133 | call hecmw_matresid_33(hecMESH, hecMAT, X, B, WW(:,R), Tcomm)
|
| 134 |
|
| 135 | !C-- compute ||{b}||
|
| 136 | call hecmw_InnerProduct_R(hecMESH, NDOF, B, B, BNRM2, Tcomm)
|
| 137 | if (BNRM2.eq.0.d0) then
|
| 138 | iter = 0
|
| 139 | MAXIT = 0
|
| 140 | RESID = 0.d0
|
| 141 | X = 0.d0
|
| 142 | endif
|
| 143 |
|
| 144 | E_TIME = HECMW_WTIME()
|
| 145 | if (TIMElog.eq.2) then
|
| 146 | call hecmw_time_statistics(hecMESH, E_TIME - S_TIME, &
|
| 147 | t_max, t_min, t_avg, t_sd)
|
| 148 | if (hecMESH%my_rank.eq.0) then
|
| 149 | write(*,*) 'Time solver setup'
|
| 150 | write(*,*) ' Max :',t_max
|
| 151 | write(*,*) ' Min :',t_min
|
| 152 | write(*,*) ' Avg :',t_avg
|
| 153 | write(*,*) ' Std Dev :',t_sd
|
| 154 | endif
|
| 155 | Tset = t_max
|
| 156 | else
|
| 157 | Tset = E_TIME - S_TIME
|
| 158 | endif
|
| 159 |
|
| 160 | Tcomm = 0.d0
|
| 161 | call hecmw_barrier(hecMESH)
|
| 162 | S1_TIME = HECMW_WTIME()
|
| 163 | !C
|
| 164 | !C************************************************* Conjugate Gradient Iteration start
|
| 165 | !C
|
| 166 | do iter = 1, MAXIT
|
| 167 |
|
| 168 | !C===
|
| 169 | !C +----------------+
|
| 170 | !C | {z}= [Minv]{r} |
|
| 171 | !C +----------------+
|
| 172 | !C===
|
| 173 | call hecmw_precond_33_apply(hecMESH, hecMAT, WW(:,R), WW(:,Z), WW(:,WK), Tcomm)
|
| 174 |
|
| 175 | !C===
|
| 176 | !C +---------------+
|
| 177 | !C | {RHO}= {r}{z} |
|
| 178 | !C +---------------+
|
| 179 | !C===
|
| 180 | call hecmw_InnerProduct_R(hecMESH, NDOF, WW(:,R), WW(:,Z), RHO, Tcomm)
|
| 181 | ! if RHO is NaN or Inf then no converge
|
| 182 | if (RHO == 0.d0) then
|
| 183 | ! converged due to RHO==0
|
| 184 | exit
|
| 185 | elseif (iter > 1 .and. RHO*RHO1 <= 0) then
|
| 186 | n_indef_precond = n_indef_precond + 1
|
| 187 | if (n_indef_precond >= 3) then
|
| 188 | ! diverged due to indefinite preconditioner
|
| 189 | ERROR = HECMW_SOLVER_ERROR_DIVERGE_PC
|
| 190 | exit
|
| 191 | endif
|
| 192 | endif
|
| 193 |
|
| 194 | !C===
|
| 195 | !C +-----------------------------+
|
| 196 | !C | {p} = {z} if ITER=1 |
|
| 197 | !C | BETA= RHO / RHO1 otherwise |
|
| 198 | !C +-----------------------------+
|
| 199 | !C===
|
| 200 | if ( ITER.eq.1 ) then
|
| 201 | do i = 1, NNDOF
|
| 202 | WW(i,P) = WW(i,Z)
|
| 203 | enddo
|
| 204 | else
|
| 205 | BETA = RHO / RHO1
|
| 206 | do i = 1, NNDOF
|
| 207 | WW(i,P) = WW(i,Z) + BETA*WW(i,P)
|
| 208 | enddo
|
| 209 | endif
|
| 210 |
|
| 211 | !C===
|
| 212 | !C +--------------+
|
| 213 | !C | {q}= [A] {p} |
|
| 214 | !C +--------------+
|
| 215 | !C===
|
| 216 | call hecmw_matvec_33(hecMESH, hecMAT, WW(:,P), WW(:,Q), Tcomm)
|
| 217 |
|
| 218 | !C===
|
| 219 | !C +---------------------+
|
| 220 | !C | ALPHA= RHO / {p}{q} |
|
| 221 | !C +---------------------+
|
| 222 | !C===
|
| 223 | call hecmw_InnerProduct_R(hecMESH, NDOF, WW(:,P), WW(:,Q), C1, Tcomm)
|
| 224 | ! if C1 is NaN or Inf, no converge
|
| 225 | if (C1 <= 0) then
|
| 226 | ! diverged due to indefinite or negative definite matrix
|
| 227 | ERROR = HECMW_SOLVER_ERROR_DIVERGE_MAT
|
| 228 | exit
|
| 229 | endif
|
| 230 |
|
| 231 | ALPHA= RHO / C1
|
| 232 |
|
| 233 | !C===
|
| 234 | !C +----------------------+
|
| 235 | !C | {x}= {x} + ALPHA*{p} |
|
| 236 | !C | {r}= {r} - ALPHA*{q} |
|
| 237 | !C +----------------------+
|
| 238 | !C===
|
| 239 | do i = 1, NNDOF
|
| 240 | X(i) = X(i) + ALPHA * WW(i,P)
|
| 241 | ! WW(i,R)= WW(i,R) - ALPHA * WW(i,Q)
|
| 242 | enddo
|
| 243 |
|
| 244 | if ( mod(ITER,N_ITER_RECOMPUTE_R)==0 ) then
|
| 245 | call hecmw_matresid_33(hecMESH, hecMAT, X, B, WW(:,R), Tcomm)
|
| 246 | else
|
| 247 | do i = 1, NNDOF
|
| 248 | WW(i,R)= WW(i,R) - ALPHA * WW(i,Q)
|
| 249 | enddo
|
| 250 | endif
|
| 251 |
|
| 252 | call hecmw_InnerProduct_R(hecMESH, NDOF, WW(:,R), WW(:,R), DNRM2, Tcomm)
|
| 253 |
|
| 254 | RESID= dsqrt(DNRM2/BNRM2)
|
| 255 |
|
| 256 | !C##### ITERATION HISTORY
|
| 257 | if (my_rank.eq.0.and.ITERLog.eq.1) write (*,'(i7, 1pe16.6)') ITER, RESID
|
| 258 | !C#####
|
| 259 |
|
| 260 | if (ESTCOND /= 0 .and. hecMESH%my_rank == 0) then
|
| 261 | if (ITER == 1) then
|
| 262 | D(1) = 1.d0 / ALPHA
|
| 263 | else
|
| 264 | D(ITER) = 1.d0 / ALPHA + BETA / ALPHA1
|
| 265 | E(ITER-1) = sqrt(BETA) / ALPHA1
|
| 266 | endif
|
| 267 | ALPHA1 = ALPHA
|
| 268 | if (mod(ITER,ESTCOND) == 0) call hecmw_estimate_condition_CG(ITER, D, E)
|
| 269 | endif
|
| 270 |
|
| 271 | if ( RESID.le.TOL ) then
|
| 272 | if ( mod(ITER,N_ITER_RECOMPUTE_R)==0 ) exit
|
| 273 | !C----- recompute R to make sure it is really converged
|
| 274 | call hecmw_matresid_33(hecMESH, hecMAT, X, B, WW(:,R), Tcomm)
|
| 275 | call hecmw_InnerProduct_R(hecMESH, NDOF, WW(:,R), WW(:,R), DNRM2, Tcomm)
|
| 276 | RESID= dsqrt(DNRM2/BNRM2)
|
| 277 | if ( RESID.le.TOL ) exit
|
| 278 | endif
|
| 279 | if ( ITER .eq.MAXIT ) ERROR = HECMW_SOLVER_ERROR_NOCONV_MAXIT
|
| 280 |
|
| 281 | RHO1 = RHO
|
| 282 |
|
| 283 | enddo
|
| 284 | !C
|
| 285 | !C************************************************* Conjugate Gradient Iteration end
|
| 286 | !C
|
| 287 | call hecmw_solver_scaling_bk_33(hecMAT)
|
| 288 | !C
|
| 289 | !C-- INTERFACE data EXCHANGE
|
| 290 | !C
|
| 291 | START_TIME= HECMW_WTIME()
|
| 292 | call hecmw_update_3_R (hecMESH, X, hecMAT%NP)
|
| 293 | END_TIME = HECMW_WTIME()
|
| 294 | Tcomm = Tcomm + END_TIME - START_TIME
|
| 295 |
|
| 296 | deallocate (WW)
|
| 297 | !call hecmw_precond_33_clear(hecMAT)
|
| 298 |
|
| 299 | IF (hecmw_mat_get_usejad(hecMAT).ne.0) THEN
|
| 300 | call hecmw_JAD_FINALIZE()
|
| 301 | ENDIF
|
| 302 |
|
| 303 | if (ESTCOND /= 0 .and. ERROR == 0 .and. hecMESH%my_rank == 0) then
|
| 304 | call hecmw_estimate_condition_CG(ITER, D, E)
|
| 305 | deallocate(D, E)
|
| 306 | endif
|
| 307 |
|
| 308 | E1_TIME = HECMW_WTIME()
|
| 309 | if (TIMElog.eq.2) then
|
| 310 | call hecmw_time_statistics(hecMESH, E1_TIME - S1_TIME, &
|
| 311 | t_max, t_min, t_avg, t_sd)
|
| 312 | if (hecMESH%my_rank.eq.0) then
|
| 313 | write(*,*) 'Time solver iterations'
|
| 314 | write(*,*) ' Max :',t_max
|
| 315 | write(*,*) ' Min :',t_min
|
| 316 | write(*,*) ' Avg :',t_avg
|
| 317 | write(*,*) ' Std Dev :',t_sd
|
| 318 | endif
|
| 319 | Tsol = t_max
|
| 320 | else
|
| 321 | Tsol = E1_TIME - S1_TIME
|
| 322 | endif
|
| 323 |
|
| 324 | end subroutine hecmw_solve_CG_33
|
| 325 |
|
| 326 | end module hecmw_solver_CG_33
|
