1!************************************************************ 2! ramg_calcPPEv 3! calculate u = PPE*v ( parallelly using lesLib calls) 4!************************************************************ 5 subroutine ramg_calcPPEv(colm,rowp,lhsK,lhsP, 6 & ilwork,BC,iBC,iper, 7 & u,v) 8 use ramg_data 9 include "common.h" 10 include "mpif.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 integer,intent(in),dimension(nlwork) :: ilwork 17 integer,intent(in),dimension(nshg) :: iBC,iper 18 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 19 20 real(kind=8),intent(in),dimension(nshg) :: v 21 real(kind=8),intent(inout),dimension(nshg) :: u 22 23 real(kind=8),dimension(nshg,3) :: utmp 24 real(kind=8),dimension(nshg,4) :: utmp4 25 integer :: i,j 26 27 call commOut(v,ilwork,1,iper,iBC,BC) 28 call fLesSparseApG(colm,rowp,lhsP,v,utmp,nshg,nnz_tot) 29 call commIn(utmp,ilwork,3,iper,iBC,BC) 30 31 do i=1,3 32 utmp(:,i) = utmp(:,i)*ramg_flowDiag%p(:,i)**2 33 enddo 34 35 do i=1,3 36 utmp4(:,i) = -utmp(:,i) 37 enddo 38 39 utmp4(:,4) = v 40 41 call commOut(utmp4,ilwork,4,iper,iBC,BC) 42 call fLesSparseApNGtC(colm,rowp,lhsP,utmp4,u,nshg,nnz_tot) 43 call commIn(u,ilwork,1,iper,iBC,BC) 44! There is a slight modify at lesSparse.f: 45! row(20*nNodes) ==> row(nnz_tot) 46 47 end subroutine ! ramg_calcPPEv 48 49 50!************************************************************ 51! ramg_jacobi 52! calculate u = D^-1 * [ -( L+U ) v + r ] 53! 54! and also the residule in and out resout = | r-Au | 55! resin = | r-Av | 56! 57!************************************************************ 58 subroutine ramg_jacobi(Acolm,Arowp,Alhs,Anshg,Annz_tot,!A 59 & r,u,v) 60 61 include "common.h" 62 integer,intent(in) :: Anshg,Annz_tot 63 integer,intent(in),dimension(Anshg+1) :: Acolm 64 integer,intent(in),dimension(Annz_tot) :: Arowp 65 real(kind=8),intent(in),dimension(Annz_tot) :: Alhs 66 real(kind=8),intent(in),dimension(Anshg) :: v,r 67 real(kind=8),intent(inout),dimension(Anshg) :: u 68 real(kind=8) :: tmp,tmpin,tmpout 69 70 integer :: i,j,k,p 71 72 real(kind=8) :: damp_jacobi 73 74 ! Useless function, will be removed 75 76 u = 0 77 damp_jacobi = 1.0/ramg_relax 78 do i=1,Anshg 79 do k = Acolm(i)+1,Acolm(i+1)-1 80 j = Arowp(k) 81 tmp = Alhs(k)*v(j) 82 u(i) = u(i) - tmp 83 enddo 84 u(i) = damp_jacobi*u(i) + r(i) 85 u(i) = u(i)/Alhs(Acolm(i)) 86 enddo 87 88 end subroutine ! ramg_jacobi 89 90!*************************************************************** 91! ramg_boundary_mls1: 92! Do polynomial smoothing on boundary nodes only. 93! remember this can always be simplified to boundary jacobi 94!************************************************************** 95 subroutine ramg_boundary_mls1(acolm,arowp,alhs, 96 & anshg,annz_tot,!A 97 & r,u,v,clevel,pflag, 98 & ilwork,BC,iBC,iper) 99 use ramg_data 100 include "common.h" 101 include "mpif.h" 102 include "auxmpi.h" 103 104 integer,intent(in) :: anshg,annz_tot 105 integer,intent(in),dimension(anshg+1) :: acolm 106 integer,intent(in),dimension(annz_tot) :: arowp 107 real(kind=8),intent(in),dimension(annz_tot) :: alhs 108 real(kind=8),intent(in),dimension(anshg) :: v,r 109 real(kind=8),intent(inout),dimension(anshg) :: u 110 integer,intent(in) :: clevel 111 integer,intent(in),dimension(nlwork) :: ilwork 112 integer,intent(in),dimension(nshg) :: iBC,iper 113 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 114 logical,dimension(anshg),intent(inout) :: pflag 115 116 integer i,j,k,p,q 117 real(kind=8),dimension(anshg) :: tmp,v2,aux 118 real(kind=8) :: cf1 119 real(kind=8) :: cpusec(2) 120 121 if (numpe.eq.1) return 122 123 call cpu_time(cpusec(1)) 124 125 tmp = 0 126 v2 = v 127 128 ! tmp = A*x 129 call ramg_calcAv_b(clevel,tmp,v2,ilwork,BC,iBC,iper,1,0) 130! call ramg_calcAv(acolm,arowp,alhs,anshg,annz_tot, 131! & tmp,v2,1) 132 ! tmp = A*x-b=r 133 tmp = tmp-r 134 ! aux = A*r ! not necessary 135 ! call ramg_calcAv_b(clevel,aux,tmp,ilwork,BC,iBC,iper,1,0) 136 137 cf1 = 1.1*mlsCf(clevel,1) 138 139 !v2 = mlsCf(clevel,1)*tmp!+mlsCf(clevel,2)*aux 140 !v2 = 1.1*v2 141 !v2 = tmp 142 143 do i=1,anshg 144 p = amg_paraext(clevel)%p(i) 145 q = amg_paramap(clevel)%p(i) 146 if (p.ne.(myrank+1)) then 147 pflag(i) = .true. 148 if ((p.gt.0).or.(p.ne.q)) then ! master boundary or extended 149 u(i) = v(i)-cf1*tmp(i) 150 endif 151 endif 152 enddo 153 154 call cpu_time(cpusec(2)) 155 ramg_time(20) = ramg_time(20)+cpusec(2)-cpusec(1) 156 157 end subroutine !ramg_boundary_mls1 158 159 160!************************************************************ 161! ramg_gauss 162! 163! forward and backward, defined in fwdbck 164! Gauss-Seidel smoothing, u = gaussseidel(M,r) 165! 166!************************************************************ 167 subroutine ramg_gauss(acolm,arowp,alhs,anshg,annz_tot,!A 168 & r,u,v,fwdbck,clevel, 169 & ilwork,BC,iBC,iper) 170 use ramg_data 171 include "common.h" 172 include "mpif.h" 173 include "auxmpi.h" 174 175 integer,intent(in) :: anshg,annz_tot 176 integer,intent(in),dimension(anshg+1) :: acolm 177 integer,intent(in),dimension(annz_tot) :: arowp 178 real(kind=8),intent(in),dimension(annz_tot) :: alhs 179 real(kind=8),intent(in),dimension(anshg) :: v,r 180 real(kind=8),intent(inout),dimension(anshg) :: u 181 integer,intent(in) :: fwdbck!1=fwd,2=bck 182 integer,intent(in) :: clevel 183 integer,intent(in),dimension(nlwork) :: ilwork 184 integer,intent(in),dimension(nshg) :: iBC,iper 185 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 186 187 188 real(kind=8) :: tmp 189 real(kind=8),dimension(anshg) :: u2 190 integer :: istr,iend,iint 191 integer :: i,j,k,p,ki,kj,kp 192 logical :: postsmooth,p2 193 logical,dimension(anshg) :: pflag 194! integer :: itag,iacc,iother,numseg,isgbeg,itkbeg 195! integer :: numtask 196 197 u = v 198 199 postsmooth = (fwdbck.eq.1) 200 if (postsmooth) then ! post smoothing : 1->nshg 201 istr = 1 202 iend = anshg 203 iint = 1 204 else ! pre smoothing : nshg->1 205 istr = anshg 206 iend = 1 207 iint = -1 208 end if 209 210 pflag = .false. 211 ! boundary 212! if (postsmooth) then 213 call ramg_boundary_mls1(acolm,arowp,alhs,anshg,annz_tot, 214 & r,u,v,clevel,pflag,ilwork,BC,iBC,iper) 215! endif 216! if (.false.) then 217 218 do i=istr,iend,iint 219 ki = CF_map(clevel)%p(i) 220 p2 = .not.(pflag(ki)) ! nodes that have not been treated 221! if (numpe.gt.1) then 222! p2=p2.and.(amg_paraext(clevel)%p(ki).eq.(myrank+1)) 223! endif 224 if (p2) then 225 tmp = 0 226 do k = acolm(ki)+1,acolm(ki+1)-1 227 j = arowp(k) 228 if ((postsmooth).or.(pflag(j))) then 229 ! save some if u(j)=0 230 tmp = tmp + alhs(k)*u(j) 231 endif 232 enddo 233 u(ki) = r(ki) - tmp 234 pflag(ki) = .true. 235 !u(ki) = u(ki)/alhs(acolm(ki)) 236 ! The above line is commented out because we know that 237 ! A(i,i)=1.0, this is the result of prescaling of the matrix. 238 endif !p2 239 enddo 240! if (.not.postsmooth) then 241 u2 = u 242 call ramg_boundary_mls1(acolm,arowp,alhs,anshg,annz_tot, 243 & r,u,u2,clevel,pflag,ilwork,BC,iBC,iper) 244! endif 245 246 end subroutine ! ramg_gauss 247 248!********************************************** 249! ramg dense direct solver routines. 250! 251! ramg_sparse2dense 252! Sparse to dense form 253! 254! ramg_direct_LU ( setup only ) 255! Outside package for direct solve sparse 256! matrix by LU 257!********************************************** 258 259 subroutine ramg_direct_LU(Acolm,Arowp,Alhs,Anshg,Annz_tot) 260 use ramg_data 261 262 integer,intent(in) :: Anshg,Annz_tot 263 integer,intent(in),dimension(Anshg+1) :: Acolm 264 integer,intent(in),dimension(Annz_tot) :: Arowp 265 real(kind=8),intent(in),dimension(Annz_tot) :: Alhs 266 267! real(kind=8),dimension(Anshg,Anshg) :: mtxA 268! integer,dimension(Anshg) :: indx 269 real(kind=8) :: d 270 271 integer :: i,j 272 273 !write(*,*) "ramg_setup_flag",ramg_setup_flag 274 if (ramg_setup_flag .ne. 0) return 275 276 !write(*,*)'again' 277 if (allocated(cmtxA)) deallocate(cmtxA) 278 if (allocated(cindx)) deallocate(cindx) 279 280 allocate(cmtxA(Anshg,Anshg)) 281 allocate(cindx(Anshg)) 282 283 if (myrank.eq.master) then 284 write(*,*)"Start to setup LU decomposition at coarsest level" 285 endif 286 call ramg_sparse2dense(Acolm,Arowp,Alhs,Anshg,Annz_tot,cmtxA) 287 288! Asol = Arhs 289 290 call ludcmp(cmtxA,Anshg,Anshg,cindx,d) 291! call lubksb(mtxA,Anshg,Anshg,indx,Asol) 292 if (myrank.eq.master) then 293 write(*,*)"End of setup LU decomposition at coarsest level" 294 endif 295 296 end subroutine ! ramg_direct_LU 297 298 subroutine ramg_sparse2dense(Acolm,Arowp,Alhs,Anshg,Annz_tot, 299 & mtxA) 300 301 integer,intent(in) :: Anshg,Annz_tot 302 integer,intent(in),dimension(Anshg+1) :: Acolm 303 integer,intent(in),dimension(Annz_tot) :: Arowp 304 real(kind=8),intent(in),dimension(Annz_tot) :: Alhs 305 real(kind=8),intent(inout),dimension(Anshg,Anshg) :: mtxA 306 307 integer :: i,j,k,ip,jp,kp 308 309 mtxA = 0 310 do i=1,Anshg 311 do j= Acolm(i),Acolm(i+1)-1 312 k = Arowp(j) 313 mtxA(i,k) = Alhs(j) 314 enddo 315 enddo 316 317 end subroutine !ramg_sparse2dense 318 319 320!********************************************** 321! ramg_direct_solve 322!********************************************** 323 subroutine ramg_direct_solve( 324 & Arhs,Asol, 325 & colm,rowp,lhsK,lhsP, 326 & ilwork,BC,iBC,iper) 327 use ramg_data 328 include "common.h" 329 include "mpif.h" 330 include "auxmpi.h" 331 real(kind=8),intent(in),dimension(amg_nshg(ramg_levelx)) 332 & ::Arhs 333 real(kind=8),intent(inout),dimension(amg_nshg(ramg_levelx)) 334 & :: Asol 335 integer,intent(in),dimension(nshg+1) :: colm 336 integer,intent(in),dimension(nnz_tot) :: rowp 337 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 338 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 339 340 integer,intent(in),dimension(nlwork) :: ilwork 341 integer,intent(in),dimension(nshg) :: iBC,iper 342 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 343 344 real(kind=8) :: mres_i,mres_n,mres_o 345 real(kind=8),dimension(amg_nshg(ramg_levelx)) :: myV 346 integer :: titer,cl 347 348! Asol=Arhs 349! return 350 cl = ramg_levelx 351 352 Asol = 0 353 call ramg_L2_norm_p(cl,Arhs,1,mres_i) 354 !mres_i=sqrt(flesDot1(Arhs,amg_nshg(cl),1)) 355 mres_o = 1e-7*mres_i 356 titer = 0 357 358 !write(*,*)'mcheck direct solve begin,',myrank 359 if (iamg_c_solver.eq.0) then 360 ! Do two smoothing steps 361 call ramg_smoother(cl,myV,Asol,Arhs, 362 & colm,rowp,lhsK,lhsP, 363 & ilwork,BC,iBC,iper,2) 364 call ramg_smoother(cl,Asol,myV,Arhs, 365 & colm,rowp,lhsK,lhsP, 366 & ilwork,BC,iBC,iper,1) 367 else if ( iamg_c_solver .eq. 1) then ! direct solve using G-S 368 IF (.FALSE.) THEN ! smoother solve to exact 369 mres_n = mres_i 370 do while ( (mres_n.gt.mres_o).and.(titer.lt.1000) ) 371 call ramg_smoother(cl,myV,Asol,Arhs, 372 & colm,rowp,lhsK,lhsP, 373 & ilwork,BC,iBC,iper,2) 374 call ramg_smoother(cl,Asol,myV,Arhs, 375 & colm,rowp,lhsK,lhsP, 376 & ilwork,BC,iBC,iper,1) 377 call ramg_calcAv_g(cl,myV,Asol,colm,rowp,lhsK,lhsP, 378 & ilwork,BC,iBC,iper,1) 379 myV = myV-Arhs 380 mres_n=sqrt(flesDot1(myV,amg_nshg(cl),1)) 381 titer = titer + 1 382 enddo 383 ELSE ! cg solve to exact 384 call ramg_CG(Asol,Arhs,cl, 385 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper,titer) 386 ENDIF 387 else ! direct solve using dense-matrix LU-decomposition 388 Asol = Arhs 389 call lubksb(cmtxA,amg_nshg(cl),amg_nshg(cl),cindx,Asol) 390 endif 391 392 end subroutine ! ramg_direct_solve 393 394 395 subroutine ramg_smoother(level,xnew,xold,xrhs, 396 & colm,rowp,lhsK,lhsP, 397 & ilwork,BC,iBC,iper,fwdbck) 398 use ramg_data 399 include "common.h" 400 include "mpif.h" 401 include "auxmpi.h" 402 403 integer level 404 real(kind=8),dimension(amg_nshg(level)),intent(in) ::xrhs,xold 405 real(kind=8),dimension(amg_nshg(level)),intent(inout) :: xnew 406 integer,intent(in),dimension(nshg+1) :: colm 407 integer,intent(in),dimension(nnz_tot) :: rowp 408 409 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 410 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 411 412 integer,intent(in),dimension(nlwork) :: ilwork 413 integer,intent(in),dimension(nshg) :: iBC,iper 414 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 415 integer fwdbck 416 417 integer i,k 418 logical :: jacobi 419 420 ! NOTICE: fwdbck=2 : pre-smoothing 421 ! NOTICE: fwdbck=1 : post-smoothing 422 423 jacobi = .false. !(ramg_relax.gt.0) 424 ! Disable Jacobi, if needed, using 1-st order MLS instead 425 426 k = 1 427 do i=1,k 428 if ( iamg_smoother .eq. 3 ) then !.or. (level.gt.1) ) then 429 call ramg_mls_apply(xnew,xold,xrhs,level, 430 & colm,rowp,lhsK,lhsP, 431 & ilwork,BC,iBC,iper,fwdbck) 432 else if ( iamg_smoother .eq. 1 ) then 433 call ramg_gauss(amg_A_colm(level)%p,amg_A_rowp(level)%p, 434 & amg_A_lhs(level)%p,amg_nshg(level), 435 & amg_nnz(level),xrhs,xnew,xold, 436 & fwdbck,level, 437 & ilwork,BC,iBC,iper) 438 else if ( ( iamg_smoother .eq. 2) ) then !.and. (level.eq.1) ) then 439 call ramg_cheby_apply(xnew,xold,xrhs,level, 440 & colm,rowp,lhsK,lhsP, 441 & ilwork,BC,iBC,iper,fwdbck) 442 else if (jacobi) then 443 call ramg_jacobi(amg_A_colm(level)%p, 444 & amg_A_rowp(level)%p,amg_A_lhs(level)%p, 445 & amg_nshg(level),amg_nnz(level), 446 & xrhs,xnew,xold) 447 endif 448 enddo 449 450 end subroutine !ramg_smoother 451