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 ! # of pre-CGC relaxation sweeps c integer pstswp ! # of post-CGC relaxation sweeps c integer ncycle ! V-cycle counter (I/O) 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 ncyc = 3 ! # of V-cycles c c ! Define rhs array c . c . c . c c ncycle = 0 ! initialize ncycle, mgv auto-increments thereafter c do icyc = 1 , ncyc ! Perform ncyc v-cycles c ! mgv returns norm of residuals after one (preswp,pstswp) v-cycle c resid = mgv(rank,shape,vh,u,rhs,q,sizeq,preswp,pstswp,ncycle) c end do 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======================================================================= real*8 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 Front end to mg1v0 which invokes mg1v0 with c c preswf = presw = preswp c pstswf = pstsw = pstswp c----------------------------------------------------------------------- real*8 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 Front end to mg2v0 which invokes mg2v0 with c c preswf = presw = preswp c pstswf = pstsw = pstswp c----------------------------------------------------------------------- real*8 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 Front end to mg3v0 which invokes mg3v0 with c c preswf = presw = preswp c pstswf = pstsw = pstswp c----------------------------------------------------------------------- real*8 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 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----------------------------------------------------------------------- real*8 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(0,*) '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(0,*) 'mg1v0: ', nl, ' levels.' write(0,*) '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 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----------------------------------------------------------------------- real*8 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(0,*) 'mg2v0: Bad nxf value: ', nxf badargs = .true. end if if( checkfacts(nyf - 1,maxn0,ny0,nfact_y) .lt. 0 ) then write(0,*) '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(0,*) 'mg2v0: ', nl, ' levels.' write(0,*) '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 if( ltrace ) then write(0,*) '+++++++++++++++++++++++++' write(0,*) 'mg2v0: ncycle =', ncycle write(0,*) 'mg2v0: resnrm =', resnrm write(0,*) '+++++++++++++++++++++++++' end if mg2v0 = resnrm return end 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----------------------------------------------------------------------- real*8 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(0,*) 'mg3v0: Bad nxf value: ', nxf badargs = .true. end if if( checkfacts(nyf - 1,maxn0,ny0,nfact_y) .lt. 0 ) then write(0,*) 'mg3v0: Bad nyf value: ', nyf badargs = .true. end if if( checkfacts(nzf - 1,maxn0,nz0,nfact_z) .lt. 0 ) then write(0,*) '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(0,*) 'mg3v0: ', nl, ' levels.' write(0,*) '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 Resets pointer into main storage array. c----------------------------------------------------------------------- subroutine qreset(ceil) implicit none integer ceil integer p, pceil common / comp / p, pceil pceil = ceil p = 1 return end c----------------------------------------------------------------------- c Returns next available location in main storage array and updates c pointer. 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 Returns number of words of memory used. c----------------------------------------------------------------------- integer function qused() integer p, pceil common / comp / p, pceil qused = p return end c----------------------------------------------------------------------- c Returns number of words of available memory. c----------------------------------------------------------------------- integer function qleft() integer p, pceil common / comp / p, pceil qleft = pceil - p + 1 return end c----------------------------------------------------------------------- c Helper routine for argument checking ... c c This version doesn't demand that c c num = prime * 2^n 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(0,*) 'checkfacts: Number being factored ', num call ivdump(vfact,nfact,'factors',0) 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