1!********************************************************* 2! ramg control 3! Control AMG preparation/solve process 4! will have dynamic control (based on solver vecters) 5!********************************************************* 6 7 subroutine ramg_control(colm,rowp,lhsK,lhsP, 8 & ilwork,BC,iBC,iper) 9 use ramg_data 10 include "common.h" 11 12 integer,intent(in),dimension(nshg+1) :: colm 13 integer,intent(in),dimension(nnz_tot) :: rowp 14 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 15 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 16 ! the boundary info 17 integer, intent(in), dimension(nlwork) :: ilwork 18 integer, intent(in),dimension(nshg) :: iBC,iper 19 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 20 21 integer i 22 23 if (iamg_init .eq. 0 ) then 24 ! The overall initialization 25 ! do init 26 ramg_winct = 0 27 ramg_setup_flag = 0 28 iamg_init = 1 29 ramg_time = 0 30 else 31 ramg_setup_flag = mod(ramg_setup_flag+1 ,iamg_setup_frez) 32 end if 33 34 ! Extract PPE 35 call ramg_extract_ppe(colm,rowp,lhsK,lhsP, 36 & ilwork,BC,iBC,iper) 37 ! Coarsening 38 call ramg_prep(ilwork,BC,iBC,iper) 39 call ramg_init_ilwork(ilwork,BC,iBC,iper) 40 ! Prepare for MLS smoothing 41 if (iamg_smoother.eq.2) then 42 call ramg_cheby_setup(colm,rowp,lhsK,lhsP, 43 & ilwork,BC,iBC,iper) 44 else 45 call ramg_mls_setup(colm,rowp,lhsK,lhsP, 46 & ilwork,BC,iBC,iper) 47 endif 48 if (iamg_c_solver.eq.2) then 49 ! Setup Coarsest Direct solver 50 call ramg_direct_LU(amg_A_colm(ramg_levelx)%p, 51 & amg_A_rowp(ramg_levelx)%p, 52 & amg_A_lhs(ramg_levelx)%p,amg_nshg(ramg_levelx), 53 & amg_nnz(ramg_levelx)) 54 endif 55 if (maxnev.gt.0) then 56 ! Setup GGB 57 call ramg_ggb_setup(colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 58 endif 59 60 61 end subroutine !ramg_control 62 63!********************************************************* 64! Ramg preparation 65!********************************************************* 66 67 subroutine ramg_prep(ilwork,BC,iBC,iper) 68 use ramg_data 69 include "common.h" 70 include "mpif.h" 71 include "auxmpi.h" 72 73 integer,intent(in),dimension(nlwork) :: ilwork 74 integer,intent(in),dimension(nshg) :: iBC,iper 75 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 76 77 logical :: maxstopsign,mxs2 78 integer :: i,p,p2 79 80 if (ramg_setup_flag.ne.0) return 81 maxstopsign = .false. 82 i = 1 83 do while ((i+1.le.iamg_nlevel) .and. (.not.maxstopsign)) 84 if (amg_nshg(i+1).ne.0) then 85 call ramg_deallocate(i+1) 86 end if 87 call ramg_coarse_setup(i,i+1,strong_eps,iamg_interp, 88 & ramg_trunc, 89 & ilwork,BC,iBC,iper,nshg,nlwork,ndofBC) 90 call ramg_calcITAI(i,i+1,maxstopsign) 91 if ((iamg_verb.gt.2).and.(myrank.eq.master)) then 92 write(*,*)'COARSEN: level:',i+1,' nshg:',amg_nshg(i+1), 93 & ' nnz:',amg_nnz(i+1) 94 endif 95 i = i+1 96 call MPI_Barrier(MPI_COMM_WORLD,ierr) 97 maxstopsign = (amg_nshg(i).eq.amg_nshg(i-1)) 98 IF (.true.) THEN 99 !call ramg_checkcoarse(i,ilwork,BC,iBC,iper,maxstopsign) 100 p = 1 101 if (maxstopsign) p = 0 102 call MPI_AllReduce(p,p2,1,MPI_INTEGER,MPI_SUM, 103 & MPI_COMM_WORLD,ierr) 104 if (p2.eq.0) then 105 maxstopsign = .true. 106 !write(*,*)'mcheck stopped' 107 else 108 maxstopsign = .false. 109 endif 110 ENDIF 111 enddo 112 ramg_levelx=i 113 if (maxstopsign) then 114 ramg_levelx = ramg_levelx-1 115 endif 116 117! if ( (irun_amg_prec.eq.1).and.(mlsDeg.gt.0) ) then 118 if (.false.) then 119 deallocate(amg_A_colm(1)%p) 120 deallocate(amg_A_rowp(1)%p) 121 deallocate(amg_A_lhs(1)%p) 122 endif 123 124 allocate(CF_map(ramg_levelx)%p(amg_nshg(ramg_levelx))) 125 allocate(CF_revmap(ramg_levelx)%p(amg_nshg(ramg_levelx))) 126 do i=1,amg_nshg(ramg_levelx) 127 CF_map(ramg_levelx)%p(i) = i 128 CF_revmap(ramg_levelx)%p(i) = i 129 enddo 130 131 !call ramg_output_coarsening 132 133 end subroutine !ramg_prep 134 135 subroutine ramg_checkcoarse(level1,ilwork,BC,iBC,iper, 136 & cfstop) 137 use ramg_data 138 include "common.h" 139 include "mpif.h" 140 include "auxmpi.h" 141 142 integer,intent(in) :: level1 143 integer,intent(in),dimension(nlwork) :: ilwork 144 integer,intent(in),dimension(nshg) :: iBC,iper 145 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 146 147 logical,intent(inout) :: cfstop 148 integer :: i,numtask,itkbeg,numseg,iacc 149 integer :: j,k,p 150 integer,allocatable,dimension(:) 151 & :: subcfstop,subnnz,subcf,subcfrev,subnei 152 real(kind=8) :: rhoratio 153 logical :: subneireduce 154 155 IF (( iamg_reduce.le.1).or.(numpe.gt.1)) THEN 156! IF (.TRUE.) THEN 157 if (numpe.ge.2) then 158 numtask = ilwork(1)+1 159 else 160 numtask = 1 161 endif 162 163 allocate(subcfrev(numpe)) 164 subcfrev = 0 165 allocate(subnei(numtask)) 166 subnei = 0 167 168 subnei(1) = myrank+1 169 subcfrev(subnei(1))=1 170 itkbeg = 1 171 do i=2,numtask 172 subnei(i)=ilwork(itkbeg+3)+1 173 subcfrev(subnei(i)) = i 174 itkbeg = itkbeg + 4 + 2*ilwork(itkbeg+4) 175 enddo 176 177 ELSE !reduced 178 numtask = rmapmax-1 179 allocate(subcfrev(numtask)) 180 allocate(subnei(numtask)) 181 do i=1,numtask 182 subcfrev(i) = i 183 enddo 184 ENDIF 185 186 allocate(subcfstop(numtask)) 187 allocate(subcf(numtask)) 188 allocate(subnnz(numtask)) 189 subcfstop = 0 190 subcf = 0 191 subnnz = 0 192 193 194 do i = 1,amg_nnz(level1) 195 k = amg_A_rowp(level1)%p(i) 196 p = iabs(amg_paramap(level1)%p(k)) 197! if (iamg_reduce.gt.1) then 198! p = 1 ! for reduced only 199! endif 200 p = subcfrev(p) 201 subnnz(p) = subnnz(p) + 1 202 enddo 203 do i = 1,amg_nshg(level1) 204 p = iabs(amg_paramap(level1)%p(i)) 205! if (iamg_reduce.gt.1) then 206! p = 1 ! for reduced only 207! endif 208 p = subcfrev(p) 209 subcf(p) = subcf(p) + 1 210 enddo 211 212 do i=1,numtask 213 IF ((iamg_reduce.le.1).or.(numpe.gt.1)) THEN 214! IF (.TRUE.) THEN 215 p = subnei(i) 216 subneireduce = (p.ne.(myrank+1)) 217 ELSE !reduced 218 subneireduce = (i.gt.iamg_reduce) 219 END IF 220 if ((subcf(i).lt.30).and.subneireduce) then 221 subcfstop(i) = 1 222 endif 223 if (.not.subneireduce) then 224 if (subcf(i).lt.200) then 225 subcfstop(i) = 1 226 endif 227 if ((subnnz(i)/(subcf(i)**2)).gt.0.6) then 228 subcfstop(i) = 1 229 endif 230 end if 231 enddo 232 233 cfstop = .true. 234 do i=1,numtask 235 if (subcfstop(i).eq.0) then 236 cfstop = .false. 237 exit 238 endif 239 enddo 240 241 242 IF ((iamg_reduce.le.1).or.(numpe.gt.1)) THEN 243 ! if interior is coarsened, boundary should be fixed whatever 244 if (subcfstop(1).eq.1) then 245 cfstop = .true. 246 endif 247 ENDIF 248 249 if (.not.cfstop) then 250 do i=1,amg_nshg(level1) 251 p = iabs(amg_paramap(level1)%p(i)) 252! if (iamg_reduce.gt.1) then 253! p = 1 ! reduced only 254! endif 255 p = subcfrev(p) 256 if (subcfstop(p).eq.1) then 257 amg_paramap(level1)%p(i) = -iabs(amg_paramap(level1)%p(i)) 258 endif 259 enddo 260 endif 261 262 deallocate(subcfrev) 263 deallocate(subnei) 264 deallocate(subcfstop) 265 deallocate(subcf) 266 deallocate(subnnz) 267 268 end subroutine ! ramg_checkcoarse 269 270!****************************************************************** 271! A bunch of code that hook to leslib's cg solve 272! Force it to restart with AMG if meets plateau 273!****************************************************************** 274 275!******** 276! sc writes: 277! irun_amg_prec is a variable controling how to use AMG wisely. 278! 0 : no run 279! 1 : always use AMG prec'd CG 280! 2 : restart CG a ( if plain CG hits plateau, do again with AMG ) 281! 3 : restart CG b ( if plain exceeds maxiter, do again with AMG ) 282! also refer to input.config for a detailed description. 283!******** 284 285 286 subroutine ramg_normcheck(tmpnorm) 287 use ramg_data 288 289 include "common.h" 290 include "mpif.h" 291 292 real(kind=8),intent(inout) :: tmpnorm 293 real(kind=8) :: sqnorm 294 295 !write(*,*)myrank,' winct:',ramg_winct,sqrt(tmpnorm) 296 !return 297 if (irun_amg_prec.ne.2) return ! No control at all 298 if (ramg_winct.eq.0) then 299 !ramg_window = sqrt(tmpnorm) 300 return ! Not relevant 301 endif 302 ramg_winct=0 303 if (ramg_flag.gt.1) return ! Second run does not need control 304 305 ! Following are for the first run 306 !write(*,*)'normcheck: ',sqnorm 307 !return; 308 309 if (ramg_redo.ge.100) then 310 ! The rest of GMRES should be terminated 311 if (ramg_redo.lt.104) then 312 ramg_winct = 1 313 ramg_redo=ramg_redo+1 314! write(*,*)'***' 315 else 316 ramg_redo = 0 317 tmpnorm = 0 318 endif 319 return 320 endif 321 sqnorm = sqrt(tmpnorm) 322 323 if (ramg_redo.le.7) then 324 ramg_redo = ramg_redo + 1 325 call ramg_winpushin(sqnorm) 326 return 327 endif 328 329 call ramg_winpushin(sqnorm) 330 331 if (ramg_window(7).gt.ramg_window(1)) then ! Stop point 332! write(*,*)'myrank:',myrank,ramg_window(7),ramg_window(1) 333! if (ramg_window(7).lt.1e-3) then 334 ! stop point in cg 335 if (myrank.eq.master) then 336 write(*,*)"Prepare to restart CG" 337 endif 338 tmpnorm = 0 339 ramg_redo = 100 340 ramg_winct = 1 341 ramg_flag = ramg_flag -1 342 endif 343 344 end subroutine ! ramg_normcheck 345 346 subroutine ramg_winpushin(sqnorm) 347 use ramg_data 348 include "common.h" 349 include "mpif.h" 350 351 real(kind=8),intent(in) :: sqnorm 352 integer i 353 integer totlen 354 355 totlen = 7 356 357 do i=1,totlen-1 358 ramg_window(i)=ramg_window(i+1) 359 enddo 360 ramg_window(totlen) = sqnorm 361 362 end subroutine ! ramg_winpushin 363 364 365!********************************************************* 366! <EOF> RAMG Control 367!********************************************************* 368