1!************************************************************* 2! RAMG Coarse: 3! do C/F spliting and setup coarse matrix 4!************************************************************** 5! ramg: ramg_coarse_setup 6! Input: amg_A matrix level1 7! Output: I_h_H matrix level2 to level1 8! I_H_h matrix level1 to level2 9!************************************************************** 10 subroutine ramg_coarse_setup(level1,level2,eps_str, 11 & interp_flag,trunc_tol, 12 & ilwork,BC,iBC,iper)!,nshg,nlwork,ndofBC) 13 use ramg_data 14 include "common.h" 15 include "mpif.h" 16 include "auxmpi.h" 17 18! implicit none 19 20 !******* parameters ******* 21 integer,intent(in) :: level1,level2 22 real(kind=8),intent(in) :: eps_str 23 integer,intent(in) :: interp_flag 24 real(kind=8),intent(in) :: trunc_tol 25! integer,intent(in) :: nshg,nlwork,ndofBC 26 integer,intent(in),dimension(nlwork) :: ilwork 27 integer,intent(in),dimension(nshg) :: iBC,iper 28 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 29 !******* temp variables ******* 30 integer,dimension(amg_nnz(level1)) :: amg_S,amg_Ip 31 integer,dimension(amg_nshg(level1)) :: amg_F,aLoc 32 real(kind=8),dimension(amg_nshg(level1)) :: amg_la 33 integer,dimension(amg_nshg(level1)) :: amg_Fn,amg_Fp 34 real(kind=8) :: alpha,beta,alphac,betac,diag,rtp 35 real(kind=8) :: tmp1,tmp2,tmp3,tmp4 36 ! AMG_F : 0: UNDECIDED, 1: COARSE, 2: FINE 37 integer :: I_nnz,ki,kj,jj 38 integer :: i,j,k,m,n,p,q 39 character :: fname*80 40 41 integer :: numtask,itkbeg,numseg,iacc 42 integer,allocatable,dimension(:) :: subcf,subcfrev,subnei 43 44 integer :: mnnz 45 integer :: cfilter 46 type(i2dd) :: amg_I_rowp 47 type(r2dd) :: amg_I 48 type(i1d) :: amg_I_colm 49 50 51 real,dimension(10) :: cpusec 52 logical :: ti_flag 53 integer :: mem_err,mem_err_s 54 !******* end of variables ******** 55 56 call cpu_time(cpusec(1)) 57 58 if (abs(trunc_tol).lt.1e-4) then 59 ti_flag = .false. 60 else 61 ti_flag = .true. 62 end if 63 64 amg_S = 0 65 amg_F = 0 66 67 if (numpe.ge.2) then 68 numtask = ilwork(1)+1 69 else 70 numtask = 1 71 endif 72 73 IF ((iamg_reduce.le.1).or.(numpe.gt.1)) THEN 74! IF (.TRUE.) THEN 75 allocate(subcfrev(numpe)) 76 subcfrev = 0 77 allocate(subcf(numtask)) 78 subcf = 0 79 allocate(subnei(numtask)) 80 subnei = 0 81 82 subnei(1) = myrank+1 83 subcfrev(subnei(1))=1 84 itkbeg = 1 85 do i=2,numtask 86 subnei(i)=ilwork(itkbeg+3)+1 ! iother+1 87 subcfrev(subnei(i)) = i 88 itkbeg = itkbeg + 4 + 2*ilwork(itkbeg+4) 89 enddo 90 91 ELSE !reduced 92 numtask = rmapmax-1 93 allocate(subcfrev(numtask)) 94 allocate(subcf(numtask)) 95 allocate(subnei(numtask)) 96 subcfrev = 0 97 subcf = 0 98 subnei = 0 99 do i=1,numtask 100 subcfrev(i) = i 101 enddo 102 END IF 103 104 call ramg_CFsplit(amg_F,amg_S,amg_nshg(level1),amg_nnz(level1), 105 & amg_A_colm(level1)%p,amg_A_rowp(level1)%p, 106 & amg_A_lhs(level1)%p,amg_paramap(level1)%p, 107 & ilwork,BC,iBC,iper, 108 & eps_str,-1,0) 109 ! read FC spliting to file 110 if (level1.eq.1) then 111 ! write(fname,'((A7)(I1))')'GTG_FC_',level1 112 ! call ramg_readin_i(amg_F,amg_nshg(level1),fname) 113 ! write FC spliting to file 114 !write(fname,'((A7)(I1))')'amg_FC_',level1 115 !call ramg_dump_i(amg_F,amg_nshg(level1),1,fname) 116 endif 117 118 call ramg_update_cfmap(amg_F,level1) 119 ! communication goes here 120 call MPI_Barrier(MPI_COMM_WORLD,ierr) 121 call commOut_i(amg_cfmap,ilwork,1,iper,iBC,BC) 122 call ramg_readin_cfmap(amg_F,level1) 123 124 IF (.true.) THEN 125 subcf = 0 126 subnei = 0 127 do i = 1,amg_nshg(level1) 128 p = iabs(amg_paramap(level1)%p(i)) 129 p = subcfrev(p) 130 subnei(p) = subnei(p) + 1 131 if (amg_F(i).eq.1) then 132 subcf(p) = subcf(p) + 1 133 endif 134 enddo 135 do i = 1,numtask 136 enddo 137 ENDIF 138 139 call cpu_time(cpusec(3)) 140 141 deallocate(subcf) 142 deallocate(subcfrev) 143 deallocate(subnei) 144 145 ! SETUP Smoothing ordering 146 allocate(CF_map(level1)%p(amg_nshg(level1))) 147 allocate(CF_revmap(level1)%p(amg_nshg(level1))) 148 149 k = 1 150 do i=1,amg_nshg(level1) 151 if (amg_F(i).eq.1) then ! fine/coarse first 152 CF_map(level1)%p(k) = i 153 CF_revmap(level1)%p(i) = k 154 k = k+1 155 end if 156 enddo 157 do i=1,amg_nshg(level1) 158 if (amg_F(i).eq.2) then 159 CF_map(level1)%p(k) = i 160 CF_revmap(level1)%p(i) = k 161 k = k+1 162 end if 163 enddo 164 165 if (level1.eq.-1) then 166 do i=1,amg_nshg(level1) 167 CF_map(level1)%p(i) = i 168 CF_revmap(level1)%p(i) = i 169 enddo 170 endif 171 172 ! generate I := I_H_h, from coarse to fine 173 ! then TRANSPOSE and get I^T, from fine to coarse 174 175 ! determine the size of second level 176 amg_nshg(level2) = 0 177 amg_Ip = 0 178 I_nnz = 0 179 do i = 1, amg_nshg(level1) 180 if (amg_F(i).eq.1) then 181 amg_nshg(level2) = amg_nshg(level2) + 1 182 I_nnz = I_nnz + 1 183 amg_Ip(i) = I_nnz 184 end if 185 enddo 186 187 allocate(amg_paramap(level2)%p(amg_nshg(level2))) 188 allocate(amg_paraext(level2)%p(amg_nshg(level2))) 189 j = 1 190 do i=1,amg_nshg(level1) 191 if (amg_F(i).eq.1) then 192 amg_paramap(level2)%p(j) = amg_paramap(level1)%p(i) 193 amg_paraext(level2)%p(j) = amg_paraext(level1)%p(i) 194 j = j + 1 195 end if 196 enddo 197 198 ! ALLOCATE I_H_h from coarse to fine 199 mem_err_s = 0 200 allocate(amg_I_colm%p(amg_nshg(level1)),stat=mem_err) 201 mem_err_s = mem_err_s + mem_err 202 allocate(amg_I_rowp%pp(amg_nshg(level1)),stat=mem_err) 203 mem_err_s = mem_err_s + mem_err 204 allocate(amg_I%pp(amg_nshg(level1)),stat=mem_err) 205 mem_err_s = mem_err_s + mem_err 206 if (mem_err_s.ne.0) then 207 write(*,*)'Memory allocation error while geting I_H_h' 208 end if 209 210 cpusec(4) = 0 211 cpusec(5) = 0 212 cpusec(8) = 0 213 mnnz = 0 214 aLoc = 0 215 loop_i: do i = 1, amg_nshg(level1) 216 call cpu_time(cpusec(6)) 217 cfilter = amg_paramap(level1)%p(i) 218 if (amg_F(i).eq.1) then ! i as a coarse node do trivial work 219 allocate(amg_I_rowp%pp(i)%p(1)) 220 allocate(amg_I%pp(i)%p(1)) 221 amg_I_colm%p(i) = 1 222 amg_I_rowp%pp(i)%p(1) = amg_Ip(i) 223 amg_I%pp(i)%p(1) = 1 224 mnnz = mnnz + 1 225 cycle loop_i 226 end if 227 n = 0 228 ! a-hat, P-hat, N-hat 229 ! amg_la ! a-hat 230 ! amg_Fn ! N-hat 231 ! amg_Fp ! P-hat 232 ! diagonal part 233 diag = -1.0d0/amg_A_lhs(level1)%p(amg_A_colm(level1)%p(i),1) 234 ! non-diagonal 235 do k = amg_A_colm(level1)%p(i)+1,amg_A_colm(level1)%p(i+1)-1 236 j = amg_A_rowp(level1)%p(k) 237 if (cfilter.eq.amg_paramap(level1)%p(j)) then 238 n = n + 1 239 if ( (amg_F(j).eq.1) .and. (mod(amg_S(k),2).eq.1) ) then 240 amg_Fp(n) = 1 241 else 242 amg_Fp(n) = 0 243 endif 244 amg_Fn(n) = j 245 amg_la(n) = amg_A_lhs(level1)%p(k,1) 246 aLoc(j) = n 247 end if 248 enddo 249 250 ! standard interpolation here ! for parallel dont use 251 if ( interp_flag .eq. 1) then 252 do k=amg_A_colm(level1)%p(i)+1,amg_A_colm(level1)%p(i+1)-1 253 j = amg_A_rowp(level1)%p(k) 254 if ( (amg_F(j).eq.2) .and. (mod(amg_S(k),2).eq.1)) then 255 rtp = amg_A_lhs(level1)%p(amg_A_colm(level1)%p(j),1) 256 if (rtp.ne.0) then 257 rtp = amg_la(aLoc(j))/rtp 258 endif 259 do kj=amg_A_colm(level1)%p(j),amg_A_colm(level1)%p(j+1)-1 260 jj = amg_A_rowp(level1)%p(kj) 261 m = aLoc(jj) 262 if (m.eq.0) then 263 m = n+1 264 aLoc(jj) = m 265 amg_Fn(m) = jj 266 amg_Fp(m) = 0 267 amg_la(m) = 0 268 n = n+1 269 end if 270 amg_la(m) = amg_la(m) - rtp*amg_A_lhs(level1)%p(kj,1) 271 if ( (amg_F(jj).eq.1) .and. (mod(amg_S(kj),2).eq.1) ) then 272 amg_Fp(m) = 1 273 else 274 amg_Fp(m) = 0 275 endif 276 enddo 277 end if 278 enddo 279 endif 280 281 call cpu_time(cpusec(7)) 282 cpusec(4) = cpusec(4) + cpusec(7)-cpusec(6) 283 ! calculate alpha, beta, nnz-in-row 284 alpha = 0 285 beta = 0 286 alphac = 0 287 betac = 0 288 I_nnz = 0 289 do k = 1,n ! all non diagonal 290 aLoc(amg_Fn(k))=0 291 ! in N-hat 292 if ( amg_la(k).lt.0 ) then ! a-hat < 0 293 alpha = alpha + amg_la(k) 294 else ! a-hat > 0 295 beta = beta + amg_la(k) 296 end if 297 ! in P-hat 298 if (amg_Fp(k).eq.1) then 299 I_nnz = I_nnz + 1 ! nonzero in I_h_H row 300 if (amg_la(k).lt.0) then ! a-hat < 0 in coarse 301 alphac = alphac + amg_la(k) 302 else ! a-hat > 0 in coarse 303 betac = betac + amg_la(k) 304 end if 305 amg_Fn(I_nnz) = amg_Fn(k) 306 amg_la(I_nnz) = amg_la(k) 307 end if 308 enddo 309 alpha = diag*alpha/alphac 310 beta = diag*beta/betac 311 rtp = 0 312 do k = 1,I_nnz 313 if ( amg_la(k).lt.0) then 314 amg_la(k) = amg_la(k)*alpha 315 else 316 amg_la(k) = amg_la(k)*beta 317 end if 318 rtp = max(rtp,abs(amg_la(k))) 319 enddo 320 if (ti_flag) then 321 n = I_nnz 322 I_nnz = 0 323 tmp1 = 0 324 tmp2 = 0 325 tmp3 = 0 326 tmp4 = 0 327 rtp = rtp*trunc_tol 328 do k=1,n 329 tmp3 = tmp3 + min(amg_la(k),zero) 330 tmp4 = tmp4 + max(amg_la(k),zero) 331 if (abs(amg_la(k)).gt.rtp) then 332 tmp1 = tmp1 + min(amg_la(k),zero) 333 tmp2 = tmp2 + max(amg_la(k),zero) 334 I_nnz = I_nnz + 1 335 amg_Fn(I_nnz) = amg_Fn(k) 336 amg_la(I_nnz) = amg_la(k) 337 end if 338 enddo 339 if ( abs(tmp1).gt.1e-5) tmp1 = tmp3/tmp1 340 if ( abs(tmp2).gt.1e-5) tmp2 = tmp4/tmp2 341 do k=1,I_nnz 342 if (amg_la(k).lt.0) then 343 amg_la(k) = amg_la(k)*tmp1 344 else 345 amg_la(k) = amg_la(k)*tmp2 346 end if 347 enddo 348 end if 349 ! put up I matrix row 350 call cpu_time(cpusec(6)) 351 cpusec(5) = cpusec(5) + cpusec(6)-cpusec(7) 352 mnnz = mnnz + I_nnz 353 amg_I_colm%p(i) = I_nnz 354 allocate(amg_I_rowp%pp(i)%p(I_nnz)) 355 allocate(amg_I%pp(i)%p(I_nnz)) 356 do k = 1,I_nnz ! scan over P-hat 357 amg_I_rowp%pp(i)%p(k) = amg_Ip(amg_Fn(k)) 358 ! put in rowp 359 if (amg_la(k).lt.0) then 360 amg_I%pp(i)%p(k) = amg_la(k) 361 else 362 amg_I%pp(i)%p(k) = amg_la(k) 363 endif 364 enddo 365 call cpu_time(cpusec(7)) 366 cpusec(8) = cpusec(8) + cpusec(7) - cpusec(6) 367 enddo loop_i ! Now we have I_H_h 368 369 ! copy the data 370 mem_err_s = 0 371 allocate(I_cf_colm(level1)%p(amg_nshg(level1)+1),stat=mem_err) 372 mem_err_s = mem_err_s + mem_err 373 allocate(I_cf_rowp(level1)%p(mnnz),stat=mem_err) 374 mem_err_s = mem_err_s + mem_err 375 allocate(I_cf(level1)%p(mnnz),stat=mem_err) 376 mem_err_s = mem_err_s + mem_err 377 if (mem_err_s .ne. 0 ) then 378 write(*,*)'allocation error in I_cf' 379 end if 380 mnnz = 0 381 do i=1,amg_nshg(level1) 382 I_cf_colm(level1)%p(i)=mnnz+1 383 do j=1,amg_I_colm%p(i) 384 I_cf_rowp(level1)%p(mnnz+j)=amg_I_rowp%pp(i)%p(j) 385 I_cf(level1)%p(mnnz+j)=amg_I%pp(i)%p(j) 386 enddo 387 mnnz = mnnz + amg_I_colm%p(i) 388 enddo 389 I_cf_colm(level1)%p(amg_nshg(level1)+1) = mnnz+1 390 391 ! dump interpolator 392! if (level1.eq.1) then 393! write(fname,'((A4)(I1))')'Icf_',level1 394! write(*,*)amg_nshg(level1),mnnz 395! call ramg_dump_matlab_A(I_cf_colm(level1)%p,I_cf_rowp(level1)%p, 396! & I_cf(level1)%p,amg_nshg(level1),mnnz,1,fname) 397! endif 398 399 mem_err_s = 0 400 do i = 1,amg_nshg(level1) 401 if (amg_I_colm%p(i).ne.0) then 402 deallocate(amg_I_rowp%pp(i)%p,stat=mem_err) 403 mem_err_s = mem_err_s + mem_err 404 deallocate(amg_I%pp(i)%p,stat=mem_err) 405 mem_err_s = mem_err_s + mem_err 406 end if 407 enddo 408 deallocate(amg_I_colm%p,stat=mem_err) 409 mem_err_s = mem_err_s + mem_err 410 deallocate(amg_I_rowp%pp,stat=mem_err) 411 mem_err_s = mem_err_s + mem_err 412 deallocate(amg_I%pp,stat=mem_err) 413 mem_err_s = mem_err_s + mem_err 414 if (mem_err_s .ne. 0) then 415 write(*,*)'deallocation error in I_cf' 416 end if 417 418 ! now build I_fc, inverse of I_cf 419 call cpu_time(cpusec(6)) 420 421 mem_err_s = 0 422 allocate(I_fc_colm(level1)%p(amg_nshg(level2)+1),stat=mem_err) 423 mem_err_s = mem_err_s + mem_err 424 allocate(I_fc_rowp(level1)%p(mnnz),stat=mem_err) 425 mem_err_s = mem_err_s + mem_err 426 allocate(I_fc(level1)%p(mnnz),stat=mem_err) 427 mem_err_s = mem_err_s + mem_err 428 if (mem_err_s .ne. 0 ) then 429 write(*,*)'allocation error in I_fc' 430 end if 431 432 I_fc_colm(level1)%p(1:amg_nshg(level2)+1) = 0 433 do i=1,mnnz 434 I_fc_colm(level1)%p(I_cf_rowp(level1)%p(i)) = 435 & I_fc_colm(level1)%p(I_cf_rowp(level1)%p(i)) + 1 436 enddo 437 mnnz = 1 438 do i=1,amg_nshg(level2) 439 j = I_fc_colm(level1)%p(i) 440 I_fc_colm(level1)%p(i) = mnnz 441 mnnz = mnnz + j 442 enddo 443 I_fc_colm(level1)%p(amg_nshg(level2)+1) = mnnz 444 445 do i=1,amg_nshg(level1) 446 do k=I_cf_colm(level1)%p(i),I_cf_colm(level1)%p(i+1)-1 447 j = I_cf_rowp(level1)%p(k) 448 kj = I_fc_colm(level1)%p(j) 449 I_fc_colm(level1)%p(j) = I_fc_colm(level1)%p(j) + 1 450 I_fc_rowp(level1)%p(kj) = i 451 I_fc(level1)%p(kj) = I_cf(level1)%p(k) 452 enddo 453 enddo 454 455 do i=amg_nshg(level2),2,-1 456 I_fc_colm(level1)%p(i) = I_fc_colm(level1)%p(i-1) 457 enddo 458 I_fc_colm(level1)%p(1) = 1 459 460 call cpu_time(cpusec(10)) 461 462 end subroutine!ramg_coarse_setup 463 464!************************************************************** 465! Heap routines 466!************************************************************** 467 subroutine ramg_initheap(heap,invMap,wght,nheaps,ilen) 468 implicit none 469 integer :: nheaps,ilen 470 integer,dimension(0:ilen-1),intent(inout) :: heap,invMap 471 real(kind=8),dimension(0:ilen-1),intent(in) :: wght 472 473 integer :: i,j,k,t 474 475 do i=0,nheaps-1 476 heap(i) = i 477 enddo 478 479 do i=1,nheaps-1 480 k = i 481 j = ishft(k-1,-1) 482 do while ( ( k.gt.0 ).and.(wght(heap(k)).gt.wght(heap(j))) ) 483 t = heap(j) 484 heap(j) = heap(k) 485 heap(k) = t 486 k = j 487 j = ishft(j-1,-1) 488 enddo 489 enddo 490 491 do i=0,nheaps-1 492 invmap(heap(i)) = i 493 enddo 494 495 end subroutine ! ramg_initheap 496 497 subroutine ramg_popheap(heap,invmap,wght,nheaps,popid,ilen) 498 implicit none 499 integer,intent(inout) :: nheaps 500 integer,intent(in) :: popid,ilen 501 real(kind=8),dimension(0:ilen-1),intent(in) :: wght 502 integer,dimension(0:ilen-1),intent(inout) :: heap,invmap 503 504 integer i,j,k 505 nheaps = nheaps-1 506 heap(popid) = heap(nheaps) 507 invmap(heap(popid)) = popid 508 call ramg_adjheap(heap,invmap,wght,nheaps,popid,ilen) 509 510 end subroutine !ramg_popheap 511 512 subroutine ramg_adjheap(heap,invmap,wght,nheaps,popid,ilen) 513 implicit none 514 integer,intent(in) :: nheaps,ilen 515 integer,dimension(0:ilen-1),intent(inout) :: heap,invmap 516 real(kind=8),dimension(0:ilen-1),intent(in) :: wght 517 integer,intent(in) :: popid 518 519 integer i,j,k,t 520 521 i = popid 522 j = ishft(i-1,-1); 523 do while ( (i.gt.0).and.(wght(heap(j)).lt.wght(heap(i)))) 524 t = heap(i) 525 heap(i) = heap(j) 526 heap(j) = t 527 invmap(heap(i)) = i 528 invmap(heap(j)) = j 529 i = j 530 j = ishft(i-1,-1) 531 enddo 532 533 i = popid 534 do while (i.lt.(nheaps/2)) 535 j = 2*i+1 536 if ((j.lt.(nheaps-1)).and.(wght(heap(j)).lt.wght(heap(j+1)))) 537 & then 538 j = j + 1 539 end if 540 if (wght(heap(i)).gt.wght(heap(j))) then 541 exit 542 end if 543 t = heap(i) 544 heap(i) = heap(j) 545 heap(j) = t 546 invmap(heap(i)) = i 547 invmap(heap(j)) = j 548 i = j 549 enddo 550 551 end subroutine !ramg_adjheap 552 553!**************************************************** 554! ramg_CFsplit: split coarse/fine nodes 555! Input: matrix in sparse form 556! parameters: (filter array)(aggressive/standard) 557! Output: C/F splitting result, Strong corelation set amg_S 558! within given filter subset 559!**************************************************** 560 561 subroutine ramg_CFsplit(amg_CF,amg_S,anshg,annz,acolm,arowp, 562 & alhs,afmap,ilwork,BC,iBC,iper, 563 & eps_str,afilter,spflag) 564 use ramg_data 565 566 include "common.h" 567 include "mpif.h" 568 include "auxmpi.h" 569 570 ! parameters 571 integer,intent(in) :: anshg,annz 572 integer,intent(inout),dimension(anshg) :: amg_CF 573 integer,intent(in),dimension(anshg+1) :: acolm 574 integer,intent(in),dimension(annz) :: arowp 575 integer,intent(inout),dimension(annz) :: amg_S 576 real(kind=8),intent(in),dimension(annz) :: alhs 577 integer,intent(in),dimension(anshg) :: afmap 578 integer,intent(in),dimension(nlwork) :: ilwork 579 integer,intent(in),dimension(nshg) :: iBC,iper 580 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 581 582 real(kind=8) :: eps_str 583 integer,intent(in) :: afilter,spflag 584 585 ! afilter: indicating which subset of afmap for coarsening 586 ! -1 for self 587 ! spflag: aggressive or not 588 589 integer :: i,j,k,n,m 590 integer,dimension(anshg) :: aheap,ainvheap 591 real(kind=8) :: Lmax 592 real(kind=8),dimension(anshg) :: amg_L 593 real(kind=8),dimension(anshg) :: rowmax,prowmax 594 real(kind=8) :: MSCALE,LSCALE,deltal 595 integer :: flagS,nheaps,cfilter 596 ! setup S, S^T matrix in amg_S 597 ! after this, S(i,j) = 1,3 : S 598 ! S(i,j) = 2,3 : S^T 599 600 integer,allocatable,dimension(:) :: oneck 601 602 !amg_S = 0 603 604 allocate(oneck(numpe+1)) 605 oneck = 0 606 607 LSCALE = 0.005 608 MSCALE = 0.1*LSCALE/anshg 609 610 do i=1,anshg 611 cfilter = afmap(i) ! each row only deal with its own group 612 amg_CF(i) = 0 613 amg_L(i) = 0 614 rowmax(i) = 1E+10 615 prowmax(i) = -1e+10 616 do j=acolm(i)+1,acolm(i+1)-1 617 k = arowp(j) 618 if ((cfilter.eq.afmap(k)).and.(alhs(j).lt.rowmax(i))) then 619 rowmax(i) = alhs(j) 620 end if 621 if ((cfilter.eq.afmap(k)).and.(alhs(j).gt.prowmax(i))) then 622 prowmax(i) = alhs(j) 623 end if 624 enddo 625 rowmax(i) = eps_str*rowmax(i) 626 prowmax(i) = eps_str*prowmax(i) 627 ! if (anshg.lt.1000) then 628 !write(*,"((I5)(3E14.5))")i,rowmax(i),alhs(acolm(i)),prowmax(i) 629 ! endif 630 enddo 631 do i=1,anshg 632 cfilter = afmap(i) 633 p = iabs(cfilter) 634 oneck(p) = oneck(p)+1 635 do j=acolm(i)+1,acolm(i+1)-1 636 k = arowp(j) 637 if (cfilter.eq.afmap(k)) then ! same subset 638 if (alhs(j).lt.0) then ! negative entries 639 if (alhs(j) .lt. rowmax(i)) then 640 amg_S(j) = amg_S(j) + 1 641 end if 642 if (alhs(j) .lt. rowmax(k)) then 643 amg_S(j) = amg_S(j) + 2 644 end if 645 else ! positive entries 646 IF (.FALSE.) THEN 647 if (alhs(j) .gt. prowmax(i)) then 648 amg_S(j) = amg_S(j) + 1 649 end if 650 if (alhs(j) .gt. prowmax(k)) then 651 amg_S(j) = amg_S(j) + 2 652 end if 653 ENDIF 654 endif ! if negative or positive 655 endif 656 enddo 657 enddo 658 ! mark isolated unknowns "coarse", S^T=0, 659 ! keep territory which should not be coarsened "coarse" 660 ! keep slave coarse temporarily 661 do i=1,anshg 662 cfilter = afmap(i) 663 flagS = 0 664 do j=acolm(i),acolm(i+1)-1 665 k = arowp(j) 666 if ((amg_S(j).ne.0).and.(afmap(k).eq.cfilter)) then 667 flagS = 1 668 exit 669 end if 670 enddo 671 if ( flagS.eq.0) then 672 p = iabs(cfilter) 673 oneck(p) = oneck(p)-1 674! if (oneck(p).gt.0) then 675 if (.false.) then 676 amg_CF(i) = 2 ! mark fine ! NO, should mark isolated coarse 677 else 678 amg_CF(i) = 1 ! last coarse ! everybody coarse 679 endif 680 end if 681 enddo 682 683 k = 0 684 do i=1,anshg 685 if (amg_CF(i).eq.2) then 686 k = k+1 687 endif 688 enddo 689 if (k.eq.anshg) then 690 write(*,*)'mcheck coarsening error here' 691 stop 692 endif 693 ! setup initial lamda value 694 Lmax = -1.0e+6 ! Lmax stores max value of lambda, m points to it 695 do i=1,anshg 696 cfilter = afmap(i) 697 if (amg_CF(i).eq.0) then 698 ! unassigned with in the filter 699 amg_L(i) = i*MSCALE+iabs(cfilter) 700 do j=acolm(i),acolm(i+1)-1 701 k = arowp(j) 702 if ((afmap(k).eq.cfilter).and.(amg_S(j).ge.2)) then ! in S^T 703 deltal = 0 704 if (amg_CF(k).eq.0) then ! undecided 705 deltal = LSCALE 706 else if (amg_CF(k).eq.2) then ! free 707 deltal = 2*LSCALE 708 end if!filter 709 if (alhs(j).gt.0) deltal=-deltal 710 amg_L(i) = amg_L(i) + deltal 711 712 !if (alhs(j).gt.0) then 713 ! write(*,*)'hooh: ',alhs(j) 714 !endif 715 ! this is to seperate nodes in each section 716 ! in the heap. e.g. 0-1000 for neighbour with proc 0, 717 ! 2000-3000 for neighbour with proc 1... 718 end if 719 enddo 720 if (amg_L(i).gt.Lmax) then 721 Lmax = amg_L(i) 722 m = i 723 end if 724 end if 725 enddo 726 ! aLoc: heap, amg_Fn: invheap 727 nheaps = anshg 728 call ramg_initheap(aheap,ainvheap,amg_L,nheaps,anshg) 729 !write(*,*)oneck 730 ! setup coarse/fine grid 731 do while (nheaps.gt.0) 732 m = aheap(1)+1 733 if (amg_CF(m).ne.0) then 734 amg_L(m) = 0 735 call ramg_popheap(aheap,ainvheap,amg_L,nheaps, 736 & ainvheap(m),anshg) 737 cycle 738 endif 739 Lmax = amg_L(m)-iabs(afmap(m))-m*MSCALE ! reuse Lmax 740 if (Lmax.lt.LSCALE) then 741 k = iabs(afmap(m)) 742 if (oneck(k).eq.1) then 743 amg_CF(m) = 1 744 else 745 amg_CF(m) = 2 746 endif 747 oneck(iabs(afmap(m)))=oneck(iabs(afmap(m)))-1 748 749 amg_L(m) = 0!afmap(m)+m*MSCALE 750 call ramg_popheap(aheap,ainvheap,amg_L,nheaps, 751 & ainvheap(m),anshg) 752 cycle 753 end if 754 ! set coarse grid 755 amg_CF(m) = 1 ! mark coarse 756 amg_L(m) = 0!afmap(m)+m*MSCALE!0 !floor(amg_L(m)) 757 call ramg_popheap(aheap,ainvheap,amg_L,nheaps,ainvheap(m), 758 & anshg) 759 ! set fine grid 760 do j=acolm(m),acolm(m+1)-1 761 k = arowp(j) 762 if (afmap(k).eq.afmap(m)) then 763 if (amg_S(j).ge.2) then ! in m's S^T if !B 764 if (amg_CF(k).eq.0) then !not defined yet !if A 765 ! update flag 766 amg_CF(k) = 2 ! mark fine 767 amg_L(k) = 0!afmap(k)+k*MSCALE!0 !floor(amg_L(k)) 768 call ramg_popheap(aheap,ainvheap,amg_L,nheaps,ainvheap(k), 769 & anshg) 770 ! update lambda in k's S 771 do i = acolm(k),acolm(k+1)-1 772 n = arowp(i) 773 if (((amg_S(i).eq.1).or.(amg_S(i).eq.3)).and. 774 & (amg_CF(n).eq.0)) then 775 ! in k's S , U==>F 776 if (alhs(i).le.0) then ! negative 777 amg_L(n) = amg_L(n) + LSCALE 778 else ! positive 779 amg_L(n) = amg_L(n) - LSCALE 780 endif 781 call ramg_adjheap(aheap,ainvheap,amg_L,nheaps,ainvheap(n), 782 & anshg) 783 end if 784 enddo 785 end if ! if A 786 end if ! if B 787 if (mod(amg_S(j),2).eq.1) then !if B2 788 if (amg_CF(k).eq.0) then 789 if (alhs(j).le.0) then 790 amg_L(k) = amg_L(k) - LSCALE 791 else 792 amg_L(k) = amg_L(k) + LSCALE 793 endif 794 call ramg_adjheap(aheap,ainvheap,amg_L,nheaps,ainvheap(k), 795 & anshg) 796 end if 797 end if ! if B2 798 end if ! subset 799 enddo 800 enddo 801 k = 0 802 do i=1,anshg 803 if (amg_CF(i).eq.2) then 804 k = k+1 805 endif 806 enddo 807 if (k.eq.anshg) then 808 write(*,*)'mcheck coarsening error here,2' 809 stop 810 endif 811 !write(*,*)oneck 812 813 ! any U left 814 j = 0 815 do i = 1,anshg 816 if (amg_CF(i).eq.0) then 817 amg_CF(i) = 1 ! Mark every node coarse if isolated 818! amg_CF(i) = 2 ! Mark every node fine if isolated 819 end if 820 if (amg_CF(i).eq.1) then 821 j = j+1 822 endif 823 enddo 824 !write(*,*)'mcheck rank2b, j=',myrank,j 825 if (.false.) then 826 k = 0 827 do i=1,anshg 828 if (amg_CF(i).eq.1) then 829 k = k+1 830 do j=acolm(i)+1,acolm(i+1)-1 831 m = arowp(j) 832 if (amg_CF(m).eq.1) then ! two coarse connected 833 if (amg_S(j).ge.1) then 834 ! write(*,*)'mcheck err',alhs(j) 835 endif 836 endif 837 enddo 838 endif 839 enddo 840 if (k.eq.anshg) then 841 write(*,*)'mcheck coarsening error here,3,' 842 stop 843 endif 844 endif 845 !write(*,*)'mcheck end of cfsplit' 846 !call ramg_dump_i(amg_CF,anshg,1,'splitcf ') 847 deallocate(oneck) 848 849 end subroutine !ramg_CFsplit 850 851 subroutine ramg_update_cfmap(amg_CF,slevel) 852 use ramg_data 853 integer,intent(in) :: slevel 854 integer,intent(in),dimension(amg_nshg(slevel)) :: amg_CF 855 integer :: i,j,p 856 integer,dimension(amg_nshg(slevel)) :: revmap 857 858 revmap = 0 859 j = 1 860 do i=1,amg_nshg(1) 861 if (amg_cfmap(i).eq.slevel) then 862 revmap(j) = i 863 j = j+1 864 end if 865 enddo 866 do i=1,j-1 867 if (amg_CF(i).eq.1) then 868 amg_cfmap(revmap(i)) = amg_cfmap(revmap(i)) + 1 869 end if 870 enddo 871 872 end subroutine ! ramg_update_cfmap 873 874 subroutine ramg_readin_cfmap(amg_CF,slevel) 875 use ramg_data 876 integer,intent(in) :: slevel 877 integer,intent(inout),dimension(amg_nshg(slevel)) :: amg_CF 878 integer :: i,j,p 879 880 j = 1 881 do i=1,amg_nshg(1) 882 if (amg_cfmap(i).eq.slevel) then 883 amg_CF(j) = 2 ! fine 884 j = j + 1 885 else if (amg_cfmap(i).eq.(slevel+1)) then 886 amg_CF(j) = 1 ! coarse 887 j = j + 1 888 endif 889 enddo 890 891 end subroutine ! ramg_readin_cfmap 892