1!*********************************************************** 2! ramg_cheby_apply: u = Cheby(v,(lhs,r)) 3! v : approx. solution before 4! u : approx. solution (after) 5! r : residual vector, r=Ax-b 6!*********************************************************** 7 8 9 subroutine ramg_cheby_apply(u,v,r,level,colm,rowp,lhsK,lhsP, 10 & ilwork,BC,iBC,iper,fwdbck) 11 use ramg_data 12 include "common.h" 13 include "mpif.h" 14 include "auxmpi.h" 15 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 integer,intent(in),dimension(nshg+1) :: colm 20 integer,intent(in),dimension(nnz_tot) :: rowp 21 22 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 23 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 24 25 26 real(kind=8),intent(in),dimension(amg_nshg(level)) :: v,r 27 real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u 28 integer,intent(in) :: fwdbck 29 30 ! Local 31 real(kind=8),dimension(amg_nshg(level)) :: pAux 32 real(kind=8),dimension(amg_nshg(level)) :: epdk,epx 33 34 real(kind=8),dimension(amg_nshg(level)) :: cy1,cz1 35 36 real(kind=8) :: cf,beta,alpha,delta,theta,s1,rhok,invtheta 37 real(kind=8) :: rhokp1,dtemp1,dtemp2,twoodelta 38 real(kind=8) :: chebyniu,chebyomega,chebyck1,chebyck2,chebygamma 39 integer :: i,j,k 40 41 beta = 1.1*mlsCf(level,1) ! 1.1 * Ev_max 42 alpha = mlsCf(level,1)*ramg_chebyratio 43 44 IF (.FALSE.) THEN 45 ! mine ( matrix computing ) 46 chebyniu = 1.0+2.0*(1.0-beta)/(beta-alpha) 47 chebygamma = 2.0/(2.0-alpha-beta) 48 chebyomega = chebyniu/(2.0*chebyniu*chebyniu-1.0) 49 chebyomega = chebyomega*2.0*(2.0-beta-alpha)/(beta-alpha) 50 if (fwdbck.eq.2) then ! initial guess x0 = 0 51 cy1 = r 52 else 53 call ramg_calcAv_g(level,pAux,v,colm,rowp,lhsK,lhsP, 54 & ilwork,BC,iBC,iper,1) 55 cy1 = v-1.1*3.0/beta*(pAux-r) 56 endif 57 58 call ramg_calcAv_g(level,pAux,cy1,colm,rowp,lhsK,lhsP, 59 & ilwork,BC,iBC,iper,1) 60 cz1 = r-pAux 61 62 if (fwdbck.eq.2) then 63 u = chebyomega*(cy1+chebygamma*cz1) 64 else 65 u = chebyomega*(cy1-v+chebygamma*cz1)+v 66 endif 67 68 ELSE 69 ! trilinos 70 delta = (beta-alpha)/2.0 71 twoodelta = 4.0/(beta-alpha) 72 theta = (beta+alpha)/2.0 73 invtheta = 2.0/(beta+alpha) 74 s1 = 2.0*theta/delta 75 rhok = delta/theta 76 77 !write(*,*)delta,twoodelta,2.0/delta 78 if (fwdbck.eq.2) then ! pre-smoothing, initial guess v=0 79 epdk = r*invtheta!/theta 80 epx = epdk 81 else 82 call ramg_calcAv_g(level,pAux,v,colm,rowp,lhsK,lhsP, 83 & ilwork,BC,iBC,iper,1) 84 epdk = r-pAux 85 epdk = epdk*invtheta!/theta 86 epx = v+epdk 87 endif 88 89 do k=1,mlsDeg-1 90 call ramg_calcAv_g(level,pAux,epx,colm,rowp,lhsK,lhsP, 91 & ilwork,BC,iBC,iper,1) 92 rhokp1 = 1.0/(s1-rhok) 93 dtemp1 = rhokp1*rhok 94 dtemp2 = rhokp1*twoodelta 95 rhok = rhokp1 96 epdk = dtemp1*epdk+dtemp2*(r-pAux) 97 epx = epx + epdk 98 enddo 99 u = epx 100 101 ENDIF 102 103 end subroutine ! ramg_cheby_apply 104 105 106 subroutine ramg_cheby_setup(colm,rowp,lhsK,lhsP, 107 & ilwork,BC,iBC,iper) 108 use ramg_data 109 include "common.h" 110 include "mpif.h" 111 include "auxmpi.h" 112 113 integer,intent(in),dimension(nshg+1) :: colm 114 integer,intent(in),dimension(nnz_tot) :: rowp 115 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 116 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 117 integer,intent(in),dimension(nlwork) :: ilwork 118 integer,intent(in),dimension(nshg) :: iBC,iper 119 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 120 121 integer :: i,j,k,lvl,ideg 122 123 if (ramg_setup_flag .ne. 0 ) return 124 125 if (myrank.eq.master) then 126 write(*,*)'Start Chebyshev setup' 127 endif 128 do lvl=1,ramg_levelx 129 call ramg_mls_eigen(mlsCf(lvl,1),lvl,1,colm,rowp,lhsK,lhsP, 130 & ilwork,BC,iBC,iper) 131! call ramg_mls_min_eigen(mlsCf(lvl,2),lvl,1,colm,rowp,lhsK,lhsP, 132! & ilwork,BC,iBC,iper) 133! mlsCf(lvl,2) = mlsCf(lvl,1) 134 if (myrank.eq.master) then 135 write(*,7000)lvl,mlsCf(lvl,1) 1367000 format('Level: ',I2,T15,'Ev_max: ',E10.4) 137 endif 138 enddo 139 if (myrank.eq.master) then 140 write(*,*)'End Chebyshev setup' 141 endif 142 143 end subroutine !ramg_cheby_setup 144