c=========================================================== c Driver/testing routine for mgv ... c=========================================================== implicit none real*8 mgv integer stderr, stdin, stdout parameter ( stderr = 0, stdin = 5, stdout = 6 ) real*8 s2r8 real*8 r8arg integer s2i4 integer i4arg character*2 itoc character*128 sarg, catsqz real*8 r8_never parameter ( r8_never = -1.0d-60 ) integer i4_never parameter ( i4_never = -2 000 000 000 ) integer iargc, qalloc, qleft, qused real*8 mgvcyc, dvnrm2, dvlinf integer p, pceil common / comp / p, pceil c----------------------------------------------------------------------- c Working storage ... c c 1 Meg = 2**20 c----------------------------------------------------------------------- integer maxq parameter ( maxq = 16 * 2**20 ) real*8 q(maxq) integer maxrank parameter ( maxrank = 3 ) integer rank integer shape(maxrank), size real*8 vh(maxrank) integer uf, uxctf, rhsf, duf, * mresf, xf integer lf, nf, ncyc, icyc, * ncycle, preswf, pstswf, presw, * pstsw, p_save real*8 resnrm, e2nrm, einfnrm, hf logical ltrace parameter ( ltrace = .true. ) logical compute_residual parameter ( compute_residual = .true. ) if( iargc() .lt. 2 ) go to 100 rank = i4arg(1,-1) if( rank .ne. 1 .and. rank .ne. 2 .and. rank .ne. 3 ) then write(0,*) 'mg0: Only rank-1, 2 and 3 implemented' stop end if lf = i4arg(2,-1) if( lf .lt. 1 ) go to 100 nf = 2 ** lf + 1 ncyc = i4arg(3,1) preswf = i4arg(4,3) pstswf = i4arg(5,3) presw = i4arg(6,3) pstsw = i4arg(7,3) if( ltrace ) then write(0,5000) rank, nf, rank, ncyc, * preswf, pstswf, presw, pstsw 5000 format(' Rank ',i1,' Problem size: ',i3,'^',i1, * ' ',i3,' (',i2,',',i2,',',i2,',',i2,') cycles'/) end if hf = 1.0d0 / (nf - 1) call ivls(shape,nf,rank) call dvls(vh,hf,rank) call qreset(maxq) size = nf**rank uf = qalloc(size) uxctf = qalloc(size) rhsf = qalloc(size) duf = qalloc(size) mresf = duf xf = duf call clcrhs(q(rhsf),rank,shape,vh) call clcu(q(uxctf),rank,shape,vh) call dvls(q(uf),0.0d0,size) ncycle = 0 p_save = p do icyc = 1 , ncyc resnrm = mgv (rank,shape,vh, * q(uf),q(rhsf), * q(p_save),maxq-p_save+1, * preswf,pstswf,ncycle) call dvvs(q(uxctf),q(uf),q(duf),size) e2nrm = dvnrm2(q(duf),size) einfnrm = dvlinf(q(duf),size) write(0,5050) ncycle, resnrm, e2nrm, einfnrm 5050 format(' Cycle: ',i2/' |rres|_2:',1pE13.6, * ' |uxct - u|_2:',1pE13.6, * ' |uxct - u|_i:',1pE13.6) if( compute_residual ) then call cmres(q(mresf),q(uf),q(rhsf),rank,shape,vh) e2nrm = dvnrm2(q(mresf),size) einfnrm = dvlinf(q(mresf),size) write(0,5075) e2nrm, einfnrm 5075 format(' |res|_2: ',1pE13.6,' |res|_i:',1pE13.6) end if write(0,*) end do if( ltrace ) then write(0,*) p_save + qused(), ' words of storage used' write(0,*) qleft(), ' words of storage available' end if stop 100 continue write(0,*) 'usage: mg0 [ '// * ' ] [ ]' stop end c=========================================================== c Front end to clcu1, clcu2 and clcu3 ... c=========================================================== subroutine clcu(u,rank,shape,vh) implicit none integer rank real*8 u(*), vh(rank) integer shape(rank) if( rank .eq. 1 ) then call clcu1(u,vh(1),shape(1)) else if( rank .eq. 2 ) then call clcu2(u,vh(1),shape(1),shape(2)) else if( rank .eq. 3 ) then call clcu3(u,vh(1),shape(1),shape(2),shape(3)) else write(0,*) 'clcu: rank = ', rank,' not implemented' end if return end c=========================================================== c Front end to clcrhs1, clcrhs2 and clcrhs3 ... c=========================================================== subroutine clcrhs(rhs,rank,shape,vh) implicit none integer rank real*8 rhs(*), vh(rank) integer shape(rank) if( rank .eq. 1 ) then call clcrhs1(rhs,vh(1),shape(1)) else if( rank .eq. 2 ) then call clcrhs2(rhs,vh(1),shape(1),shape(2)) else if( rank .eq. 3 ) then call clcrhs3(rhs,vh(1),shape(1),shape(2),shape(3)) else write(0,*) 'clcrhs: rank = ', rank,' not implemented' end if return end c=========================================================== c Calculates exact solution u(x) = sin(l_x pi x) c and corresponding right hand side for laplacian(u) = rhs. c=========================================================== subroutine clcu1(u,h,nx) implicit none include 'mgprob.inc' integer nx real*8 u(nx) real*8 h integer i, j do i = 1 , nx u(i) = sin(pi * lx * (i - 1) * h) end do call dvlbs(u,0.0d0,nx) return end c=========================================================== subroutine clcrhs1(rhs,h,nx) implicit none include 'mgprob.inc' integer nx real*8 rhs(nx) real*8 h real*8 ui integer i do i = 1 , nx ui = sin(pi * lx * (i - 1) * h) rhs(i) = -lx * lx * pi * pi * ui end do call dvlbs(rhs,0.0d0,nx) return end c=========================================================== c Calculates exact solution u(x,y) = sin(l_x pi x) sin(l_y pi y) c and corresponding right hand side for laplacian(u) = rhs. c=========================================================== subroutine clcu2(u,h,nx,ny) implicit none include 'mgprob.inc' integer nx, ny real*8 u(nx,ny) real*8 h integer i, j do j = 1 , ny do i = 1 , nx u(i,j) = sin(pi * lx * (i - 1) * h) * * sin(pi * ly * (j - 1) * h) end do end do call dmlbs(u,0.0d0,nx,ny) return end c=========================================================== subroutine clcrhs2(rhs,h,nx,ny) implicit none include 'mgprob.inc' integer nx, ny real*8 rhs(nx,ny) real*8 h real*8 uij integer i, j do j = 1 , ny do i = 1 , nx uij = * sin(pi * lx * (i - 1) * h) * * sin(pi * ly * (j - 1) * h) rhs(i,j) = -(lx * lx + ly * ly) * * pi * pi * uij end do end do call dmlbs(rhs,0.0d0,nx,ny) return end c=========================================================== c Calculates exact solution c u(x,y) = sin(l_x pi x) sin(l_y pi y) sin(l_z pi z) c and corresponding right hand side for c laplacian(u) = rhs. c=========================================================== subroutine clcu3(u,h,nx,ny,nz) implicit none include 'mgprob.inc' integer nx, ny, nz real*8 u(nx,ny,nz) real*8 h integer i, j, k do k = 1 , nz do j = 1 , ny do i = 1 , nx u(i,j,k) = sin(pi * lx * (i - 1) * h) * * sin(pi * ly * (j - 1) * h) * * sin(pi * lz * (k - 1) * h) end do end do end do call d3lbs(u,0.0d0,nx,ny,nz) return end c=========================================================== subroutine clcrhs3(rhs,h,nx,ny,nz) implicit none include 'mgprob.inc' integer nx, ny, nz real*8 rhs(nx,ny,nz) real*8 h real*8 uijk integer i, j, k do k = 1 , nz do j = 1 , ny do i = 1 , nx uijk = * sin(pi * lx * (i - 1) * h) * * sin(pi * ly * (j - 1) * h) * * sin(pi * lz * (k - 1) * h) rhs(i,j,k) = -(lx * lx + ly * ly + lz * lz) * * pi * pi * uijk end do end do end do call d3lbs(rhs,0.0d0,nx,ny,nz) return end