1!************************************ 2! RAMG Extract: 3! extract ppe and scale ppe 4!********************************************************** 5 6!*********************************************************** 7! ramg_extract_ppe 8! prepare PPE matrix 9! input: boundary, lhsK,lhsP 10! output: PPE & RHS at level 1 11!********************************************************** 12 subroutine ramg_extract_ppe( 13 & colm,rowp,lhsK,lhsP, 14 & ilwork,BC,iBC,iper) 15 use ramg_data 16 include "common.h" 17 include "mpif.h" 18 include "auxmpi.h" 19! implicit none 20 21 22 !***********parameters************** 23 !the matrix 24! integer,intent(in) :: nshg 25! integer,intent(in) :: nnz_tot 26! integer,intent(in) :: nflow 27 !the matrix 28 integer,intent(in),dimension(nshg+1) :: colm 29 integer,intent(in),dimension(nnz_tot) :: rowp 30 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 31 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 32! real(kind=8),dimension(:,:),allocatable :: lhsGP 33! integer,dimension(:),allocatable :: lhsGProwp 34! integer,dimension(:),allocatable :: lhsGPcolm 35 type(r2d) :: lhsGP 36 type(i1d) :: lhsGProwp,lhsGPcolm 37 ! the boundary info 38 integer, intent(in), dimension(nlwork) :: ilwork 39 integer, intent(in),dimension(nshg) :: iBC,iper 40 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 41 42 !*********parameters end************** 43 44 !****** local variables ********** 45 real(kind=8),dimension(nshg,4) :: mflowDiag,pflowDiag 46 real(kind=8) :: rtemp,rt,rtp,rtn 47 real(kind=8),dimension(nshg) :: rtest 48 integer :: i,j,k,m,n,p,q,ki,kj,ni,nj,qq 49 integer :: cki,ckj,ckk 50 integer,dimension(nshg) :: temprow 51 integer :: mnnz 52 53 logical :: extentri 54 55 character :: fname*80 56 !****** end local variables ****** 57 58 if (ramg_setup_flag .ne. 0) return 59 60 !*** calculate memory for PPE *** 61 !*** using the fact that matrix is symmetric *** 62 !*** lhs = GT_ik * KI_kk * G_kj + C_ij *** 63 !*** rhs = -GT_ik * KI_kk * Rm_k - Rc_k *** 64 !*** where G and GT has same sparsity pattern if 65 !*** only consider the main matrix (M,3) 66 67 if ( amg_nshg(1) .gt.0 ) then 68 call ramg_deallocate(1) 69 deallocate(ramg_flowDiag%p) 70 deallocate(amg_cfmap) 71 end if 72 73 ! colm 74 call ramg_allocate(1,nshg,0,1) 75 allocate(amg_paramap(1)%p(nshg)) 76 allocate(amg_paraext(1)%p(nshg)) 77 78 mflowDiag(:,:)=0 79 80! if ( (numpe.gt.1) ) then !.and.(iamg_reduce.gt.0) ) then 81! call ramg_prep_reduce_map 82 !call ramg_dump_i(ncorp_map,nshg,1,'map_nproc ') 83 !deallocate(ncorp_map) 84! endif 85 86 call drvLesPrepDiag(mflowDiag,ilwork,iBC,BC,iper, 87 & rowp,colm,lhsK,lhsP) 88 89 !*** Complete diagonal values in mflowdiag,pflowdiag *** 90 !*** lhs K, G, C should have "global" values on diagonal *** 91 !call ramg_dump_rn_map(mflowDiag,nshg,4,'mflowDiagb') 92 call commIn(mflowdiag,ilwork,4,iper,iBC,BC) 93 call commOut(mflowdiag,ilwork,4,iper,iBC,BC) 94 !call ramg_dump_rn_map(mflowDiag,nshg,4,'mflowDiaga') 95 ! call ramg_dump(mflowDiag,nshg,'flowdiag ') 96 97 allocate(ramg_flowDiag%p(nshg,4)) 98 allocate(amg_cfmap(nshg)) 99 100 amg_cfmap = 1 101 ! amg_cfmap is such a variable that it keeps 102 ! coarsening information through the coarsest level 103 ! in a finest level array to enable communication. 104 ! The structure is : 111223322113231122... 105 ! obviously the nodes that kept into next level got 106 ! an extra 1 107 108 do i=1,4 109 ramg_flowDiag%p(:,i) = mflowDiag(:,i) 110 do j=1,nshg 111 if (mflowdiag(j,4).eq.0) then 112 write(*,*)'mflowdiag zero at ',j,i 113 stop 114 end if 115 enddo 116 enddo 117 118 call ramg_initBCflag(amg_paramap(1)%p,ilwork,BC,iBC,iper) 119 120 call ramg_global_lhs(colm,rowp,lhsP,nnz_tot, 121 & lhsGPcolm,lhsGProwp,lhsGP,4, 122 & ilwork,BC,iBC,iper) 123 124 ! prepare extended boundary in amg_paraext 125 ni = 0 126 nj = 0 127 amg_paraext(1)%p = amg_paramap(1)%p 128 do i=1,nshg 129 if (amg_paramap(1)%p(i).ne.(myrank+1)) then 130 ni = ni+1 131 nj = nj+1 132 do m=lhsGPcolm%p(i),lhsGPcolm%p(i+1)-1 133 k = lhsGProwp%p(m) 134 if (amg_paraext(1)%p(k).eq.(myrank+1)) then 135 nj = nj+1 136 amg_paraext(1)%p(k)=amg_paramap(1)%p(i) 137 endif 138 enddo 139 endif 140 enddo 141 !write(*,*)'proc',myrank,ni,nj,nshg 142 143 ! end of preparing extended boundary 144 145 extentri = .false. ! HAVE extra entries 146 !extentri = .true. ! NO extra entries 147 148 ! output K^{hat}^-1 and lhsP (G,C) to external file 149 !write(fname,*)'Khatinv' 150 !call ramg_dump_rn(mflowdiag,nshg,4,fname) 151 !write(fname,*)'lhsP' 152 !call ramg_dump_matlab_A(colm,rowp,lhsGP,nshg,nnz_tot,4,fname) 153 154 mnnz = 1 155 temprow = 0 156 do i = 1,nshg ! each row 157 amg_A_colm(1)%p(i) = mnnz 158 !q = iabs(amg_paramap(1)%p(i)) ! used only in reduced serial 159 do m = lhsGPcolm%p(i),lhsGPcolm%p(i+1)-1 160 k = lhsGProwp%p(m) 161 do j = lhsGPcolm%p(k),lhsGPcolm%p(k+1)-1 162 p = lhsGProwp%p(j) 163 !qq = iabs(amg_paramap(1)%p(p)) ! used only in reduced 164 !serial 165 if (temprow(p).ne.i) then 166! if ((iamg_reduce.gt.1).and.(numpe.eq.1).and.(extentri))then 167! if ( (rmap1d(1,qq).eq.rmap1d(1,q)).or. 168! & (rmap1d(2,qq).eq.rmap1d(2,q)).or. 169! & (rmap1d(1,qq).eq.rmap1d(2,q)).or. 170! & (rmap1d(2,qq).eq.rmap1d(1,q)) ) then 171! temprow(p) = i 172! mnnz = mnnz + 1 173! endif 174! else 175 temprow(p) = i 176 mnnz = mnnz + 1 177! endif 178 end if 179 end do 180 end do 181 end do 182 amg_A_colm(1)%p(nshg+1) = mnnz 183 mnnz = mnnz - 1 184 !*** now we have amg_nnz(1) as nnz_tot for PPE 185 call ramg_allocate(1,0,mnnz,1) 186 !allocate(levelbaseflag(mnnz)) 187 !levelbaseflag = 0 188 189 if (iamg_verb.gt.5) then 190 write(6,7003) nshg,amg_nnz(1) 1917003 format(/'nshg='i10', nnz='i10) 192 endif 193 194 !*** end of PPE memory alloc *** 195 196 !**** begin putting rowp and values into PPE ** 197 mnnz = 0 198 temprow = 0 199 200 201 ! rowp 202 do i = 1,nshg 203 mnnz = mnnz + 1 204 amg_A_rowp(1)%p(mnnz) = i 205 temprow(i) = i 206 !levelbaseflag(mnnz) = 1 207 !q = iabs(amg_paramap(1)%p(i)) !used in reduced serial 208 do m = lhsGPcolm%p(i),lhsGPcolm%p(i+1)-1 209 j = lhsGProwp%p(m) 210 do n = lhsGPcolm%p(j),lhsGPcolm%p(j+1)-1 211 k = lhsGProwp%p(n) 212 !qq = iabs(amg_paramap(1)%p(k)) !reduced serial 213 if ( temprow(k) .ne. i) then 214! if ((iamg_reduce.gt.1).and.(numpe.eq.1).and.(.true.))then 215! if ( (rmap1d(1,qq).eq.rmap1d(1,q)).or. 216! & (rmap1d(2,qq).eq.rmap1d(2,q)).or. 217! & (rmap1d(1,qq).eq.rmap1d(2,q)).or. 218! & (rmap1d(2,qq).eq.rmap1d(1,q)) ) then 219! mnnz = mnnz + 1 220! amg_A_rowp(1)%p(mnnz) = k 221! temprow(k) = i 222! levelbaseflag(mnnz+1) = 1 223! endif 224! endif 225! else 226 mnnz = mnnz + 1 227 amg_A_rowp(1)%p(mnnz) = k 228 temprow(k) = i 229! endif 230 end if 231 enddo 232 enddo 233 do m=amg_A_colm(1)%p(i)+1,amg_A_colm(1)%p(i+1)-1 234 ki = amg_A_rowp(1)%p(m) 235 do n=m+1,amg_A_colm(1)%p(i+1)-1 236 kj = amg_A_rowp(1)%p(n) 237 if (kj.lt.ki) then 238 amg_A_rowp(1)%p(m) = kj 239 amg_A_rowp(1)%p(n) = ki 240 ki = kj 241! p = levelbaseflag(n) 242! levelbaseflag(n) = levelbaseflag(m) 243! levelbaseflag(m) = p 244 end if 245 enddo 246 enddo 247 enddo 248 249! j = 0 250! do i=1,mnnz 251! j = j + levelbaseflag(i) 252! enddo 253! write(*,*)'mcheck incomplete: ',j 254 255 rt = 0 256 257 ! matrix value 258 ! cki,ckj,ckk, this is to avoid double summation on PPE, 259 ! since both master and slave have complete lhsP, parts of it 260 ! should be removed when constructing PPE. if (i,j) are both from 261 ! slave 262 do i=1,nshg 263 do m = amg_A_colm(1)%p(i),amg_A_colm(1)%p(i+1)-1 264 j = amg_A_rowp(1)%p(m) 265 !if (levelbaseflag(m).eq.0) then 266 ! amg_A_lhs(1)%p(m,1) = 0 267 ! cycle 268 !endif 269 ki = lhsGPcolm%p(i) 270 ni = lhsGPcolm%p(i+1) 271 kj = lhsGPcolm%p(j) 272 nj = lhsGPcolm%p(j+1) 273 rtemp = 0 274 ! question: the original lhsK,lhsP, are they symmetric? 275 ! symmetric in nnz structure but not in value? 276 kloop: do while ( ( ki.lt.ni ).and.(kj.lt.nj) ) 277 do while ((ki.lt.ni).and. 278 & (lhsGProwp%p(ki).lt.lhsGProwp%p(kj)) ) 279 ki = ki + 1 280 enddo 281 if (ki.eq.ni) exit kloop 282 do while ( (kj.lt.nj) .and. 283 & (lhsGProwp%p(kj).lt.lhsGProwp%p(ki)) ) 284 kj = kj + 1 285 enddo 286 if (kj.eq.nj) exit kloop 287 p = lhsGProwp%p(ki) 288 q = lhsGProwp%p(kj) 289 if (p.eq.q) then 290c if (amg_paramap(1)%p(p).eq.amg_paramap(1)%p(q)) then 291 k = q 292 cki = amg_paramap(1)%p(i) 293 ckj = amg_paramap(1)%p(j) 294 ckk = amg_paramap(1)%p(k) 295 if (.not.((iabs(cki).eq.iabs(ckj)).and. 296 & (iabs(ckj).eq.iabs(ckk)).and. 297 & (cki.lt.0).and.(ckj.lt.0).and.(ckk.lt.0))) then 298 rtemp = rtemp 299 & + lhsGP%p(1,ki)*lhsGP%p(1,kj)*mflowDiag(k,1)**2 300 & + lhsGP%p(2,ki)*lhsGP%p(2,kj)*mflowDiag(k,2)**2 301 & + lhsGP%p(3,ki)*lhsGP%p(3,kj)*mflowDiag(k,3)**2 302 endif 303 ki = ki+1 304 kj = kj+1 305c endif 306 end if 307 enddo kloop 308 amg_A_lhs(1)%p(m,1)=rtemp 309 enddo 310 enddo 311 do i=1,nshg 312 ki = amg_A_colm(1)%p(i)+1 313 mloop: do m=lhsGPcolm%p(i),lhsGPcolm%p(i+1)-1 314 j = lhsGProwp%p(m) 315 if (j.eq.i) then 316 cki = amg_paramap(1)%p(i) 317 ckj = amg_paramap(1)%p(j) 318 if (.not.((iabs(cki).eq.iabs(ckj)).and. 319 & (cki.lt.0).and.(ckj.lt.0))) then 320 amg_A_lhs(1)%p(amg_A_colm(1)%p(i),1) = 321 & + amg_A_lhs(1)%p(amg_A_colm(1)%p(i),1) + lhsGP%p(4,m) 322 endif 323 cycle mloop 324 end if 325 do while(amg_A_rowp(1)%p(ki).lt.j) 326 ki = ki+1 327 enddo 328 cki = amg_paramap(1)%p(i) 329 ckj = amg_paramap(1)%p(j) 330 ckk = amg_paramap(1)%p(amg_A_rowp(1)%p(ki)) 331 if (.not.((iabs(cki).eq.iabs(ckj)).and. 332 & (iabs(ckj).eq.iabs(ckk)).and. 333 & (cki.lt.0).and.(ckj.lt.0).and.(ckk.lt.0))) then 334 amg_A_lhs(1)%p(ki,1) = amg_A_lhs(1)%p(ki,1) + lhsGP%p(4,m) 335 endif 336 ki = ki+1 337 enddo mloop 338 enddo 339 340 !gtg: lhsgp(4,m), mflowdiag ignored 341 342 deallocate(lhsGP%p) 343 deallocate(lhsGProwp%p) 344 deallocate(lhsGPcolm%p) 345 346 if (.true.) then 347 348 call ramg_global_lhs(amg_A_colm(1)%p,amg_A_rowp(1)%p, 349 & amg_A_lhs(1)%p,amg_nnz(1), 350 & lhsGPcolm,lhsGProwp,lhsGP,1, 351 & ilwork,BC,iBC,iper) 352 353 amg_A_colm(1)%p = lhsGPcolm%p 354 amg_nnz(1) = lhsGPcolm%p(amg_nshg(1)+1)-1 355 deallocate(amg_A_rowp(1)%p) 356 deallocate(amg_A_lhs(1)%p) 357 allocate(amg_A_rowp(1)%p(amg_nnz(1))) 358 allocate(amg_A_lhs(1)%p(amg_nnz(1),1)) 359 amg_A_rowp(1)%p = lhsGProwp%p 360 amg_A_lhs(1)%p = lhsGP%p 361 deallocate(lhsGP%p) 362 deallocate(lhsGProwp%p) 363 deallocate(lhsGPcolm%p) 364 365 endif 366 367 if (.false.) then 368 ! Dump for Matlab 369 call ramg_dump_matlab_map(amg_A_colm(1)%p,amg_A_rowp(1)%p, 370 & amg_A_lhs(1)%p, 371 & amg_nshg(1),amg_nnz(1),1, 372 & 'A0 ') 373! call ramg_dump_matlab_map(lhsGPcolm%p,lhsGProwp%p, 374! & lhsGP%p, 375! & amg_nshg(1),lhsGPcolm%p(amg_nshg(1)+1)-1,1, 376! & 'A0 ') 377 !call ramg_dump(amg_ppeDiag(1)%p,nshg,'ppediag ') 378 end if 379 380 381 if (.false.) then ! Check Diagonal Scaling, M-matrix 382 open(264,file='mtcheck.dat',status='unknown') 383 do i=1,amg_nshg(1) 384 m = 0 385 rtemp=amg_A_lhs(1)%p(amg_A_colm(1)%p(i),1) 386 rtp = 0 387 rtn = 0 388 do j=amg_A_colm(1)%p(i)+1,amg_A_colm(1)%p(i+1)-1 389 rtemp = rtemp + amg_A_lhs(1)%p(j,1) 390 rt = amg_A_lhs(1)%p(j,1) 391 if (rt.gt.0) then 392 if ( rt.gt.rtp ) rtp = rt 393 end if 394 if (rt.lt.0) then 395 if (rt.le.rtn) rtn = rt 396 end if 397 enddo 398 rtp = rtp/ (amg_A_lhs(1)%p(amg_A_colm(1)%p(i),1)) 399 rtn = rtn/amg_A_lhs(1)%p(amg_A_colm(1)%p(i),1) 400 write(264,"((I5)(A)(E9.3)(A)(E9.3)(A)(E9.3))") 401 & i,' ',rtp,' ',rtn,' ',rtemp 402 enddo 403 close(264) 404 end if 405 406 ! scaled by diag(PPE) 1...1 407 do i=1,nshg 408 amg_ppeDiag(1)%p(i) = 409 & 1.0/sqrt(amg_A_lhs(1)%p(amg_A_colm(1)%p(i),1)) 410 enddo 411 412 do i=1,nshg 413 do m=amg_A_colm(1)%p(i),amg_A_colm(1)%p(i+1)-1 414 j = amg_A_rowp(1)%p(m) 415 amg_A_lhs(1)%p(m,1) = amg_A_lhs(1)%p(m,1) * 416 & amg_ppeDiag(1)%p(i)*amg_ppeDiag(1)%p(j) 417 end do 418 amg_A_rhs(1)%p(i) = amg_A_rhs(1)%p(i)*amg_ppeDiag(1)%p(i) 419 end do 420 421 ! If solve heat/scalar, should disable following 3 lines 422 ! because we only have 1 scaling. 423 do i=1,nshg 424 amg_ppeDiag(1)%p(i) = amg_ppeDiag(1)%p(i)/mflowDiag(i,4) 425 enddo 426 427 if (.false.) then 428 ! Dump for Matlab 429 call ramg_dump_matlab_map(amg_A_colm(1)%p,amg_A_rowp(1)%p, 430 & amg_A_lhs(1)%p, 431 & amg_nshg(1),amg_nnz(1),1, 432 & 'A0scaled ') 433! call ramg_dump_matlab_map(lhsGPcolm%p,lhsGProwp%p, 434! & lhsGP%p, 435! & amg_nshg(1),lhsGPcolm%p(amg_nshg(1)+1)-1,1, 436! & 'A0 ') 437 !call ramg_dump(amg_ppeDiag(1)%p,nshg,'ppediag ') 438 end if 439 440 441 return 442 443 end subroutine !ramg_extract_ppe 444 445