c======================================================================= c ROUTINE: mgv c c Solves Poisson equation c c Laplacian[u] = rhs c c in 1-, 2- and 3-D using fixed V-cycle LCS c (linear correction scheme) MG. c c NOTE: This routine is a real*8 FUNCTION. Be sure to declare c it as such in the invoking program unit. c c TYPICAL USAGE: (I) and (O) denote input and ouptut arguments c respectively c c real*8 mgv c c integer rank c parameter ( rank = 3) ! Problem dimension (I) c integer nx, ny, nz c parameter ( nx = 33, ny = 33, nz = 33 ) c real*8 u(nx,ny,nz) ! Soln (O) c real*8 rhs(nx,ny,nz) ! Rhs (I) c integer sizeq c parameter ( sizeq = 3 * nx * ny * nz ) ! Size of work array (I) c real*8 q(sizeq) ! Work array (I) c real*8 vh(rank) ! Mesh spacings (I) c integer shape(rank) ! Mesh dimensions (I) c integer preswp, pstswp, ncycle ! Multi-grid parameters (I) c c c vh(1) = h ! This example assumes mesh spacings in c ! all coordinate direction = h c vh(2) = h c vh(3) = h c shape(1) = nx c shape(2) = ny c shape(3) = nz c c preswp = 3 ! # of pre-CGC sweeps c pstswp = 3 ! # of post-CGC sweeps c ncycle = 1 ! # of V-cycles c c ! Define rhs array c . c . c . c c ! mgv returns norm of residuals after ncycle (preswp,pstswp) c ! v-cycles. c c resid = mgv(rank,shape,vh,u,rhs,q,sizeq,preswp,pstswp,ncycle) c . c . c . c c c IMPLEMENTATION NOTES: This routine is a front end to mg1v0, c mg2v0 and mg3v0. It invokes those routines with c c preswf = presw = preswp c pstswf = pstsw = pstswp c======================================================================= double precision function mgv * (rank,shape,vh,uf,rhsf, * q,sizeq, * preswp,pstswp,ncycle) implicit none real*8 mg1v0, mg2v0, mg3v0 integer rank, sizeq integer shape(rank) real*8 vh(rank), uf(*), rhsf(*), * q(sizeq) integer preswp, pstswp, ncycle logical ltrace parameter ( ltrace = .false. ) if( ltrace ) then write(0,*) 'mgv: Parameter dump' write(0,*) 'rank, preswp, pstswp, ncycle' write(0,*) rank, preswp, pstswp, ncycle call ivdump(shape,rank,'h',0) call dvdump(vh,rank,'h',0) end if if( rank .eq. 1 ) then mgv = mg1v0(uf,rhsf,shape(1), vh(1), * q,sizeq, * preswp,pstswp,preswp,pstswp,ncycle) else if( rank .eq. 2 ) then mgv = mg2v0(uf,rhsf,shape(1),shape(2), vh(1), * q,sizeq, * preswp,pstswp,preswp,pstswp,ncycle) else if( rank .eq. 3 ) then mgv = mg3v0(uf,rhsf,shape(1),shape(2),shape(3),vh(1), * q,sizeq, * preswp,pstswp,preswp,pstswp,ncycle) else write(0,*) 'mgv: Invalid rank: ', rank mgv = -1.0d0 end if return end c----------------------------------------------------------------------- c c Front end to mg1v0 which invokes mg1v0 with c c preswf = presw = preswp c pstswf = pstsw = pstswp c c----------------------------------------------------------------------- double precision function mg1v * (uf,rhsf,nxf, hf, q,sizeq, * preswp,pstswp,ncycle) implicit none real*8 mg1v0 integer nxf, sizeq real*8 uf(nxf), rhsf(nxf), * q(sizeq) real*8 hf integer preswp, pstswp, ncycle mg1v = mg1v0(uf,rhsf,nxf, hf,q,sizeq, * preswp,pstswp,preswp,pstswp,ncycle) return end c----------------------------------------------------------------------- c c Front end to mg2v0 which invokes mg2v0 with c c preswf = presw = preswp c pstswf = pstsw = pstswp c c----------------------------------------------------------------------- double precision function mg2v * (uf,rhsf,nxf,nyf, hf, q,sizeq, * preswp,pstswp,ncycle) implicit none real*8 mg2v0 integer nxf, nyf, sizeq real*8 uf(nxf,nyf), rhsf(nxf,nyf), * q(sizeq) real*8 hf integer preswp, pstswp, ncycle mg2v = mg2v0(uf,rhsf,nxf,nyf,hf,q,sizeq, * preswp,pstswp,preswp,pstswp,ncycle) return end c----------------------------------------------------------------------- c c Front end to mg3v0 which invokes mg3v0 with c c preswf = presw = preswp c pstswf = pstsw = pstswp c c----------------------------------------------------------------------- double precision function mg3v * (uf,rhsf,nxf,nyf,nzf, hf, q,sizeq, * preswp,pstswp,ncycle) implicit none real*8 mg3v0 integer nxf, nyf, nzf, sizeq real*8 uf(nxf,nyf,nzf), rhsf(nxf,nyf,nzf), * q(sizeq) real*8 hf integer preswp, pstswp, ncycle mg3v = mg3v0(uf,rhsf,nxf,nyf,nzf,hf,q,sizeq, * preswp,pstswp,preswp,pstswp,ncycle) return end c----------------------------------------------------------------------- c c Basic 1d v-cycle multi grid solver for c c Laplacian[u(x)] = rhs(x) c c This version allows nxf, nyf which are of form c 2^p * + 1 c c where .le. maxn0 defined in the parameter c statement below c c----------------------------------------------------------------------- double precision function mg1v0 * (uf,rhsf,nxf, hf, q, sizeq, * preswf,pstswf,presw,pstsw,ncycle) implicit none real*8 relaxc1 integer qalloc, lilog2 logical lp2np1 integer nxf, sizeq real*8 uf(nxf), rhsf(nxf), * q(sizeq) real*8 hf integer preswf, pstswf, presw, pstsw, * ncycle c c Multi-grid parameters, data structures and pointers c real*8 epsic parameter ( epsic = 1.0d-10 ) integer lmin, lmax parameter ( lmin = 1, lmax = 20 ) integer nx(lmin:lmax) real*8 h(lmin:lmax) integer u(lmin:lmax), rhs(lmin:lmax), * mres(lmin:lmax), du(lmin:lmax) integer imposs parameter ( imposs = -2 000 000 000 ) integer p, pceil common / comp / p, pceil c c For parameter (nxf) checking c integer checkfacts logical badargs integer nfact_x, nx0 integer maxn0 parameter ( maxn0 = 17 ) c c Locals ... c integer lxf, * l, lf, lc, nl, * sweep real*8 resnrm logical ltrace parameter ( ltrace = .false. ) save mg1v0 = -1.0d0 c c Parameter check ... c if( ltrace ) then write(0,*) '>>> mg1v0: Parameter dump: ' write(0,*) nxf, hf write(0,*) preswf, pstswf, presw, pstsw, ncycle end if badargs = .false. if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then write(*,*) 'mg1v0: Bad nxf value: ', nxf badargs = .true. end if if( badargs) then return end if nl = nfact_x + 1 lc = 1 lf = lc + nl - 1 nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1 if( ltrace ) then write(*,*) 'mg1v0: ', nl, ' levels.' write(*,*) 'mg1v0: Coarse grid ', nx0 end if if( ncycle .eq. 0 ) then c c Set up grid structure, allocate storage ... c call qreset(sizeq) nx(lf) = nxf h(lf) = hf u(lf) = imposs rhs(lf) = imposs mres(lf) = qalloc(nx(lf)) du(lf) = mres(lf) do 1000 l = lf - 1 , lc , -1 nx(l) = (nx(l+1) - 1) / 2 + 1 h(l) = 2.0d0 * h(l+1) u(l) = qalloc(nx(l)) rhs(l) = qalloc(nx(l)) mres(l) = qalloc(nx(l)) du(l) = mres(l) 1000 continue if( ltrace ) then call ivdump(nx,nl,'nx',0) call dvdump(h,nl,'h',0) call ivdump(u,nl,'u pointers',0) call ivdump(rhs,nl,'rhs pointers',0) call ivdump(mres,nl,'mres pointers',0) call ivdump(du,nl,'du pointers',0) end if end if c c Cycle from fine to coarse ... c do l = lf , lc , -1 if( l .ne. lf .or. ncycle .eq. 0 ) then if( l .eq. lc ) then 100 continue resnrm = relaxc1(q(u(l)),q(rhs(l)), * nx(l),h(l)) if( ltrace ) then call trace_res('mg1v0: pre sweep',l,resnrm) call cres1_trace('mg1v0: pre lc',60, & q(mres(l)),q(u(l)),q(rhs(l)), & nx(l),h(l)) end if if( resnrm .gt. epsic ) go to 100 else if( l .eq. lf ) then do sweep = 1 , preswf resnrm = relaxc1(uf,rhsf, * nx(l),h(l)) call cmres1(q(mres(l)),uf,rhsf, * nx(l),h(l)) if( ltrace ) then call trace_res('mg1v0: pre sweep',l,resnrm) call cres1_trace('mg1v0: pre lf',60, & q(mres(l)),uf,rhsf, & nx(l),h(l)) end if end do else do sweep = 1 , presw resnrm = relaxc1(q(u(l)),q(rhs(l)), * nx(l),h(l)) if( ltrace ) then call trace_res('mg1v0: pre sweep',l,resnrm) call cres1_trace('mg1v0: pre inter',60, & q(mres(l)),q(u(l)),q(rhs(l)), & nx(l),h(l)) end if end do end if end if c c Set up coarse grid problem ... c if( l .eq. lf ) then call cmres1(q(mres(l)),uf,rhsf,nx(l),h(l)) else if( l .ne. lc ) then call cmres1(q(mres(l)),q(u(l)),q(rhs(l)), * nx(l),h(l)) end if if( l .ne. lc ) then call dvrshw(q(mres(l)),q(rhs(l-1)), * nx(l-1)) call dvls(q(u(l-1)),0.0d0, * nx(l-1)) end if end do c c Cycle from coarse to fine ... c do l = lc + 1 , lf c c Interpolate correction ... c call dvlint(q(u(l-1)),q(du(l)),nx(l-1)) if( l .eq. lf ) then call dvva(uf,q(du(l)),uf,nx(l)) else call dvva(q(u(l)),q(du(l)),q(u(l)),nx(l)) end if if( l .eq. lf ) then do sweep = 1 , pstswf resnrm = relaxc1(uf,rhsf,nx(l),h(l)) if( ltrace ) then call trace_res('mg1v0: post sweep',l,resnrm) call cres1_trace('mg1v0: post lf',60, & q(mres(l)),uf,rhsf, & nx(l),h(l)) end if end do else do sweep = 1 , pstsw resnrm = relaxc1(q(u(l)),q(rhs(l)), * nx(l),h(l)) if( ltrace ) then call trace_res('mg1v0: post sweep',l,resnrm) call cres1_trace('mg1v0: post other',60, & q(mres(l)),q(u(l)),q(rhs(l)), & nx(l),h(l)) end if end do end if end do ncycle = ncycle + 1 mg1v0 = resnrm return end c----------------------------------------------------------------------- c c Basic 2d v-cycle multi grid solver for c c Laplacian[u(x,y)] = rhs(x,y) c c on rectangular domain. c c This version allows nxf, nyf which are of form c 2^p * + 1 c c where .le. maxn0 defined in the parameter c statement below c c----------------------------------------------------------------------- double precision function mg2v0 * (uf,rhsf,nxf,nyf, hf, q, sizeq, * preswf,pstswf,presw,pstsw,ncycle) implicit none real*8 relaxc2 integer qalloc, lilog2 logical lp2np1 integer nxf, nyf, sizeq real*8 uf(nxf,nyf), rhsf(nxf,nyf), * q(sizeq) real*8 hf integer preswf, pstswf, presw, pstsw, * ncycle c c Multi-grid parameters, data structures and pointers c real*8 epsic parameter ( epsic = 1.0d-10 ) integer lmin, lmax parameter ( lmin = 1, lmax = 20 ) integer nx(lmin:lmax), ny(lmin:lmax) real*8 h(lmin:lmax) integer u(lmin:lmax), rhs(lmin:lmax), * mres(lmin:lmax), du(lmin:lmax) integer imposs parameter ( imposs = -2 000 000 000 ) integer p, pceil common / comp / p, pceil c c For parameter (nxf, nyf) checking c integer checkfacts logical badargs integer nfact_x, nx0 integer nfact_y, ny0 integer maxn0 parameter ( maxn0 = 17 ) c c Locals ... c integer lxf, lyf, * l, lf, lc, nl, * sweep real*8 resnrm logical ltrace parameter ( ltrace = .false. ) save mg2v0 = -1.0d0 c c Parameter check ... c if( ltrace ) then write(0,*) '>>> mg2v0: Parameter dump: ' write(0,*) nxf, nyf, hf write(0,*) preswf, pstswf, presw, pstsw, ncycle end if badargs = .false. if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then write(*,*) 'mg2v0: Bad nxf value: ', nxf badargs = .true. end if if( checkfacts(nyf - 1,maxn0,ny0,nfact_y) .lt. 0 ) then write(*,*) 'mg2v0: Bad nyf value: ', nyf badargs = .true. end if if( badargs) then return end if nl = min(nfact_x,nfact_y) + 1 lc = 1 lf = lc + nl - 1 nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1 ny0 = (nyf - 1) / (2 ** (nl - 1)) + 1 if( ltrace ) then write(*,*) 'mg2v0: ', nl, ' levels.' write(*,*) 'mg2v0: Coarse grid ', nx0, ' x ', ny0 end if if( ncycle .eq. 0 ) then c c Set up grid structure, allocate storage ... c call qreset(sizeq) nx(lf) = nxf ny(lf) = nyf h(lf) = hf u(lf) = imposs rhs(lf) = imposs mres(lf) = qalloc(nx(lf)*ny(lf)) du(lf) = mres(lf) do 1000 l = lf - 1 , lc , -1 nx(l) = (nx(l+1) - 1) / 2 + 1 ny(l) = (ny(l+1) - 1) / 2 + 1 h(l) = 2.0d0 * h(l+1) u(l) = qalloc(nx(l)*ny(l)) rhs(l) = qalloc(nx(l)*ny(l)) mres(l) = qalloc(nx(l)*ny(l)) du(l) = mres(l) 1000 continue if( ltrace ) then call ivdump(nx,nl,'nx',0) call ivdump(ny,nl,'ny',0) call dvdump(h,nl,'h',0) call ivdump(u,nl,'u pointers',0) call ivdump(rhs,nl,'rhs pointers',0) call ivdump(mres,nl,'mres pointers',0) call ivdump(du,nl,'du pointers',0) end if end if c c Cycle from fine to coarse ... c do l = lf , lc , -1 if( l .ne. lf .or. ncycle .eq. 0 ) then if( l .eq. lc ) then 100 continue resnrm = relaxc2(q(u(l)),q(rhs(l)), * nx(l),ny(l),h(l)) if( ltrace ) then call trace_res('mg2v0: pre sweep',l,resnrm) end if if( resnrm .gt. epsic ) go to 100 else if( l .eq. lf ) then do sweep = 1 , preswf resnrm = relaxc2(uf,rhsf, * nx(l),ny(l),h(l)) call cmres2(q(mres(l)),uf,rhsf, * nx(l),ny(l),h(l)) if( ltrace ) then call trace_res('mg2v0: pre sweep',l,resnrm) end if end do else do sweep = 1 , presw resnrm = relaxc2(q(u(l)),q(rhs(l)), * nx(l),ny(l),h(l)) if( ltrace ) then call trace_res('mg2v0: pre sweep',l,resnrm) end if end do end if end if c c Set up coarse grid problem ... c if( l .eq. lf ) then call cmres2(q(mres(l)),uf,rhsf,nx(l),ny(l),h(l)) else if( l .ne. lc ) then call cmres2(q(mres(l)),q(u(l)),q(rhs(l)), * nx(l),ny(l),h(l)) end if if( l .ne. lc ) then call dmrshw(q(mres(l)),q(rhs(l-1)), * nx(l-1),ny(l-1)) call dmls(q(u(l-1)),0.0d0, * nx(l-1),ny(l-1)) end if end do c c Cycle from coarse to fine ... c do l = lc + 1 , lf c c Interpolate correction ... c call dmlint(q(u(l-1)),q(du(l)),nx(l-1),ny(l-1)) if( l .eq. lf ) then call dmma(uf,q(du(l)),uf,nx(l),ny(l)) else call dmma(q(u(l)),q(du(l)),q(u(l)),nx(l),ny(l)) end if if( l .eq. lf ) then do sweep = 1 , pstswf resnrm = relaxc2(uf,rhsf,nx(l),ny(l),h(l)) if( ltrace ) then call trace_res('mg2v0: post sweep',l,resnrm) end if end do else do sweep = 1 , pstsw resnrm = relaxc2(q(u(l)),q(rhs(l)), * nx(l),ny(l),h(l)) if( ltrace ) then call trace_res('mg2v0: post sweep',l,resnrm) end if end do end if end do ncycle = ncycle + 1 mg2v0 = resnrm return end c----------------------------------------------------------------------- c c Basic 3d v-cycle multi grid solver for c c Laplacian[u(x,y,z)] = rhs(x,y,z) c c on rectangular domain. c c This version allows nxf, nyf, nzy which are of form c 2^p * + 1 c c where .le. maxn0 defined in the parameter c statement below c c----------------------------------------------------------------------- double precision function mg3v0 * (uf,rhsf,nxf,nyf,nzf, hf, q, sizeq, * preswf,pstswf,presw,pstsw,ncycle) implicit none real*8 d3nrm2, relaxc3 integer qalloc, lilog2 logical lp2np1 integer nxf, nyf, nzf, sizeq real*8 uf(nxf,nyf,nzf), rhsf(nxf,nyf,nzf), * q(sizeq) real*8 hf integer preswf, pstswf, presw, pstsw, * ncycle c c Multi-grid parameters, data structures and pointers c real*8 epsic parameter ( epsic = 1.0d-10 ) integer lmin, lmax parameter ( lmin = 1, lmax = 20 ) integer nx(lmin:lmax), ny(lmin:lmax), * nz(lmin:lmax) real*8 h(lmin:lmax) integer u(lmin:lmax), rhs(lmin:lmax), * mres(lmin:lmax), du(lmin:lmax) integer imposs parameter ( imposs = -2 000 000 000 ) integer p, pceil common / comp / p, pceil c c For parameter (nxf, nyf, nzf) checking c integer checkfacts logical badargs integer nfact_x, nx0 integer nfact_y, ny0 integer nfact_z, nz0 integer maxn0 parameter ( maxn0 = 17 ) c c Locals ... c integer lxf, lyf, lzf, * l, lf, lc, nl, * sweep real*8 resnrm logical ltrace parameter ( ltrace = .false. ) save mg3v0 = -1.0d0 c c Parameter check ... c if( ltrace ) then write(0,*) '>>> mg3v0: Parameter dump: ' write(0,*) nxf, nyf, nzf, hf write(0,*) preswf, pstswf, presw, pstsw, ncycle end if badargs = .false. if( checkfacts(nxf - 1,maxn0,nx0,nfact_x) .lt. 0 ) then write(*,*) 'mg3v0: Bad nxf value: ', nxf badargs = .true. end if if( checkfacts(nyf - 1,maxn0,ny0,nfact_y) .lt. 0 ) then write(*,*) 'mg3v0: Bad nyf value: ', nyf badargs = .true. end if if( checkfacts(nzf - 1,maxn0,nz0,nfact_z) .lt. 0 ) then write(*,*) 'mg3v0: Bad nzf value: ', nzf badargs = .true. end if if( badargs) then return end if nl = min(nfact_x,nfact_y,nfact_z) + 1 lc = 1 lf = lc + nl - 1 nx0 = (nxf - 1) / (2 ** (nl - 1)) + 1 ny0 = (nyf - 1) / (2 ** (nl - 1)) + 1 nz0 = (nzf - 1) / (2 ** (nl - 1)) + 1 if( ltrace ) then write(*,*) 'mg3v0: ', nl, ' levels.' write(*,*) 'mg3v0: Coarse grid ', nx0, ' x ', ny0, * ' x ', nz0 end if if( ncycle .eq. 0 ) then c c Set up grid structure, allocate storage ... c call qreset(sizeq) nx(lf) = nxf ny(lf) = nyf nz(lf) = nzf h(lf) = hf u(lf) = imposs rhs(lf) = imposs mres(lf) = qalloc(nx(lf)*ny(lf)*nz(lf)) du(lf) = mres(lf) do 1000 l = lf - 1 , lc , -1 nx(l) = (nx(l+1) - 1) / 2 + 1 ny(l) = (ny(l+1) - 1) / 2 + 1 nz(l) = (nz(l+1) - 1) / 2 + 1 h(l) = 2.0d0 * h(l+1) u(l) = qalloc(nx(l)*ny(l)*nz(l)) rhs(l) = qalloc(nx(l)*ny(l)*nz(l)) mres(l) = qalloc(nx(l)*ny(l)*nz(l)) du(l) = mres(l) 1000 continue if( ltrace ) then call ivdump(nx,nl,'nx',0) call ivdump(ny,nl,'ny',0) call ivdump(nz,nl,'nz',0) call dvdump(h,nl,'h',0) call ivdump(u,nl,'u pointers',0) call ivdump(rhs,nl,'rhs pointers',0) call ivdump(mres,nl,'mres pointers',0) call ivdump(du,nl,'du pointers',0) end if end if c c Cycle from fine to coarse ... c do l = lf , lc , -1 if( l .ne. lf .or. ncycle .eq. 0 ) then if( l .eq. lc ) then 100 continue resnrm = relaxc3(q(u(l)),q(rhs(l)), * nx(l),ny(l),nz(l),h(l)) if( ltrace ) then call trace_res('mg3v0: pre sweep',l,resnrm) end if if( resnrm .gt. epsic ) go to 100 else if( l .eq. lf ) then do sweep = 1 , preswf resnrm = relaxc3(uf,rhsf, * nx(l),ny(l),nz(l),h(l)) call cmres3(q(mres(l)),uf,rhsf, * nx(l),ny(l),nz(l),h(l)) if( ltrace ) then call trace_res('mg3v0: pre sweep',l,resnrm) end if end do else do sweep = 1 , presw resnrm = relaxc3(q(u(l)),q(rhs(l)), * nx(l),ny(l),nz(l),h(l)) if( ltrace ) then call trace_res('mg3v0: pre sweep',l,resnrm) end if end do end if end if c c Set up coarse grid problem ... c if( l .eq. lf ) then call cmres3(q(mres(l)),uf,rhsf,nx(l),ny(l),nz(l),h(l)) else if( l .ne. lc ) then call cmres3(q(mres(l)),q(u(l)),q(rhs(l)), * nx(l),ny(l),nz(l),h(l)) end if if( l .ne. lc ) then call d3rshw(q(mres(l)),q(rhs(l-1)), * nx(l-1),ny(l-1),nz(l-1)) call d3ls(q(u(l-1)),0.0d0, * nx(l-1),ny(l-1),nz(l-1)) end if end do c c Cycle from coarse to fine ... c do l = lc + 1 , lf c c Interpolate correction ... c call d3lint(q(u(l-1)),q(du(l)),nx(l-1),ny(l-1),nz(l-1)) if( l .eq. lf ) then call d33a(uf,q(du(l)),uf,nx(l),ny(l),nz(l)) else call d33a(q(u(l)),q(du(l)),q(u(l)),nx(l),ny(l),nz(l)) end if if( l .eq. lf ) then do sweep = 1 , pstswf resnrm = relaxc3(uf,rhsf,nx(l),ny(l),nz(l),h(l)) if( ltrace ) then call trace_res('mg3v0: post sweep',l,resnrm) end if end do else do sweep = 1 , pstsw resnrm = relaxc3(q(u(l)),q(rhs(l)), * nx(l),ny(l),nz(l),h(l)) if( ltrace ) then call trace_res('mg3v0: post sweep',l,resnrm) end if end do end if end do ncycle = ncycle + 1 mg3v0 = resnrm return end c----------------------------------------------------------------------- c c Resets pointer into main storage array. c c----------------------------------------------------------------------- subroutine qreset(ceil) implicit none integer ceil integer p, pceil common / comp / p, pceil pceil = ceil p = 1 return end c----------------------------------------------------------------------- c c Returns next available location in main storage array and updates c pointer. c c----------------------------------------------------------------------- integer function qalloc(size) integer size integer p, pceil common / comp / p, pceil if( size .gt. 0 ) then if( p + size .le. pceil ) then qalloc = p p = p + size else write(0,1000) int(size / 1.0d6 + 0.5d0), & int(pceil / 1.0d6 + 0.5d0) 1000 format(' >>> qalloc: Out of memory'/ & ' >>> qalloc: Requested block size:', i5, & ' Mbytes'/ & ' >>> qalloc: Total available: ', i5, & ' Mbytes') stop end if else write(0,*) '>>> qalloc: size ', size, ' <= 0.' qalloc = -1 end if return end c----------------------------------------------------------------------- c c Returns number of words of memory used. c c----------------------------------------------------------------------- integer function qused() integer p, pceil common / comp / p, pceil qused = p return end c----------------------------------------------------------------------- c c Returns number of words of available memory. c c----------------------------------------------------------------------- integer function qleft() integer p, pceil common / comp / p, pceil qleft = pceil - p + 1 return end c----------------------------------------------------------------------- c c Tracing routine c c----------------------------------------------------------------------- subroutine trace_res(string,l,resnrm) implicit none character*(*) string integer l real*8 resnrm write(0,1000) string, l, resnrm 1000 format(a,' Level: ',i2,' |res.|: ',1pe12.3) return end c----------------------------------------------------------------------- c c Does one-step 1d Point-GS relaxation and returns norm of c running ('instantaneous') residuals---checkerboard ordering. c c----------------------------------------------------------------------- double precision function relaxc1(u,rhs,nx,h) implicit none integer nx real*8 u(nx), rhs(nx) real*8 h real*8 hm2, jacel, rresl integer i, pass relaxc1 = 0.0d0 hm2 = 1.0d0 / (h * h) jacel = -2.0d0 * hm2 do pass = 1 , 2 do i = pass + 1 , nx - 1 , 2 rresl = hm2 * (u(i+1) + u(i-1) - 2.0d0*u(i)) - * rhs(i) u(i) = u(i) - rresl / jacel relaxc1 = relaxc1 + rresl * rresl end do end do relaxc1 = sqrt(relaxc1 / nx) return end c----------------------------------------------------------------------- c c Computes negative of 1d residuals ... c c----------------------------------------------------------------------- subroutine cmres1(mres,u,rhs,nx,h) implicit none integer nx real*8 mres(nx), u(nx), * rhs(nx) real*8 h integer i real*8 hm2 hm2 = 1.0d0 / (h * h) do i = 2 , nx - 1 mres(i) = rhs(i) - * (hm2 * (u(i-1) - 2.0d0 * u(i) + u(i+1))) end do call dvlbs(mres,0.0d0,nx) return end c----------------------------------------------------------------------- c c Does one-step 2d Point-GS relaxation and returns norm of c running ('instantaneous') residuals---checkerboard ordering. c c----------------------------------------------------------------------- double precision function relaxc2(u,rhs,nx,ny,h) implicit none integer nx, ny real*8 u(nx,ny), rhs(nx,ny) real*8 h real*8 hm2, mjelm1, rresl integer i, j, * isw, jsw, pass relaxc2 = 0.0d0 hm2 = 1.0d0 / (h * h) mjelm1 = h * h / 4.0d0 jsw = 1 do pass = 1 , 2 isw = jsw do j = 2 , ny - 1 do i = isw + 1 , nx - 1, 2 rresl = hm2 * (u(i+1,j) + u(i-1,j) + * u(i,j+1) + u(i,j-1) - * 4.0d0 * u(i,j)) - * rhs(i,j) u(i,j) = u(i,j) + mjelm1 * rresl relaxc2 = relaxc2 + rresl * rresl end do isw = 3 - isw end do jsw = 3 - jsw end do relaxc2 = sqrt(relaxc2 / (nx * ny)) return end c----------------------------------------------------------------------- c c Computes negative of 2d residuals ... c c----------------------------------------------------------------------- subroutine cmres2(mres,u,rhs,nx,ny,h) implicit none integer nx, ny real*8 mres(nx,ny), u(nx,ny), * rhs(nx,ny) real*8 h integer i, j real*8 hm2 hm2 = 1.0d0 / (h * h) do j = 2 , ny - 1 do i = 2 , nx - 1 mres(i,j) = rhs(i,j) - * hm2 * (u(i+1,j) + u(i-1,j) + * u(i,j+1) + u(i,j-1) - * 4.0d0*u(i,j)) end do end do call dmlbs(mres,0.0d0,nx,ny) return end c----------------------------------------------------------------------- c c Does one-step 3d Point-GS relaxation and returns norm of c running ('instantaneous') residuals---checkerboard ordering. c c----------------------------------------------------------------------- double precision function relaxc3(u,rhs,nx,ny,nz,h) implicit none integer nx, ny, nz real*8 u(nx,ny,nz), rhs(nx,ny,nz) real*8 h real*8 hm2, mjelm1, rresl integer i, j, k, * isw, ksw, pass relaxc3 = 0.0d0 hm2 = 1.0d0 / (h * h) mjelm1 = h * h / 6.0d0 ksw = 2 do pass = 1 , 2 isw = ksw do k = 2 , nz - 1 do j = 2 , ny - 1 do i = isw + 1 , nx - 1, 2 rresl = hm2 * (u(i+1,j,k) + u(i-1,j,k) + * u(i,j+1,k) + u(i,j-1,k) + * u(i,j,k+1) + u(i,j,k-1) - * 6.0d0 * u(i,j,k)) - * rhs(i,j,k) u(i,j,k) = u(i,j,k) + mjelm1 * rresl relaxc3 = relaxc3 + rresl * rresl end do isw = 3 - isw end do end do ksw = 3 - ksw end do relaxc3 = sqrt(relaxc3 / (nx * ny * nz)) return end c----------------------------------------------------------------------- c c Computes negative of 3d residuals ... c c----------------------------------------------------------------------- subroutine cmres3(mres,u,rhs,nx,ny,nz,h) implicit none integer nx, ny, nz real*8 mres(nx,ny,nz), u(nx,ny,nz), * rhs(nx,ny,nz) real*8 h integer i, j, k real*8 hm2 hm2 = 1.0d0 / (h * h) do k = 2 , nz - 1 do j = 2 , ny - 1 do i = 2 , nx - 1 mres(i,j,k) = rhs(i,j,k) - * hm2 * (u(i+1,j,k) + u(i-1,j,k) + * u(i,j+1,k) + u(i,j-1,k) + * u(i,j,k+1) + u(i,j,k-1) - * 6.0d0*u(i,j,k)) end do end do end do call d3lbs(mres,0.0d0,nx,ny,nz) return end c----------------------------------------------------------------------- c c Helper routine for argument checking ... c c This version doesn't demand that c c num = prime * 2^n c c----------------------------------------------------------------------- integer function checkfacts(num,maxbase,base,npow2) implicit none integer ivprimedecomp integer num, maxbase, base, npow2 integer maxnfact parameter ( maxnfact = 100 ) integer vfact(maxnfact) integer nfact, i logical ltrace parameter ( ltrace = .false. ) if( ivprimedecomp(num,vfact,nfact) .gt. 0 ) then npow2 = 0 do i = 1 , nfact - 1 if( vfact(i) .ne. 2 ) goto 100 npow2 = npow2 + 1 end do 100 continue base = 1 do i = npow2 + 1 , nfact base = base * vfact(i) end do if( base .le. maxbase ) then checkfacts = npow2 else checkfacts = -1 end if else checkfacts = -1 end if if( ltrace ) then write(*,*) 'checkfacts: Number being factored ', num call ivdump(vfact,nfact,'factors',6) end if return end c----------------------------------------------------------------------- c Computes and dumps true residuals c----------------------------------------------------------------------- subroutine cres1_trace(label,unit,res,u,rhs,nx,h) implicit none character*(*) label integer nx, unit real*8 res(*), u(*), rhs(*) real*8 h integer shape(1) real*8 vh(1) shape(1) = nx vh(1) = h call cres(res,u,rhs,1,shape,vh) write(unit,1000) label 1000 format('cres1_trace: ',a) call dvdump(u,nx,'u',unit) call dvdump(res,nx,'residual',unit) return end c----------------------------------------------------------------------- c Computes true (non-negated) residual. c----------------------------------------------------------------------- subroutine cres(res,u,rhs,rank,shape,vh) implicit none integer rank real*8 res(*), u(*), rhs(*), * vh(rank) integer shape(rank) call cmres(res,u,rhs,rank,shape,vh) if( rank .eq. 1 ) then call dvsm(res,-1.0d0,res,shape(1)) else if( rank .eq. 2 ) then call dvsm(res,-1.0d0,res,shape(1)*shape(2)) else if( rank .eq. 3 ) then call dvsm(res,-1.0d0,res,shape(1)*shape(2)*shape(3)) else write(0,*) 'cmres: rank = ', rank,' not implemented' end if return end c----------------------------------------------------------------------- c Front end to cmres1, cmres2, cmres3 ... c----------------------------------------------------------------------- subroutine cmres(mres,u,rhs,rank,shape,vh) implicit none integer rank real*8 mres(*), u(*), rhs(*), * vh(rank) integer shape(rank) if( rank .eq. 1 ) then call cmres1(mres,u,rhs,shape(1), vh(1)) else if( rank .eq. 2 ) then call cmres2(mres,u,rhs,shape(1),shape(2), vh(1)) else if( rank .eq. 3 ) then call cmres3(mres,u,rhs,shape(1),shape(2),shape(3),vh(1)) else write(0,*) 'cmres: rank = ', rank,' not implemented' end if return end c----------------------------------------------------------------------- c Eventually should be incorporated into dmatlib.f c c Multiply matrix with scalar ... c----------------------------------------------------------------------- subroutine dmsm(a1,s1,a2,d1,d2) implicit none integer d1, d2 real*8 a1(d1,d2), a2(d1,d2) real*8 s1 call dvsm(a1,s1,a2,d1*d2) return end