1!************************************************************** 2! mls_eigen: get maximum eigenvalue of A 3! mls_setup: set coefficients of A^n 4! 5! These routine uses ARPACK and auxilary routines in 6! ramg_ggb.f 7!************************************************************** 8 subroutine ramg_mls_eigen(evmax,level,flagfb, 9 & colm,rowp,lhsK,lhsP, 10 & ilwork,BC,iBC,iper) 11 use ramg_data 12 include "common.h" 13 include "mpif.h" 14 include "auxmpi.h" 15 include 'debug.h' 16 17 real(kind=8),intent(inout) :: evmax 18 integer,intent(in) :: level 19 integer,intent(in) :: flagfb 20 integer,intent(in),dimension(nshg+1) :: colm 21 integer,intent(in),dimension(nnz_tot) :: rowp 22 23 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 24 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 25 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 30! ! parameter for PPE Ap-product 31! real(kind=8),dimension(nshg,4) :: diag 32 ! parameter for eigenvalue 33 34 integer iparam(11),ipntr(14) 35 integer gcomm 36 37! logical selectncv(maxncv) 38! real(kind=8) :: d(maxncv,3),dr(maxncv),di(maxncv), 39! & workev(3*maxncv), 40! & workl(3*maxncv*maxncv+6*maxncv) 41 logical,dimension(:),allocatable :: selectncv 42 real(kind=8),dimension(:,:),allocatable :: d 43 real(kind=8),dimension(:),allocatable :: dr,di,workev,workl 44 45 real(kind=8),allocatable,dimension(:) :: resid,workd 46 47 character bmat*1,which*2 48 integer :: ido,n,nx,nev,ncv,lworkl,info,ierr, 49 & i,j,ishifts,maxitr,model,nconv 50 real(kind=8) :: tol,sigmar,sigmai 51 logical :: first,rvec 52 53 real(kind=8) :: res_i,res_o,rtemp,tnorm,tnorm1 54 real(kind=8),allocatable,dimension(:,:):: tramg_ev 55 real(kind=8),dimension(amg_nshg(level)) :: twork1,twork2 56 57 integer,allocatable,dimension(:) :: gmap,grevmap 58 59 integer :: p,q,nsize,step,asize,gsize 60 61! ! if ( freeze parameter ) return 62! call drvLesPrepDiag(diag,ilwork,iBC,BC,iper,rowp,colm,lhsK,lhsP) 63 64 asize = amg_nshg(level) 65 allocate(gmap(asize)) 66 allocate(grevmap(asize)) 67 68 call ramg_generate_gmap(ilwork,asize,nsize,gmap,grevmap,level) 69 70 call MPI_AllReduce(nsize,gsize,1,MPI_INTEGER,MPI_MAX, 71 & MPI_COMM_WORLD,ierr) 72 73 ncv = min(maxncv,gsize) 74 !write(*,*)'mcheck ncv=',ncv 75 76 if (ncv.eq.1) then 77 evmax = 1 78 return 79 endif 80 81 allocate(selectncv(ncv)) 82 allocate(d(ncv,3)) 83 allocate(dr(ncv)) 84 allocate(di(ncv)) 85 allocate(workev(3*ncv)) 86 allocate(workl(3*ncv*ncv+6*ncv)) 87 88 allocate(resid(gsize)) 89 allocate(workd(3*gsize)) 90 !*********** Start of Parameter setting ******* 91 92 ndigit = -3 93 logfil = 6 94 mnaitr = 0 95 mnapps = 0 96 mnaupd = 0 ! controls the output (verbose) mode 97 mnaup2 = 0 98 mneigh = 0 99 mneupd = 0 100 101 nev = 1 102 bmat = 'I' 103 which = 'LM' 104 105 lworkl = 3*ncv**2+6*ncv 106 tol = 1.0E-3 107 ido = 0 108 info = 0 109 110 iparam(1) = 1 111 iparam(3) = 500 112 iparam(7) = 1 113 114 !********** End of parameter setting ******** 115 116 allocate(tramg_ev(nsize,maxncv)) 117 118 step = 1 119 tramg_ev = 0 120 gcomm = MPI_COMM_WORLD 121 122 call pdnaupd(gcomm, 123 & ido,bmat,nsize,which,nev,tol,resid,ncv, 124 & tramg_ev, 125 & nsize,iparam,ipntr,workd,workl,lworkl, 126 & info) 127 do while ((ido.eq.1).or.(ido.eq.-1)) 128 call ramg_ggb_G2A(twork1,workd(ipntr(1)),asize,nsize,1, 129 & gmap,grevmap) 130! call ramg_PPEAp(twork2,twork1,diag, 131! & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 132! call ramg_calcAv(amg_A_colm(level)%p,amg_A_rowp(level)%p, 133! & amg_A_lhs(level)%p,amg_nshg(level), 134! & amg_nnz(level),twork2,twork1,1) 135 if (flagfb.eq.1) then 136 call ramg_calcAv_g(level,twork2,twork1,colm,rowp,lhsK,lhsP, 137 & ilwork,BC,iBC,iper,1) 138 else 139 call ramg_mls_calcPAv(level,twork2,twork1, 140 & colm,rowp,lhsK,lhsP, 141 & ilwork,BC,iBC,iper) 142 endif 143 call ramg_ggb_A2G(twork2,workd(ipntr(2)),asize,nsize,1, 144 & gmap,grevmap) 145 call pdnaupd(gcomm, 146 & ido,bmat,nsize,which,nev,tol,resid,ncv, 147 & tramg_ev, 148 & nsize,iparam,ipntr,workd,workl,lworkl, 149 & info) 150 step = step + 1 151 enddo 152 if (info.lt.0) then 153 if (myrank.eq.master) then 154 write(*,*)'mcheck: mls: info:',info,iparam(5) 155 endif 156 ramg_levelx = level-1 157 else 158 rvec = .false. 159 call pdneupd (gcomm,rvec,'A',selectncv,dr,di, 160 & tramg_ev,nsize, 161 & sigmar,sigmai,workev,bmat,nsize,which,nev,tol, 162 & resid,ncv,tramg_ev, 163 & nsize,iparam,ipntr,workd,workl, 164 & lworkl,ierr) 165 end if 166 evmax = dr(1) 167! if (myrank.eq.master) then 168! write(*,*)'MLS: largest eigenvalue of A at level ',level, 169! & ' is ',evmax 170! endif 171 !write(*,*)'mcheck: ggb over pdnaupd' 172 173 deallocate(tramg_ev) 174 deallocate(gmap) 175 deallocate(grevmap) 176 deallocate(resid) 177 deallocate(workd) 178 179 deallocate(selectncv) 180 deallocate(d) 181 deallocate(dr) 182 deallocate(di) 183 deallocate(workev) 184 deallocate(workl) 185 186 end subroutine ! ramg_mls_eigen 187 188 subroutine ramg_mls_min_eigen(evmax,level,flagfb, 189 & colm,rowp,lhsK,lhsP, 190 & ilwork,BC,iBC,iper) 191 use ramg_data 192 include "common.h" 193 include "mpif.h" 194 include "auxmpi.h" 195 include 'debug.h' 196 197 real(kind=8),intent(inout) :: evmax 198 integer,intent(in) :: level 199 integer,intent(in) :: flagfb 200 integer,intent(in),dimension(nshg+1) :: colm 201 integer,intent(in),dimension(nnz_tot) :: rowp 202 203 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 204 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 205 206 integer,intent(in),dimension(nlwork) :: ilwork 207 integer,intent(in),dimension(nshg) :: iBC,iper 208 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 209 210! ! parameter for PPE Ap-product 211! real(kind=8),dimension(nshg,4) :: diag 212 ! parameter for eigenvalue 213 214 integer iparam(11),ipntr(14) 215 integer gcomm 216 217! logical selectncv(maxncv) 218! real(kind=8) :: d(maxncv,3),dr(maxncv),di(maxncv), 219! & workev(3*maxncv), 220! & workl(3*maxncv*maxncv+6*maxncv) 221 logical,dimension(:),allocatable :: selectncv 222 real(kind=8),dimension(:,:),allocatable :: d 223 real(kind=8),dimension(:),allocatable :: dr,di,workev,workl 224 225 real(kind=8),allocatable,dimension(:) :: resid,workd 226 227 character bmat*1,which*2 228 integer :: ido,n,nx,nev,ncv,lworkl,info,ierr, 229 & i,j,ishifts,maxitr,model,nconv 230 real(kind=8) :: tol,sigmar,sigmai 231 logical :: first,rvec 232 233 real(kind=8) :: res_i,res_o,rtemp,tnorm,tnorm1 234 real(kind=8),allocatable,dimension(:,:):: tramg_ev 235 real(kind=8),dimension(amg_nshg(level)) :: twork1,twork2 236 237 integer,allocatable,dimension(:) :: gmap,grevmap 238 239 integer :: p,q,nsize,step,asize,gsize 240 241! ! if ( freeze parameter ) return 242! call drvLesPrepDiag(diag,ilwork,iBC,BC,iper,rowp,colm,lhsK,lhsP) 243 244 asize = amg_nshg(level) 245 allocate(gmap(asize)) 246 allocate(grevmap(asize)) 247 248 call ramg_generate_gmap(ilwork,asize,nsize,gmap,grevmap,level) 249 250 call MPI_AllReduce(nsize,gsize,1,MPI_INTEGER,MPI_MAX, 251 & MPI_COMM_WORLD,ierr) 252 253 ncv = min(maxncv,gsize) 254 !write(*,*)'mcheck ncv=',ncv 255 256 if (ncv.eq.1) then 257 evmax = 1 258 return 259 endif 260 261 allocate(selectncv(ncv)) 262 allocate(d(ncv,3)) 263 allocate(dr(ncv)) 264 allocate(di(ncv)) 265 allocate(workev(3*ncv)) 266 allocate(workl(3*ncv*ncv+6*ncv)) 267 268 allocate(resid(gsize)) 269 allocate(workd(3*gsize)) 270 !*********** Start of Parameter setting ******* 271 272 ndigit = -3 273 logfil = 6 274 mnaitr = 0 275 mnapps = 0 276 mnaupd = 0 ! controls the output (verbose) mode 277 mnaup2 = 0 278 mneigh = 0 279 mneupd = 0 280 281 nev = 1 282 bmat = 'I' 283 which = 'SM' 284 285 lworkl = 3*ncv**2+6*ncv 286 tol = 1.0E-3 287 ido = 0 288 info = 0 289 290 iparam(1) = 1 291 iparam(3) = 500 292 iparam(7) = 1 293 294 !********** End of parameter setting ******** 295 296 allocate(tramg_ev(nsize,maxncv)) 297 298 step = 1 299 tramg_ev = 0 300 gcomm = MPI_COMM_WORLD 301 302 call pdnaupd(gcomm, 303 & ido,bmat,nsize,which,nev,tol,resid,ncv, 304 & tramg_ev, 305 & nsize,iparam,ipntr,workd,workl,lworkl, 306 & info) 307 do while ((ido.eq.1).or.(ido.eq.-1)) 308 call ramg_ggb_G2A(twork1,workd(ipntr(1)),asize,nsize,1, 309 & gmap,grevmap) 310! call ramg_PPEAp(twork2,twork1,diag, 311! & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 312! call ramg_calcAv(amg_A_colm(level)%p,amg_A_rowp(level)%p, 313! & amg_A_lhs(level)%p,amg_nshg(level), 314! & amg_nnz(level),twork2,twork1,1) 315 if (flagfb.eq.1) then 316 call ramg_calcAv_g(level,twork2,twork1,colm,rowp,lhsK,lhsP, 317 & ilwork,BC,iBC,iper,1) 318 else 319 call ramg_mls_calcPAv(level,twork2,twork1, 320 & colm,rowp,lhsK,lhsP, 321 & ilwork,BC,iBC,iper) 322 endif 323 call ramg_ggb_A2G(twork2,workd(ipntr(2)),asize,nsize,1, 324 & gmap,grevmap) 325 call pdnaupd(gcomm, 326 & ido,bmat,nsize,which,nev,tol,resid,ncv, 327 & tramg_ev, 328 & nsize,iparam,ipntr,workd,workl,lworkl, 329 & info) 330 step = step + 1 331 enddo 332 if (info.lt.0) then 333 if (myrank.eq.master) then 334 write(*,*)'mcheck: mls: info:',info,iparam(5) 335 endif 336 ramg_levelx = level-1 337 else 338 rvec = .false. 339 call pdneupd (gcomm,rvec,'A',selectncv,dr,di, 340 & tramg_ev,nsize, 341 & sigmar,sigmai,workev,bmat,nsize,which,nev,tol, 342 & resid,ncv,tramg_ev, 343 & nsize,iparam,ipntr,workd,workl, 344 & lworkl,ierr) 345 end if 346 evmax = dr(1) 347! if (myrank.eq.master) then 348! write(*,*)'MLS: smallest eigenvalue of A is ',evmax 349! endif 350 !write(*,*)'mcheck: ggb over pdnaupd' 351 352 deallocate(tramg_ev) 353 deallocate(gmap) 354 deallocate(grevmap) 355 deallocate(resid) 356 deallocate(workd) 357 358 deallocate(selectncv) 359 deallocate(d) 360 deallocate(dr) 361 deallocate(di) 362 deallocate(workev) 363 deallocate(workl) 364 365 end subroutine ! ramg_mls_min_eigen 366 367 subroutine ramg_mls_coeff(level,ilwork,BC,iBC,iper) 368 end subroutine !ramg_mls_coeff 369 370 371 subroutine ramg_mls_setup(colm,rowp,lhsK,lhsP, 372 & ilwork,BC,iBC,iper) 373 use ramg_data 374 include "common.h" 375 include "mpif.h" 376 include "auxmpi.h" 377 378 integer,intent(in),dimension(nshg+1) :: colm 379 integer,intent(in),dimension(nnz_tot) :: rowp 380 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 381 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 382 integer,intent(in),dimension(nlwork) :: ilwork 383 integer,intent(in),dimension(nshg) :: iBC,iper 384 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 385 386 387 real(kind=8) :: evmax,ddeg,aux1,aux0,auxom,rho,rho2 388 real(kind=8),dimension(10) :: omloc 389 real(kind=8) :: mlsOver ! magic number! 390 integer :: i,j,k,lvl,ideg 391 integer :: nSample,nGrid 392 real(kind=8) :: gridStep,coord,samplej 393 394 if (ramg_setup_flag .ne. 0 ) return 395 if ((iamg_smoother.eq.1).and.(numpe.eq.1)) return 396 ! do not setup if serial Gauss-Seidel 397 ! magic numbers.... 398 mlsOver = 1.1 399 400 do lvl=1,ramg_levelx 401 call ramg_mls_eigen(evmax,lvl,1,colm,rowp,lhsK,lhsP, 402 & ilwork,BC,iBC,iper) 403 404 mlsCf(lvl,:) = 0 405 omloc = 0 406 407 if (iamg_smoother.eq.1) then !Gauss-Seidel in parallel 408 ideg = 1 409 else 410 ideg = iabs(mlsDeg) 411 endif 412 ddeg = ideg 413 aux1 = 1.0/(2.0*ddeg+1.0) 414 rho = evmax!*mlsOver 415 416 do i=1,ideg 417 aux0 = 2.0*i*PI 418 auxom = rho/2.0*(1.0-cos(aux0*aux1)) 419 omloc(i) = 1.0/auxom 420 !write(*,*)'mcheck: ',omloc(i) 421 enddo 422 423 ! 4th order maximum 424 mlsCf(lvl,1) = omloc(1)+omloc(2)+omloc(3)+omloc(4)+omloc(5) 425 426 mlsCf(lvl,2) = -( omloc(1)*omloc(2)+omloc(1)*omloc(3) 427 & +omloc(1)*omloc(4)+omloc(1)*omloc(5) 428 & +omloc(2)*omloc(3)+omloc(2)*omloc(4) 429 & +omloc(2)*omloc(5)+omloc(3)*omloc(4) 430 & +omloc(3)*omloc(5)+omloc(4)*omloc(5) ) 431 432 mlsCf(lvl,3) = ( omloc(1)*omloc(2)*omloc(3) 433 & +omloc(1)*omloc(2)*omloc(4) 434 & +omloc(1)*omloc(2)*omloc(5) 435 & +omloc(1)*omloc(3)*omloc(4) 436 & +omloc(1)*omloc(3)*omloc(5) 437 & +omloc(1)*omloc(4)*omloc(5) 438 & +omloc(2)*omloc(3)*omloc(4) 439 & +omloc(2)*omloc(3)*omloc(5) 440 & +omloc(2)*omloc(4)*omloc(5) 441 & +omloc(3)*omloc(4)*omloc(5) ) 442 443 mlsCf(lvl,4) = -( omloc(1)*omloc(2)*omloc(3)*omloc(4) 444 & +omloc(1)*omloc(2)*omloc(3)*omloc(5) 445 & +omloc(1)*omloc(2)*omloc(4)*omloc(5) 446 & +omloc(1)*omloc(3)*omloc(4)*omloc(5) 447 & +omloc(2)*omloc(3)*omloc(4)*omloc(5) ) 448 449 mlsCf(lvl,5) = omloc(1)*omloc(2)*omloc(3)*omloc(4)*omloc(5) 450 451 do i=1,ideg 452 mlsOm(lvl,i) = omloc(i) 453 enddo 454 455 ! setup coefficients for post smoothing 456 ! ref: trilinos ml 457 458 nSample=20000 459 gridStep = rho/nSample 460 rho2=zero 461 do j=1,nSample 462 coord=j*gridStep 463 samplej=1.0-omloc(1)*coord 464 do i=2,ideg 465 samplej=samplej*(1.0-omloc(i)*coord) 466 enddo 467 samplej=samplej*(samplej*coord) 468 if (samplej.gt.rho2) rho2=samplej 469 enddo 470 smlsOm(lvl,1) = 2.0/(1.025*rho2)!*mlsCf(lvl,1) 471 !write(*,*)'mcheckpost:',rho2,smlsOm(lvl,1) 472 473 !do i=1,mlsDeg 474 ! write(*,*)'mcheck: ',lvl,i,mlsCf(lvl,i) 475 !enddo 476 if (myrank.eq.master) then 477 write(*,*)'MLS: Setup finished at level:',lvl,mlsCf(lvl,1) 478 endif 479 enddo !lvl 480 481 end subroutine ! ramg_mls_setup 482 483!*********************************************************** 484! Construct Sfc = Ifc*(I-\lambda_1 A-\lambda_2 A^2) 485! 486!*********************************************************** 487! Additional note: This is useless which results in 488! a super dense S instead of I, using this interpolator 489! is not effective won't save time ignore! 490 491!*********************************************************** 492! u = (S^2*A) v 493! S = polynomial smoothing 494! 495!********************************************************** 496 subroutine ramg_mls_calcPAv(level,u,v,colm,rowp,lhsK,lhsP, 497 & ilwork,BC,iBC,iper) 498 use ramg_data 499 include "common.h" 500 include "mpif.h" 501 include "auxmpi.h" 502 503 integer,intent(in),dimension(nlwork) :: ilwork 504 integer,intent(in),dimension(nshg) :: iBC,iper 505 integer,intent(in),dimension(nshg+1) :: colm 506 integer,intent(in),dimension(nnz_tot) :: rowp 507 508 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 509 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 510 511 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 512 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: v 513 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 514 515 real(kind=8),dimension(amg_nshg(level)) :: aux,aux2,r0 516 517! r0 = 0 518 call ramg_calcAv_g(level,aux,v,colm,rowp,lhsK,lhsP, 519 & ilwork,BC,iBC,iper,1) 520 call ramg_mls_sandw_pre(aux2,aux,level,colm,rowp,lhsK,lhsP, 521 & ilwork,BC,iBC,iper)!, 522 call ramg_mls_sandw_post(u,aux2,level,colm,rowp,lhsK,lhsP, 523 & ilwork,BC,iBC,iper)!, 524 525 end subroutine ! ramg_mls_calcPav 526 527 subroutine ramg_mls_setupPost(colm,rowp,lhsK,lhsP, 528 & ilwork,BC,iBC,iper) 529 use ramg_data 530 include "common.h" 531 include "mpif.h" 532 include "auxmpi.h" 533 534 integer,intent(in),dimension(nshg+1) :: colm 535 integer,intent(in),dimension(nnz_tot) :: rowp 536 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 537 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 538 integer,intent(in),dimension(nlwork) :: ilwork 539 integer,intent(in),dimension(nshg) :: iBC,iper 540 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 541 542 real(kind=8) :: evmax,ddeg,aux1,aux0,auxom,rho 543 real(kind=8),dimension(10) :: omloc 544 integer :: i,j,k,lvl 545 546 do lvl = 1,ramg_levelx 547 call ramg_mls_eigen(evmax,lvl,2,colm,rowp,lhsK,lhsP, 548 & ilwork,BC,iBC,iper) 549 smlsOm(lvl,2) = 1.0/(1.025*evmax)!*mlsCf(lvl,1) 550 !mlsOm(lvl,2) = -1.0*evmax*mlsCf(lvl,1)*mlsCf(lvl,1) 551 write(*,*)'mcheck: post,',evmax,smlsOm(lvl,2) 552 enddo 553 554 end subroutine ! ramg_mls_setupPost 555 556!*********************************************************** 557! ramg_mls_apply: u = MLS(v,(lhs,r)) 558! v : approx. solution before 559! u : approx. solution (after) 560! r : residual vector, r=Ax-b 561!*********************************************************** 562 563 subroutine ramg_mls_apply(u,v,r,level,colm,rowp,lhsK,lhsP, 564 & ilwork,BC,iBC,iper,fwdbck) 565 use ramg_data 566 include "common.h" 567 include "mpif.h" 568 include "auxmpi.h" 569 570 integer,intent(in),dimension(nlwork) :: ilwork 571 integer,intent(in),dimension(nshg) :: iBC,iper 572 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 573 integer,intent(in),dimension(nshg+1) :: colm 574 integer,intent(in),dimension(nnz_tot) :: rowp 575 576 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 577 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 578 579 580 real(kind=8),intent(in),dimension(amg_nshg(level)) :: v,r 581 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 582 integer,intent(in) :: fwdbck 583 584 real(kind=8),dimension(amg_nshg(level)) :: aux 585 586! calculate residual first if necessary 587! do not calculate residual at pre-smoothing, 588! only calculate it at post-smoothing. 589 if (fwdbck.eq.2) then 590 aux = -r 591 else 592 call ramg_calcAv_g(level,aux,v,colm,rowp,lhsK,lhsP, 593 & ilwork,BC,iBC,iper,1) 594 aux = aux-r 595 endif 596! if (.false.) then ! always post smoother 597 if (.true.) then ! always pre smoother 598 call ramg_mls_apply_fwd(u,v,aux,level,colm,rowp,lhsK,lhsP, 599 & ilwork,BC,iBC,iper) 600 else 601 call ramg_mls_apply_post(u,v,aux,level,colm,rowp,lhsK,lhsP, 602 & ilwork,BC,iBC,iper) 603 endif 604 605 u = v-1.1*u 606 607 end subroutine ! ramg_mls_apply 608 609 610 subroutine ramg_mls_apply_fwd(u,v,r,level,colm,rowp,lhsK,lhsP, 611 & ilwork,BC,iBC,iper) 612 use ramg_data 613 include "common.h" 614 include "mpif.h" 615 include "auxmpi.h" 616 617 integer,intent(in),dimension(nlwork) :: ilwork 618 integer,intent(in),dimension(nshg) :: iBC,iper 619 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 620 integer,intent(in),dimension(nshg+1) :: colm 621 integer,intent(in),dimension(nnz_tot) :: rowp 622 623 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 624 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 625 626 627 real(kind=8),intent(in),dimension(amg_nshg(level)) :: v,r 628 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 629 630 ! Local 631 real(kind=8),dimension(amg_nshg(level)) :: pAux,res 632 real(kind=8) :: cf 633 integer :: i,j,k 634 635 pAux = r 636 ! c_1*A*r 637 638 cf = mlsCf(level,1) 639 u = cf*pAux 640 do i=2,mlsDeg 641 ! c_i*A*(A*...*r) 642 call ramg_calcAv_g(level,res,pAux,colm,rowp,lhsK,lhsP, 643 & ilwork,BC,iBC,iper,1) 644 pAux = res 645 cf = mlsCf(level,i) 646 u = u+cf*res 647 enddo 648 649 end subroutine ! ramg_mls_apply_fwd 650 651 subroutine ramg_mls_apply_post(u,v,r,level,colm,rowp,lhsK,lhsP, 652 & ilwork,BC,iBC,iper) 653 use ramg_data 654 include "common.h" 655 include "mpif.h" 656 include "auxmpi.h" 657 658 integer,intent(in),dimension(nlwork) :: ilwork 659 integer,intent(in),dimension(nshg) :: iBC,iper 660 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 661 integer,intent(in),dimension(nshg+1) :: colm 662 integer,intent(in),dimension(nnz_tot) :: rowp 663 664 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 665 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 666 667 668 real(kind=8),intent(in),dimension(amg_nshg(level)) :: v,r 669 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 670 671 ! Local 672 real(kind=8),dimension(amg_nshg(level)) :: pAux,y,res,y1 673 real(kind=8) :: cf 674 integer :: i,j,k 675 676 call ramg_mls_apply_fwd(y1,v,r,level,colm,rowp,lhsK,lhsP, 677 & ilwork,BC,iBC,iper) 678! y1=1.1*y1 679 call ramg_calcAv_g(level,pAux,y1,colm,rowp,lhsK,lhsP, 680 & ilwork,BC,iBC,iper,1) 681 res = b-pAux!b-pAux ! rhs to feed in post smoothing, keep y 682 683 call ramg_mls_sandw_post(pAux,res,level,colm,rowp,lhsK,lhsP, 684 & ilwork,BC,iBC,iper) 685 call ramg_mls_sandw_pre(y,pAux,level,colm,rowp,lhsK,lhsP, 686 & ilwork,BC,iBC,iper) 687 u = smlsOm(level,1)*y+y1 688 689 end subroutine ! ramg_mls_apply_post 690 691 subroutine ramg_mls_sandw_pre(u,v,level,colm,rowp,lhsK,lhsP, 692 & ilwork,BC,iBC,iper) 693 use ramg_data 694 include "common.h" 695 include "mpif.h" 696 include "auxmpi.h" 697 698 integer,intent(in),dimension(nlwork) :: ilwork 699 integer,intent(in),dimension(nshg) :: iBC,iper 700 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 701 integer,intent(in),dimension(nshg+1) :: colm 702 integer,intent(in),dimension(nnz_tot) :: rowp 703 704 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 705 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 706 707 708 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: v 709 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 710 711 ! Local 712 real(kind=8),dimension(amg_nshg(level)) :: res 713 real(kind=8) :: cf 714 integer :: i,j,k 715 716 u = v 717 ! c_1*A*r 718 719 do i=mlsDeg,1,-1 720 call ramg_calcAv_g(level,res,u,colm,rowp,lhsK,lhsP, 721 & ilwork,BC,iBC,iper,1) 722 u = u-mlsOm(level,i)*res 723 enddo 724 725 end subroutine ! ramg_mls_sandw_pre 726 727 subroutine ramg_mls_sandw_post(u,v,level,colm,rowp,lhsK,lhsP, 728 & ilwork,BC,iBC,iper) 729 use ramg_data 730 include "common.h" 731 include "mpif.h" 732 include "auxmpi.h" 733 734 integer,intent(in),dimension(nlwork) :: ilwork 735 integer,intent(in),dimension(nshg) :: iBC,iper 736 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 737 integer,intent(in),dimension(nshg+1) :: colm 738 integer,intent(in),dimension(nnz_tot) :: rowp 739 740 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 741 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 742 743 744 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: v 745 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 746 747 ! Local 748 real(kind=8),dimension(amg_nshg(level)) :: res 749 real(kind=8) :: cf 750 integer :: i,j,k 751 752 u = v 753 ! c_1*A*r 754 755 do i=1,mlsDeg,1 756 call ramg_calcAv_g(level,res,u,colm,rowp,lhsK,lhsP, 757 & ilwork,BC,iBC,iper,1) 758 u = u-mlsOm(level,i)*res 759 enddo 760 761 end subroutine ! ramg_mls_sandw_post 762 763