1 subroutine ramg_ggb_setup(colm,rowp,lhsK,lhsP, 2 & ilwork,BC,iBC,iper) 3 4 use ramg_data 5 include "common.h" 6 include "mpif.h" 7 include "auxmpi.h" 8 include 'debug.h' 9 10 integer,intent(in),dimension(nshg+1) :: colm 11 integer,intent(in),dimension(nnz_tot) :: rowp 12 13 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 14 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 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 iparam(11),ipntr(14) 20 21 integer gcomm 22 23 logical selectncv(maxncv) 24 real(kind=8) :: d(maxncv,3),dr(maxncv),di(maxncv), 25 & workev(3*maxncv), 26 & workl(3*maxncv*maxncv+6*maxncv) 27 28 real(kind=8),allocatable,dimension(:) :: resid,workd 29 30 character bmat*1,which*2 31 integer :: ido,n,nx,nev,ncv,lworkl,info,ierr, 32 & i,j,ishifts,maxitr,model,nconv 33 real(kind=8) :: tol,sigmar,sigmai 34 logical :: first,rvec 35 36 real(kind=8) :: rtemp,tnorm,tnorm1 37 real(kind=8),allocatable,dimension(:) :: tv,tw 38 real(kind=8),allocatable,dimension(:,:):: teA,tramg_ev 39 40 integer,allocatable,dimension(:) :: gmap,grevmap 41 42 integer :: p,q,nsize,step,asize,gsize 43 character fname*20 44 45 if (maxnev.le.0) return 46 if (ramg_setup_flag.ne.0) return 47 48 if (numpe.gt.1) then 49 call MPI_Barrier(MPI_COMM_WORLD,ierr) 50 end if 51 52 asize = amg_nshg(1) 53 allocate(gmap(asize)) 54 allocate(grevmap(asize)) 55 56 call ramg_generate_gmap(ilwork,asize,nsize,gmap,grevmap,1) 57 58 call MPI_AllReduce(nsize,gsize,1,MPI_INTEGER,MPI_MAX, 59 & MPI_COMM_WORLD,ierr) 60 61 !write(*,*)'mcheck ggb:',myrank,asize,nsize,gsize 62 63 allocate(resid(gsize)) 64 allocate(workd(3*gsize)) 65 !*********** Start of Parameter setting ******* 66 67 ndigit = -3 68 logfil = 6 69 mnaitr = 0 70 mnapps = 0 71 mnaupd = 1 ! controls the output (verbose) mode 72 mnaup2 = 0 73 mneigh = 0 74 mneupd = 1 75 76 nev = maxnev 77 ncv = maxncv 78 bmat = 'I' 79 which = 'LM' 80 81 lworkl = 3*ncv**2+6*ncv 82 tol = 1.0E-3 83 ido = 0 84 info = 0 85 86 iparam(1) = 1 87 iparam(3) = 500 88 iparam(7) = 1 89 90 !********** End of parameter setting ******** 91 92 allocate(tramg_ev(nsize,maxncv)) 93 allocate(tv(asize)) 94 allocate(tw(asize)) 95 96 step = 1 97 if (myrank.eq.master) then 98 write(*,*)'AMG GGB: calculating eigenvalue/vectors' 99 endif 100 tramg_ev = 0 101 gcomm = MPI_COMM_WORLD 102 103 call pdnaupd(gcomm, 104 & ido,bmat,nsize,which,nev,tol,resid,ncv, 105 & tramg_ev, 106 & nsize,iparam,ipntr,workd,workl,lworkl, 107 & info) 108 do while ((ido.eq.1).or.(ido.eq.-1)) 109 call ramg_ggb_av(workd(ipntr(2)),workd(ipntr(1)), 110 & asize,nsize,gmap,grevmap, 111 & 1, 112 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 113 call pdnaupd(gcomm, 114 & ido,bmat,nsize,which,nev,tol,resid,ncv, 115 & tramg_ev, 116 & nsize,iparam,ipntr,workd,workl,lworkl, 117 & info) 118 step = step + 1 119 enddo 120 !write(*,*)'mcheck: ggb over pdnaupd' 121 122 if (info.lt.0) then 123 write(*,*)'mcheck: ggb: info:',info 124 else 125 rvec = .false. 126 call pdneupd (gcomm,rvec,'A',selectncv,dr,di, 127 & tramg_ev,nsize, 128 & sigmar,sigmai,workev,bmat,nsize,which,nev,tol, 129 & resid,ncv,tramg_ev, 130 & nsize,iparam,ipntr,workd,workl, 131 & lworkl,ierr) 132 !write(*,*)'mcheck:ggb:',ierr,'iconv:',iparam(5) 133 end if 134 135 allocate(ramg_ggb_ev(asize,maxnev)) 136 do i=1,maxnev 137 call ramg_ggb_G2A(ramg_ggb_ev(:,i),tramg_ev(:,i), 138 & asize,nsize,1,gmap,grevmap) 139 enddo 140 141 142 ! Set matrix eA for GGB V-cycle 143 144 allocate(ramg_ggb_eA(maxnev,maxnev)) 145 allocate(teA(asize,maxnev)) 146 147 ramg_ggb_eA = 0 148 teA = 0 149 150 do i=1,maxnev 151 call ramg_calcAv_g(1,teA(:,i),ramg_ggb_ev(:,i), 152 & colm,rowp,lhsK,lhsP, 153 & ilwork,BC,iBC,iper,1) 154 enddo 155 do i=1,maxnev 156 do j=1,maxnev 157 do p=1,asize 158 ramg_ggb_eA(j,i) = ramg_ggb_eA(j,i) + 159 & ramg_ggb_ev(p,j)*teA(p,i) 160 enddo 161 enddo 162 enddo 163 164 ! communicate to complete GGB matrix 165 if (numpe.gt.1) then 166 call MPI_Barrier(MPI_COMM_WORLD,ierr) 167 do i=1,maxnev 168 do j=1,maxnev 169 call MPI_Allreduce(ramg_ggb_eA(j,i),rtemp,1, 170 & MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr) 171 ramg_ggb_eA(j,i) = rtemp 172 enddo 173 enddo 174 end if 175 176 !call MPI_Barrier(MPI_COMM_WORLD,ierr) 177 178 deallocate(teA) 179 180 deallocate(tramg_ev) 181 deallocate(tv) 182 deallocate(tw) 183 deallocate(gmap) 184 deallocate(grevmap) 185 deallocate(resid) 186 deallocate(workd) 187 188 ramg_time(11:30) = 0 189 !write(*,*)'mcheck: over at g-cycle' 190 191 end subroutine !ramg_ggb_setup 192 193 subroutine ramg_generate_gmap(ilwork,asize,nsize,gmap,grevmap, 194 & level) 195 196 use ramg_data 197 include "common.h" 198 include "mpif.h" 199 include "auxmpi.h" 200 201 integer,intent(inout) :: nsize 202 integer,intent(in),dimension(nlwork) :: ilwork 203 integer,intent(in) :: asize 204 integer,intent(inout),dimension(asize) :: gmap,grevmap 205 integer,intent(in) :: level 206 207 integer :: numtask,i,j,k,m,p,newn 208 integer :: itag,iacc,iother,numseg,isgbeg,itkbeg 209 integer,allocatable,dimension(:) :: neimap 210 211 do i=1,asize 212 gmap(i) = i 213 grevmap(i) = i 214 enddo 215 216 if (numpe.le.1) then 217 nsize = asize 218 return 219 end if 220 221 numtask = ilwork(1) 222 m = 0 223 itkbeg = 1 224 allocate(neimap(numpe)) 225 neimap = 0 226 227 newn = 0 228 229 do i=1,numtask 230 m = m+1 231 itag = ilwork(itkbeg+1) 232 iacc = ilwork(itkbeg+2) 233 iother = ilwork(itkbeg+3) 234 numseg = ilwork(itkbeg+4) 235 isgbeg = ilwork(itkbeg+5) 236 if (iacc.eq.0) then !slave 237 neimap(iother+1) = 1 238! do j=1,numseg 239! k = ilwork(itkbeg+3+2*j) 240! gmap(k) = -iother 241! enddo 242 end if 243 m = m + 1 244 itkbeg = itkbeg + 4 + 2*numseg 245 enddo 246 247 k = 1 248 do i=1,asize 249 m = iabs(amg_paramap(level)%p(i)) 250 if (neimap(m).eq.0) then ! self or master 251 gmap(i) = k 252 grevmap(k) = i 253 k = k+1 254 end if 255 enddo 256 nsize = k-1 257 deallocate(neimap) 258 259 !write(*,*)'mcheck gmap,',myrank,asize,nsize 260 261 !call ramg_dump_i(gmap,asize,1,'gmapdump ') 262 !call ramg_dump_i(grevmap,asize,1,'grevmap ') 263 264 end subroutine !ramg_generate_gmap 265 266 subroutine ramg_ggb_av( 267 & gramg_sol,gramg_rhs, 268 & asize,nsize,gmap,grevmap, 269 & clevel, 270 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper 271 & ) 272 use ramg_data 273 include "common.h" 274 include "mpif.h" 275 276 !arguments 277 integer, intent(in) :: asize,nsize 278 integer, intent(in) :: clevel 279 real(kind=8),intent(inout),dimension(nsize) 280 & :: gramg_sol 281 real(kind=8),intent(in),dimension(nsize) :: gramg_rhs 282 integer, intent(in),dimension(asize) :: gmap,grevmap 283 284 integer,intent(in),dimension(nshg+1) :: colm 285 integer,intent(in),dimension(nnz_tot) :: rowp 286 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 287 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 288 integer,intent(in),dimension(nlwork) :: ilwork 289 integer,intent(in),dimension(nshg) :: iBC,iper 290 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 291 292 !local 293 real(kind=8),dimension(amg_nshg(clevel)) :: vF 294 real(kind=8),dimension(asize) :: ramg_sol,ramg_rhs 295 296 call ramg_ggb_G2A(ramg_rhs,gramg_rhs,asize,nsize,1,gmap,grevmap) 297 298 vF = 0 299 300 call ramg_calcAv_g(clevel,vF,ramg_rhs,colm,rowp,lhsK,lhsP, 301 & ilwork,BC,iBC,iper,1) 302 303 call ramg_V_cycle(ramg_sol,vF,1, 304 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper) 305 ramg_sol = ramg_rhs-ramg_sol 306 307 call ramg_ggb_A2G(ramg_sol,gramg_sol,asize,nsize,1, 308 & gmap,grevmap) 309 310 end subroutine !ramg_ggb_av 311 312 313 subroutine ramg_G_cycle( 314 & ramg_sol,ramg_rhs, 315 & clevel, 316 & colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper 317 & ) 318 use ramg_data 319 include "common.h" 320 include "mpif.h" 321 322 !arguments 323 integer, intent(in) :: clevel 324 real(kind=8),intent(inout),dimension(amg_nshg(clevel)) 325 & :: ramg_sol,ramg_rhs 326 327 integer,intent(in),dimension(nshg+1) :: colm 328 integer,intent(in),dimension(nnz_tot) :: rowp 329 real(kind=8),intent(in),dimension(9,nnz_tot) :: lhsK 330 real(kind=8),intent(in),dimension(4,nnz_tot) :: lhsP 331 integer,intent(in),dimension(nlwork) :: ilwork 332 integer,intent(in),dimension(nshg) :: iBC,iper 333 real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC 334 335 !local 336 real(kind=8),dimension(maxnev) :: vR,rtemp 337 real(kind=8),dimension(amg_nshg(clevel)) :: ggberr 338 real(kind=8),dimension(maxnev,maxnev) :: mtxggbA 339 integer,dimension(maxnev) :: indx 340 real(kind=8) :: d 341 integer :: i,j 342 343 if (maxnev.le.0) return 344 345! ggberr = 0 346 call ramg_calcAv_g(clevel,ggberr,ramg_sol,colm,rowp,lhsK,lhsP, 347 & ilwork,BC,iBC,iper,1) 348 ggberr = ggberr - ramg_rhs 349 350 ! restriction 351 vR = 0 352 do i=1,maxnev 353 do j=1,amg_nshg(clevel) 354 vR(i) = vR(i)+ramg_ggb_eV(j,i)*ggberr(j) 355 enddo 356 enddo 357 rtemp = 0 358 if (numpe.gt.1) then 359 call MPI_Allreduce(vR,rtemp,maxnev, 360 & MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr) 361 362 vR = rtemp 363 end if 364 365 ! solve maxnev by maxnev matrix 366 mtxggbA = ramg_ggb_eA 367 call ludcmp(mtxggbA,maxnev,maxnev,indx,d) 368 call lubksb(mtxggbA,maxnev,maxnev,indx,vR) 369 370 ! prolongation 371 ggberr = 0 372 do i=1,maxnev 373 do j=1,amg_nshg(clevel) 374 ggberr(j) = ggberr(j)+ramg_ggb_eV(j,i)*vR(i) 375 enddo 376 enddo 377 378 ramg_sol = ramg_sol - ggberr 379 380 end subroutine !ramg_G_cycle 381 382 subroutine ramg_ggb_G2A(avec,gvec,asize,gsize,redun, 383 & gmap,grevmap) 384 integer,intent(in) :: asize, gsize 385 integer,intent(in) :: redun 386 integer,intent(in),dimension(asize) :: gmap,grevmap 387 real(kind=8),intent(inout),dimension(asize,redun) :: avec 388 real(kind=8),intent(in),dimension(gsize,redun) :: gvec 389 integer :: i 390 391 avec = 0 392 do i=1,gsize 393 avec(grevmap(i),:) = gvec(i,:) 394 enddo 395 396 end subroutine ! ramg_ggb_G2A 397 398 399 subroutine ramg_ggb_A2G(avec,gvec,asize,gsize,redun, 400 & gmap,grevmap) 401 integer,intent(in) :: asize, gsize 402 integer,intent(in) :: redun 403 integer,intent(in),dimension(asize) :: gmap,grevmap 404 real(kind=8),intent(in),dimension(asize,redun) :: avec 405 real(kind=8),intent(inout),dimension(gsize,redun) :: gvec 406 integer :: i 407 408 gvec = 0 409 do i=1,asize 410 if (gmap(i).gt.0) then 411 gvec(gmap(i),:) = avec(i,:) 412 end if 413 enddo 414 415 end subroutine ! ramg_ggb_A2G 416 417