! fishpack.f 23 February 1998 ! ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! * * ! * f i s h p a k * ! * * ! * a package of fortran subprograms for the solution of * ! * separable elliptic partial differential equations * ! * * ! * (version 3.1 , october 1980) * ! * by * ! * john adams, paul swarztrauber and roland sweet * ! * of * ! * the national center for atmospheric research * ! * boulder, colorado (80307) u.s.a. * ! * which is sponsored by * ! * the national science foundation * ! * * ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! function bcrh(xll,xrr,iz,c,a,bh,f,sgn) ! !*********************************************************************** ! dimension a(1) ,c(1) ,bh(1) ! external f ! common /ccblk/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik xl = xll xr = xrr dx = .5*abs(xr-xl) 101 x = .5*(xl+xr) if (sgn*f(x,iz,c,a,bh)) 103,105,102 102 xr = x go to 104 103 xl = x 104 dx = .5*dx if (dx-cnv) 105,105,101 105 bcrh = .5*(xl+xr) return end subroutine blktr1 ( n, an, bn, cn, m, am, bm, cm, idimy, y, b, & w1, w2, w3, wd, ww, wu, prdct, cprdct ) ! !*********************************************************************** ! ! ! WARNING: BN is nowhere used. Is this an error? ! !! BLKTR1 solves the linear system ! ! b contains the roots of all the b polynomials ! ! w1,w2,w3,wd,ww,wu are all working arrays ! ! prdct is either prodp or prod depending on whether the boundary ! conditions in the m direction are periodic or not ! ! cprdct is either cprodp or cprod which are the complex versions ! of prodp and prod. these are called in the event that some ! of the roots of the b sub p polynomial are complex ! integer idimy integer m integer n ! real am(m) real an(n) real b(*) real bm(m) real bn(n) real cm(m) real cn(n) real cnv real eps integer ik integer k integer ncmplx integer nm integer npp real w1(1) real w2(1) real w3(1) real wd(1) real wu(1) real ww(1) real y(idimy,n) ! external cprdct external prdct ! common /cblkt/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik ! ! ! Begin reduction phase. ! kdo = k-1 do l=1,kdo ir = l-1 i2 = 2**ir i1 = i2/2 i3 = i2+i1 i4 = i2+i2 irm1 = ir-1 call indxb (i2,ir,im2,nm2) call indxb (i1,irm1,im3,nm3) call indxb (i3,irm1,im1,nm1) call prdct (nm2,b(im2),nm3,b(im3),nm1,b(im1),0,dum,y(1,i2),w3, & m,am,bm,cm,wd,ww,wu) if = 2**k do 108 i=i4,if,i4 if (i-nm) 101,101,108 101 ipi1 = i+i1 ipi2 = i+i2 ipi3 = i+i3 call indxc (i,ir,idxc,nc) if (i-if) 102,108,108 102 continue call indxa (i,ir,idxa,na) call indxb (i-i1,irm1,im1,nm1) call indxb (ipi2,ir,ip2,np2) call indxb (ipi1,irm1,ip1,np1) call indxb (ipi3,irm1,ip3,np3) call prdct (nm1,b(im1),0,dum,0,dum,na,an(idxa),w3,w1,m,am, & bm,cm,wd,ww,wu) if (ipi2-nm) 105,105,103 103 do j=1,m w3(j) = 0. w2(j) = 0. end do go to 106 105 continue call prdct (np2,b(ip2),np1,b(ip1),np3,b(ip3),0,dum, & y(1,ipi2),w3,m,am,bm,cm,wd,ww,wu) call prdct (np1,b(ip1),0,dum,0,dum,nc,cn(idxc),w3,w2,m,am, & bm,cm,wd,ww,wu) 106 continue do j=1,m y(j,i) = w1(j)+w2(j)+y(j,i) end do 108 continue end do if (npp) 132,110,132 ! ! The periodic case is treated using the capacitance matrix method ! 110 if = 2**k i = if/2 i1 = i/2 call indxb (i-i1,k-2,im1,nm1) call indxb (i+i1,k-2,ip1,np1) call indxb (i,k-1,iz,nz) call prdct (nz,b(iz),nm1,b(im1),np1,b(ip1),0,dum,y(1,i),w1,m,am, & bm,cm,wd,ww,wu) izr = i do j=1,m w2(j) = w1(j) end do do 113 ll=2,k l = k-ll+1 ir = l-1 i2 = 2**ir i1 = i2/2 i = i2 call indxc (i,ir,idxc,nc) call indxb (i,ir,iz,nz) call indxb (i-i1,ir-1,im1,nm1) call indxb (i+i1,ir-1,ip1,np1) call prdct (np1,b(ip1),0,dum,0,dum,nc,cn(idxc),w1,w1,m,am,bm, & cm,wd,ww,wu) do j=1,m w1(j) = y(j,i)+w1(j) end do call prdct (nz,b(iz),nm1,b(im1),np1,b(ip1),0,dum,w1,w1,m,am, & bm,cm,wd,ww,wu) 113 continue do 118 ll=2,k l = k-ll+1 ir = l-1 i2 = 2**ir i1 = i2/2 i4 = i2+i2 ifd = if-i2 do 117 i=i2,ifd,i4 if (i-i2-izr) 117,114,117 114 if (i-nm) 115,115,118 115 call indxa (i,ir,idxa,na) call indxb (i,ir,iz,nz) call indxb (i-i1,ir-1,im1,nm1) call indxb (i+i1,ir-1,ip1,np1) call prdct (nm1,b(im1),0,dum,0,dum,na,an(idxa),w2,w2,m,am, & bm,cm,wd,ww,wu) do j=1,m w2(j) = y(j,i)+w2(j) end do call prdct (nz,b(iz),nm1,b(im1),np1,b(ip1),0,dum,w2,w2,m, & am,bm,cm,wd,ww,wu) izr = i if (i-nm) 117,119,117 117 continue 118 continue 119 do j=1,m y(j,nm+1) = y(j,nm+1)-cn(nm+1)*w1(j)-an(nm+1)*w2(j) end do call indxb (if/2,k-1,im1,nm1) call indxb (if,k-1,ip,np) if (ncmplx) 121,122,121 121 call cprdct (nm+1,b(ip),nm1,b(im1),0,dum,0,dum,y(1,nm+1), & y(1,nm+1),m,am,bm,cm,w1,w3,ww) go to 123 122 call prdct (nm+1,b(ip),nm1,b(im1),0,dum,0,dum,y(1,nm+1), & y(1,nm+1),m,am,bm,cm,wd,ww,wu) 123 do j=1,m w1(j) = an(1)*y(j,nm+1) w2(j) = cn(nm)*y(j,nm+1) y(j,1) = y(j,1)-w1(j) y(j,nm) = y(j,nm)-w2(j) end do do 126 l=1,kdo ir = l-1 i2 = 2**ir i4 = i2+i2 i1 = i2/2 i = i4 call indxa (i,ir,idxa,na) call indxb (i-i2,ir,im2,nm2) call indxb (i-i2-i1,ir-1,im3,nm3) call indxb (i-i1,ir-1,im1,nm1) call prdct (nm2,b(im2),nm3,b(im3),nm1,b(im1),0,dum,w1,w1,m,am, & bm,cm,wd,ww,wu) call prdct (nm1,b(im1),0,dum,0,dum,na,an(idxa),w1,w1,m,am,bm, & cm,wd,ww,wu) do j=1,m y(j,i) = y(j,i)-w1(j) end do 126 continue izr = nm do 131 l=1,kdo ir = l-1 i2 = 2**ir i1 = i2/2 i3 = i2+i1 i4 = i2+i2 irm1 = ir-1 do 130 i=i4,if,i4 ipi1 = i+i1 ipi2 = i+i2 ipi3 = i+i3 if (ipi2-izr) 127,128,127 127 if (i-izr) 130,131,130 128 call indxc (i,ir,idxc,nc) call indxb (ipi2,ir,ip2,np2) call indxb (ipi1,irm1,ip1,np1) call indxb (ipi3,irm1,ip3,np3) call prdct (np2,b(ip2),np1,b(ip1),np3,b(ip3),0,dum,w2,w2,m, & am,bm,cm,wd,ww,wu) call prdct (np1,b(ip1),0,dum,0,dum,nc,cn(idxc),w2,w2,m,am, & bm,cm,wd,ww,wu) do j=1,m y(j,i) = y(j,i)-w2(j) end do izr = i go to 131 130 continue 131 continue ! ! begin back substitution phase ! 132 do 144 ll=1,k l = k-ll+1 ir = l-1 irm1 = ir-1 i2 = 2**ir i1 = i2/2 i4 = i2+i2 ifd = if-i2 do 143 i=i2,ifd,i4 if (i-nm) 133,133,143 133 imi1 = i-i1 imi2 = i-i2 ipi1 = i+i1 ipi2 = i+i2 call indxa (i,ir,idxa,na) call indxc (i,ir,idxc,nc) call indxb (i,ir,iz,nz) call indxb (imi1,irm1,im1,nm1) call indxb (ipi1,irm1,ip1,np1) if (i-i2) 134,134,136 134 do j=1,m w1(j) = 0. end do go to 137 136 call prdct (nm1,b(im1),0,dum,0,dum,na,an(idxa),y(1,imi2), & w1,m,am,bm,cm,wd,ww,wu) 137 if (ipi2-nm) 140,140,138 138 do j=1,m w2(j) = 0. end do go to 141 140 call prdct (np1,b(ip1),0,dum,0,dum,nc,cn(idxc),y(1,ipi2), & w2,m,am,bm,cm,wd,ww,wu) 141 do j=1,m w1(j) = y(j,i)+w1(j)+w2(j) end do call prdct (nz,b(iz),nm1,b(im1),np1,b(ip1),0,dum,w1,y(1,i), & m,am,bm,cm,wd,ww,wu) 143 continue 144 continue return end subroutine blktri ( iflg, np, n, an, bn, cn, mp, m, am, bm, cm, & idimy, y, ierror, w ) ! !*********************************************************************** ! !! BLKTRI solves a system of linear equations of the form ! ! an(j)*x(i,j-1) + am(i)*x(i-1,j) + (bn(j)+bm(i))*x(i,j) ! ! + cn(j)*x(i,j+1) + cm(i)*x(i+1,j) = y(i,j) ! ! for i = 1,2,...,m and j = 1,2,...,n. ! ! i+1 and i-1 are evaluated modulo m and j+1 and j-1 modulo n, i.e., ! ! x(i,0) = x(i,n), x(i,n+1) = x(i,1), ! x(0,j) = x(m,j), x(m+1,j) = x(1,j). ! ! These equations usually result from the discretization of ! separable elliptic equations. boundary conditions may be ! dirichlet, neumann, or periodic. ! ! ! on input ! ! iflg ! = 0 initialization only. certain quantities that depend on np, ! n, an, bn, and cn are computed and stored in the work ! array w. ! = 1 the quantities that were computed in the initialization are ! used to obtain the solution x(i,j). ! ! note a call with iflg=0 takes approximately one half the time ! time as a call with iflg = 1. however, the ! initialization does not have to be repeated unless np, n, ! an, bn, or cn change. ! ! np ! = 0 if an(1) and cn(n) are not both zero, which corresponds to ! periodic boundary conditions. ! = 1 if an(1) = cn(n) = 0. ! ! n ! the number of unknowns in the j-direction. n must be greater ! than 4. the operation count is proportional to mnlog2(n), hence ! n should be selected less than or equal to m. ! ! an,bn,cn ! one-dimensional arrays of length n that specify the coefficients ! in the linear equations given above. ! ! mp ! = 0 if am(1) and cm(m) are not both zero, which corresponds to ! periodic boundary conditions. ! = 1 if am(1) = cm(m) = 0 . ! ! m ! the number of unknowns in the i-direction. m must be greater ! than 4. ! ! am,bm,cm ! one-dimensional arrays of length m that specify the coefficients ! in the linear equations given above. ! ! idimy ! the row (or first) dimension of the two-dimensional array y as ! it appears in the program calling blktri. this parameter is ! used to specify the variable dimension of y. idimy must be at ! least m. ! ! y ! a two-dimensional array that specifies the values of the right ! side of the linear system of equations given above. y must be ! dimensioned at least m*n. ! ! w ! a one-dimensional array that must be provided by the user for ! work space. ! if np=1 define k=int(log2(n))+1 and set l=2**(k+1) then ! w must have dimension (k-2)*l+k+5+max(2n,6m) ! ! if np=0 define k=int(log2(n-1))+1 and set l=2**(k+1) then ! w must have dimension (k-2)*l+k+5+2n+max(2n,6m) ! ! **important** for purposes of checking, the required dimension ! of w is computed by blktri and stored in w(1) ! in floating point format. ! ! on output ! ! y ! contains the solution x. ! ! ierror ! an error flag that indicates invalid input parameters. except ! for number zero, a solution is not attempted. ! ! = 0 no error. ! = 1 m is less than 5 ! = 2 n is less than 5 ! = 3 idimy is less than m. ! = 4 blktri failed while computing results that depend on the ! coefficient arrays an, bn, cn. check these arrays. ! = 5 an(j)*cn(j-1) is less than 0 for some j. possible reasons ! for this condition are ! 1. the arrays an and cn are not correct ! 2. too large a grid spacing was used in the discretization ! of the elliptic equation ! 3. the linear equations resulted from a partial ! differential equation which was not elliptic ! ! w ! contains intermediate values that must not be destroyed if ! blktri will be called again with iflg=1. w(1) contains the ! number of locations required by w in floating point format. ! ! * * * * * * * program specifications * * * * * * * * * * * * ! ! dimension of an(n),bn(n),cn(n),am(m),bm(m),cm(m),y(idimy,n) ! arguments w(see argument list) ! ! latest june 1979 ! revision ! ! required blktri,blktri,prod,prodp,cprod,cprodp,compb,indxa, ! subprograms indxb,indxc,ppadd,psgf,ppsgf,ppspf,bsrh,tevls, ! epmach,store ! ! special the algorithm may fail if abs(bm(i)+bn(j)) is less ! conditions than abs(am(i))+abs(an(j))+abs(cm(i))+abs(cn(j)) ! for some i and j. the algorithm will also fail if ! an(j)*cn(j-1) is less than zero for some j ! see the discription of the output parameter ierror. ! ! common cblkt,value ! blocks ! ! specialist paul swarztrauber ! ! history version 1 september 1973 ! version 2 april 1976 ! version 3 june 1979 ! ! algorithm generalized cyclic reduction (see reference below) ! ! references swarztrauber,p. and r. sweet, 'efficient fortran ! subprograms for the solution of elliptic equations' ! ncar tn/ia-109, july, 1975, 138 pp. ! ! swarztrauber p. n.,a direct method for the discrete ! solution of separable elliptic equations, ! SIAM Journal on Numerical Analysis,11(1974) pp. 1136-1150. ! integer idimy integer m integer n ! real am(m) real an(n) real bm(m) real bn(n) real cm(m) real cn(n) real w(*) real y(idimy,n) ! external cprod external cprodp external prod external prodp ! common /cblkt/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik ! nm = n ierror = 0 ! ! Test input. ! if ( m .lt. 5 ) then ierror = 1 write(*,*)' ' write(*,*)'BLKTRI - Fatal error!' write(*,*)' M must be at least 5.' write(*,*)' Your input value is M = ', m return end if if ( n .lt. 5 ) then ierror = 2 write(*,*)' ' write(*,*)'BLKTRI - Fatal error!' write(*,*)' N must be at least 5.' write(*,*)' Your input value is N = ', n return end if if ( idimy .lt. m ) then ierror = 3 write(*,*)' ' write(*,*)'BLKTRI - Fatal error!' write(*,*)' IDIMY must be at least M.' write(*,*)' Your input values are IDIMY = ', idimy write(*,*)' M = ', m return end if ! ! Set up array indices. ! nh = n npp = np if ( npp .ne. 0 ) then nh = nh+1 end if ik = 2 k = 1 10 continue ik = ik+ik k = k+1 if ( nh .gt. ik ) go to 10 nl = ik ik = ik+ik nl = nl-1 iwah = (k-2)*ik+k+6 ! ! Divide W into working subarrays. ! if ( npp .ne. 0 ) then iw1 = iwah iwbh = iw1+nm w(1) = float(iw1-1+max(2*nm,6*m) ) else iwbh = iwah+nm+nm iw1 = iwbh w(1) = float(iw1-1+max(2*nm,6*m)) nm = nm-1 end if ! ! Compute the roots of the B polynomials. ! iw2 = iw1+m iw3 = iw2+m iwd = iw3+m iww = iwd+m iwu = iww+m if ( iflg .eq. 0 ) then call compb ( nl, ierror, an, bn, cn, w(2), w(iwah), w(iwbh) ) return end if ! ! Solve the linear system. ! if ( mp .ne. 0 ) then call blktr1 ( nl, an, bn, cn, m, am, bm, cm, idimy, y, w(2), & w(iw1), w(iw2), w(iw3), w(iwd), w(iww), w(iwu), prod, & cprod ) else call blktr1 ( nl, an, bn, cn, m, am, bm, cm, idimy, y, w(2), & w(iw1), w(iw2), w(iw3), w(iwd), w(iww), w(iwu), prodp, & cprodp ) end if return end function bsrh (xll,xrr,iz,c,a,bh,f,sgn) ! !*********************************************************************** ! dimension a(1) ,c(1) ,bh(1) ! external f ! common /cblkt/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik ! xl = xll xr = xrr dx = .5*abs(xr-xl) 101 continue x = .5*(xl+xr) if (sgn*f(x,iz,c,a,bh)) 103,105,102 102 xr = x go to 104 103 xl = x 104 dx = .5*dx if (dx.gt.cnv) go to 101 105 bsrh = .5*(xl+xr) return end subroutine cblkt1 (n,an,bn,cn,m,am,bm,cm,idimy,y,b,w1,w2,w3,wd, & ww,wu,prdct,cprdct) ! !*********************************************************************** ! ! cblkt1 solves the linear system c c b contains the roots of all the b polynomials c w1,w2,w3,wd,ww,wu are all working arrays c prdct is either procp or proc depending on whether the boundary c conditions in the m direction are periodic or not c cprdct is either cprocp or cproc which are called if some of the zeros c of the b polynomials are complex. c c dimension an(1) ,bn(1) ,cn(1) ,am(1) , & bm(1) ,cm(1) ,b(1) ,w1(1) , & w2(1) ,w3(1) ,wd(1) ,ww(1) , & wu(1) ,y(idimy,1) c external cprdct external prdct c common /ccblk/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik complex am ,bm ,cm ,y , & w1 ,w2 ,w3 ,wd , & ww ,wu c c begin reduction phase c kdo = k-1 do 109 l=1,kdo ir = l-1 i2 = 2**ir i1 = i2/2 i3 = i2+i1 i4 = i2+i2 irm1 = ir-1 call inxcb (i2,ir,im2,nm2) call inxcb (i1,irm1,im3,nm3) call inxcb (i3,irm1,im1,nm1) call prdct (nm2,b(im2),nm3,b(im3),nm1,b(im1),0,dum,y(1,i2),w3, & m,am,bm,cm,wd,ww,wu) if = 2**k do 108 i=i4,if,i4 if (i-nm) 101,101,108 101 ipi1 = i+i1 ipi2 = i+i2 ipi3 = i+i3 call inxcc (i,ir,idxc,nc) if (i-if) 102,108,108 102 call inxca (i,ir,idxa,na) call inxcb (i-i1,irm1,im1,nm1) call inxcb (ipi2,ir,ip2,np2) call inxcb (ipi1,irm1,ip1,np1) call inxcb (ipi3,irm1,ip3,np3) call prdct (nm1,b(im1),0,dum,0,dum,na,an(idxa),w3,w1,m,am, & bm,cm,wd,ww,wu) if (ipi2-nm) 105,105,103 103 do j=1,m w3(j) = (0.,0.) w2(j) = (0.,0.) end do go to 106 105 call prdct (np2,b(ip2),np1,b(ip1),np3,b(ip3),0,dum, & y(1,ipi2),w3,m,am,bm,cm,wd,ww,wu) call prdct (np1,b(ip1),0,dum,0,dum,nc,cn(idxc),w3,w2,m,am, & bm,cm,wd,ww,wu) 106 do 107 j=1,m y(j,i) = w1(j)+w2(j)+y(j,i) 107 continue 108 continue 109 continue if (npp) 132,110,132 c c the periodic case is treated using the capacitance matrix method c 110 if = 2**k i = if/2 i1 = i/2 call inxcb (i-i1,k-2,im1,nm1) call inxcb (i+i1,k-2,ip1,np1) call inxcb (i,k-1,iz,nz) call prdct (nz,b(iz),nm1,b(im1),np1,b(ip1),0,dum,y(1,i),w1,m,am, & bm,cm,wd,ww,wu) izr = i do 111 j=1,m w2(j) = w1(j) 111 continue do 113 ll=2,k l = k-ll+1 ir = l-1 i2 = 2**ir i1 = i2/2 i = i2 call inxcc (i,ir,idxc,nc) call inxcb (i,ir,iz,nz) call inxcb (i-i1,ir-1,im1,nm1) call inxcb (i+i1,ir-1,ip1,np1) call prdct (np1,b(ip1),0,dum,0,dum,nc,cn(idxc),w1,w1,m,am,bm, & cm,wd,ww,wu) do j=1,m w1(j) = y(j,i)+w1(j) end do call prdct (nz,b(iz),nm1,b(im1),np1,b(ip1),0,dum,w1,w1,m,am, & bm,cm,wd,ww,wu) 113 continue do 118 ll=2,k l = k-ll+1 ir = l-1 i2 = 2**ir i1 = i2/2 i4 = i2+i2 ifd = if-i2 do 117 i=i2,ifd,i4 if (i-i2-izr) 117,114,117 114 if (i-nm) 115,115,118 115 call inxca (i,ir,idxa,na) call inxcb (i,ir,iz,nz) call inxcb (i-i1,ir-1,im1,nm1) call inxcb (i+i1,ir-1,ip1,np1) call prdct (nm1,b(im1),0,dum,0,dum,na,an(idxa),w2,w2,m,am, & bm,cm,wd,ww,wu) do 116 j=1,m w2(j) = y(j,i)+w2(j) 116 continue call prdct (nz,b(iz),nm1,b(im1),np1,b(ip1),0,dum,w2,w2,m, & am,bm,cm,wd,ww,wu) izr = i if (i-nm) 117,119,117 117 continue 118 continue 119 do 120 j=1,m y(j,nm+1) = y(j,nm+1)-cn(nm+1)*w1(j)-an(nm+1)*w2(j) 120 continue call inxcb (if/2,k-1,im1,nm1) call inxcb (if,k-1,ip,np) if (ncmplx) 121,122,121 121 call cprdct (nm+1,b(ip),nm1,b(im1),0,dum,0,dum,y(1,nm+1), & y(1,nm+1),m,am,bm,cm,w1,w3,ww) go to 123 122 call prdct (nm+1,b(ip),nm1,b(im1),0,dum,0,dum,y(1,nm+1), & y(1,nm+1),m,am,bm,cm,wd,ww,wu) 123 do 124 j=1,m w1(j) = an(1)*y(j,nm+1) w2(j) = cn(nm)*y(j,nm+1) y(j,1) = y(j,1)-w1(j) y(j,nm) = y(j,nm)-w2(j) 124 continue do 126 l=1,kdo ir = l-1 i2 = 2**ir i4 = i2+i2 i1 = i2/2 i = i4 call inxca (i,ir,idxa,na) call inxcb (i-i2,ir,im2,nm2) call inxcb (i-i2-i1,ir-1,im3,nm3) call inxcb (i-i1,ir-1,im1,nm1) call prdct (nm2,b(im2),nm3,b(im3),nm1,b(im1),0,dum,w1,w1,m,am, & bm,cm,wd,ww,wu) call prdct (nm1,b(im1),0,dum,0,dum,na,an(idxa),w1,w1,m,am,bm, & cm,wd,ww,wu) do 125 j=1,m y(j,i) = y(j,i)-w1(j) 125 continue 126 continue c izr = nm do 131 l=1,kdo ir = l-1 i2 = 2**ir i1 = i2/2 i3 = i2+i1 i4 = i2+i2 irm1 = ir-1 do 130 i=i4,if,i4 ipi1 = i+i1 ipi2 = i+i2 ipi3 = i+i3 if (ipi2-izr) 127,128,127 127 if (i-izr) 130,131,130 128 call inxcc (i,ir,idxc,nc) call inxcb (ipi2,ir,ip2,np2) call inxcb (ipi1,irm1,ip1,np1) call inxcb (ipi3,irm1,ip3,np3) call prdct (np2,b(ip2),np1,b(ip1),np3,b(ip3),0,dum,w2,w2,m, & am,bm,cm,wd,ww,wu) call prdct (np1,b(ip1),0,dum,0,dum,nc,cn(idxc),w2,w2,m,am, & bm,cm,wd,ww,wu) do 129 j=1,m y(j,i) = y(j,i)-w2(j) 129 continue izr = i go to 131 130 continue 131 continue c c begin back substitution phase c 132 do 144 ll=1,k l = k-ll+1 ir = l-1 irm1 = ir-1 i2 = 2**ir i1 = i2/2 i4 = i2+i2 ifd = if-i2 do 143 i=i2,ifd,i4 if (i-nm) 133,133,143 133 imi1 = i-i1 imi2 = i-i2 ipi1 = i+i1 ipi2 = i+i2 call inxca (i,ir,idxa,na) call inxcc (i,ir,idxc,nc) call inxcb (i,ir,iz,nz) call inxcb (imi1,irm1,im1,nm1) call inxcb (ipi1,irm1,ip1,np1) if (i-i2) 134,134,136 134 do 135 j=1,m w1(j) = (0.,0.) 135 continue go to 137 136 call prdct (nm1,b(im1),0,dum,0,dum,na,an(idxa),y(1,imi2), & w1,m,am,bm,cm,wd,ww,wu) 137 if (ipi2-nm) 140,140,138 138 do 139 j=1,m w2(j) = (0.,0.) 139 continue go to 141 140 call prdct (np1,b(ip1),0,dum,0,dum,nc,cn(idxc),y(1,ipi2), & w2,m,am,bm,cm,wd,ww,wu) 141 do 142 j=1,m w1(j) = y(j,i)+w1(j)+w2(j) 142 continue call prdct (nz,b(iz),nm1,b(im1),np1,b(ip1),0,dum,w1,y(1,i), & m,am,bm,cm,wd,ww,wu) 143 continue 144 continue return end subroutine cblktr (iflg,np,n,an,bn,cn,mp,m,am,bm,cm,idimy,y, & ierror,w) c c cblktr is a complex version of blktri c both routines solve a system of linear equations of the form c c an(j)*x(i,j-1) + am(i)*x(i-1,j) + (bn(j)+bm(i))*x(i,j) c c + cn(j)*x(i,j+1) + cm(i)*x(i+1,j) = y(i,j) c c for i = 1,2,...,m and j = 1,2,...,n. c c i+1 and i-1 are evaluated modulo m and j+1 and j-1 modulo n, i.e., c c x(i,0) = x(i,n), x(i,n+1) = x(i,1), c x(0,j) = x(m,j), x(m+1,j) = x(1,j). c c these equations usually result from the discretization of c separable elliptic equations. boundary conditions may be c dirichlet, neumann, or periodic. c c c * * * * * * * * * * on input * * * * * * * * * * c c iflg c = 0 initialization only. certain quantities that depend on np, c n, an, bn, and cn are computed and stored in the work c array w. c = 1 the quantities that were computed in the initialization are c used to obtain the solution x(i,j). c c note a call with iflg=0 takes approximately one half the time c time as a call with iflg = 1 . however, the c initialization does not have to be repeated unless np, n, c an, bn, or cn change. c c np c = 0 if an(1) and cn(n) are not zero, which corresponds to c periodic bounary conditions. c = 1 if an(1) and cn(n) are zero. c c n c the number of unknowns in the j-direction. n must be greater c than 4. the operation count is proportional to mnlog2(n), hence c n should be selected less than or equal to m. c c an,bn,cn c real one-dimensional arrays of length n that specify the c coefficients in the linear equations given above. c c mp c = 0 if am(1) and cm(m) are not zero, which corresponds to c periodic boundary conditions. c = 1 if am(1) = cm(m) = 0 . c c m c the number of unknowns in the i-direction. m must be greater c than 4. c c am,bm,cm c complex one-dimensional arrays of length m that specify the c coefficients in the linear equations given above. c c idimy c the row (or first) dimension of the two-dimensional array y as c it appears in the program calling blktri. this parameter is c used to specify the variable dimension of y. idimy must be at c least m. c c y c a complex two-dimensional array that specifies the values of the c right side of the linear system of equations given above. y must c be dimensioned y(idimy,n) with idimy .ge. m. c c w c a one-dimensional array that must be provided by the user for c work space. c if np=1 define k=int(log2(n))+1 and set l=2**(k+1) then c w must have dimension (k-2)*l+k+5+max(2n,12m) c c if np=0 define k=int(log2(n-1))+1 and set l=2**(k+1) then c w must have dimension (k-2)*l+k+5+2n+max(2n,12m) c c **important** for purposes of checking, the required dimension c of w is computed by blktri and stored in w(1) c in floating point format. c c * * * * * * * * * * on output * * * * * * * * * * c c y c contains the solution x. c c ierror c an error flag that indicates invalid input parameters. except c for number zero, a solution is not attempted. c c = 0 no error. c = 1 m is less than 5 c = 2 n is less than 5 c = 3 idimy is less than m. c = 4 blktri failed while computing results that depend on the c coefficient arrays an, bn, cn. check these arrays. c = 5 an(j)*cn(j-1) is less than 0 for some j. possible reasons c for this condition are c 1. the arrays an and cn are not correct c 2. too large a grid spacing was used in the discretization c of the elliptic equation c 3. the linear equations resulted from a partial c differential equation which was not elliptic c c w c contains intermediate values that must not be destroyed if c cblktr will be called again with iflg=1. w(1) contains the c number of locations required by w in floating point format. c c * * * * * * * program specifications * * * * * * * * * * * * c c dimension of an(n),bn(n),cn(n),am(m),bm(m),cm(m),y(idimy,n) c arguments w(see argument list) c c latest june 1979 c revision c c required cblktr,cblkt1,proc,procp,cproc,cprocp,ccmpb,inxca, c subprograms inxcb,inxcc,cpadd,pgsf,ppgsf,pppsf,bcrh,tevlc, c epmach,store c c special the algorithm may fail if abs(bm(i)+bn(j)) is less c conditions than abs(am(i))+abs(an(j))+abs(cm(i))+abs(cn(j)) c for some i and j. the algorithm will also fail if c an(j)*cn(j-1) is less than zero for some j c see the discription of the output parameter ierror. c c common ccblk c blocks c c i/o none c c precision single c c specialist paul swarztrauber c c language fortran c c history cblktr is a complex version of blktri (version 3) c c algorithm generalized cyclic reduction (see reference below) c c space c required control data 7600 c c portability american national standards institute fortran. c the approximate machine accuracy is computed in c function epmach c c required none c resident c routines c c references swarztrauber,p. and r. sweet, 'efficient fortran c subprograms for the solution of elliptic equations' c ncar tn/ia-109, july, 1975, 138 pp. c c swarztrauber p. n.,a direct method for the discrete c solution of separable elliptic equations, c SIAM Journal on Numerical Analysis,11(1974) pp. 1136-1150. c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c dimension an(1) ,bn(1) ,cn(1) ,am(1) , & bm(1) ,cm(1) ,y(idimy,1) ,w(1) external proc ,procp ,cproc ,cprocp common /ccblk/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik complex am ,bm ,cm ,y c c test m and n for the proper form c nm = n m2 = m+m ierror = 0 if (m-5) 101,102,102 101 ierror = 1 go to 119 102 if (nm-5) 103,104,104 103 ierror = 2 go to 119 104 if (idimy-m) 105,106,106 105 ierror = 3 go to 119 106 nh = n npp = np if (npp) 107,108,107 107 nh = nh+1 108 ik = 2 k = 1 109 ik = ik+ik k = k+1 if (nh-ik) 110,110,109 110 nl = ik ik = ik+ik nl = nl-1 iwah = (k-2)*ik+k+6 if (npp) 111,112,111 c c divide w into working sub arrays c 111 iw1 = iwah iwbh = iw1+nm w(1) = float(iw1-1+max0(2*nm,12*m)) go to 113 112 iwbh = iwah+nm+nm iw1 = iwbh w(1) = float(iw1-1+max0(2*nm,12*m)) nm = nm-1 c c ccmpb computes the roots of the b polynomials c 113 continue ! 113 if (ierror) 119,114,119 114 iw2 = iw1+m2 iw3 = iw2+m2 iwd = iw3+m2 iww = iwd+m2 iwu = iww+m2 if (iflg) 116,115,116 115 call ccmpb (nl,ierror,an,bn,cn,w(2),w(iwah),w(iwbh)) go to 119 116 if (mp) 117,118,117 c c cblkt1 solves the linear system c 117 call cblkt1 (nl,an,bn,cn,m,am,bm,cm,idimy,y,w(2),w(iw1),w(iw2), & w(iw3),w(iwd),w(iww),w(iwu),proc,cproc) go to 119 118 call cblkt1 (nl,an,bn,cn,m,am,bm,cm,idimy,y,w(2),w(iw1),w(iw2), & w(iw3),w(iwd),w(iww),w(iwu),procp,cprocp) 119 continue return end subroutine ccmpb (n,ierror,an,bn,cn,b,ah,bh) c c ccmpb computes the roots of the b polynomials using routine c tevlc which is a modification the eispack program tqlrat. c ierror is set to 4 if either tevlc fails or if a(j+1)*c(j) is c less than zero for some j. ah,bh are temporary work arrays. c dimension an(1) ,bn(1) ,cn(1) ,b(1) , & ah(1) ,bh(1) common /ccblk/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik eps = epmach() bnorm = abs(bn(1)) do 102 j=2,nm bnorm = amax1(bnorm,abs(bn(j))) arg = an(j)*cn(j-1) if (arg) 119,101,101 101 b(j) = sign(sqrt(arg),an(j)) 102 continue cnv = eps*bnorm if = 2**k kdo = k-1 do 108 l=1,kdo ir = l-1 i2 = 2**ir i4 = i2+i2 ipl = i4-1 ifd = if-i4 do 107 i=i4,ifd,i4 call inxcb (i,l,ib,nb) if (nb) 108,108,103 103 js = i-ipl jf = js+nb-1 ls = 0 do 104 j=js,jf ls = ls+1 bh(ls) = bn(j) ah(ls) = b(j) 104 continue call tevlc (nb,bh,ah,ierror) if (ierror) 118,105,118 105 lh = ib-1 do 106 j=1,nb lh = lh+1 b(lh) = -bh(j) 106 continue 107 continue 108 continue do 109 j=1,nm b(j) = -bn(j) 109 continue if (npp) 117,110,117 110 nmp = nm+1 nb = nm+nmp do 112 j=1,nb l1 = mod(j-1,nmp)+1 l2 = mod(j+nm-1,nmp)+1 arg = an(l1)*cn(l2) if (arg) 119,111,111 111 bh(j) = sign(sqrt(arg),-an(l1)) ah(j) = -bn(l1) 112 continue call tevlc (nb,ah,bh,ierror) if (ierror) 118,113,118 113 call inxcb (if,k-1,j2,lh) call inxcb (if/2,k-1,j1,lh) j2 = j2+1 lh = j2 n2m2 = j2+nm+nm-2 114 d1 = abs(b(j1)-b(j2-1)) d2 = abs(b(j1)-b(j2)) d3 = abs(b(j1)-b(j2+1)) if ((d2 .lt. d1) .and. (d2 .lt. d3)) go to 115 b(lh) = b(j2) j2 = j2+1 lh = lh+1 if (j2-n2m2) 114,114,116 115 j2 = j2+1 j1 = j1+1 if (j2-n2m2) 114,114,116 116 b(lh) = b(n2m2+1) call inxcb (if,k-1,j1,j2) j2 = j1+nmp+nmp call cpadd (nm+1,ierror,an,cn,b(j1),b(j1),b(j2)) 117 return 118 ierror = 4 return 119 ierror = 5 return end subroutine chkpr4(iorder,a,b,m,mbdcnd,c,d,n,nbdcnd,cofx,idmn, &ierror) external cofx c c this program checks the input parameters for errors c c c c check definition of solution region c ierror = 1 if (a.ge.b .or. c.ge.d) return c c check boundary switches c ierror = 2 if (mbdcnd.lt.0 .or. mbdcnd.gt.4) return ierror = 3 if (nbdcnd.lt.0 .or. nbdcnd.gt.4) return c c check first dimension in calling routine c ierror = 5 if (idmn .lt. 7) return c c check m c ierror = 6 if (m.gt.(idmn-1) .or. m.lt.6) return c c check n c ierror = 7 if (n .lt. 5) return c c check iorder c ierror = 8 if (iorder.ne.2 .and. iorder.ne.4) return c c check intl c c c check that equation is elliptic c dlx = (b-a)/float(m) do 30 i=2,m xi = a+float(i-1)*dlx call cofx (xi,ai,bi,ci) if (ai.gt.0.0) go to 10 ierror=10 return 10 continue 30 continue c c no error found c ierror = 0 return end subroutine chkprm (intl,iorder,a,b,m,mbdcnd,c,d,n,nbdcnd,cofx, & cofy,idmn,ierror) c c this program checks the input parameters for errors c external cofx ,cofy c c check definition of solution region c ierror = 1 if (a.ge.b .or. c.ge.d) return c c check boundary switches c ierror = 2 if (mbdcnd.lt.0 .or. mbdcnd.gt.4) return ierror = 3 if (nbdcnd.lt.0 .or. nbdcnd.gt.4) return c c check first dimension in calling routine c ierror = 5 if (idmn .lt. 7) return c c check m c ierror = 6 if (m.gt.(idmn-1) .or. m.lt.6) return c c check n c ierror = 7 if (n .lt. 5) return c c check iorder c ierror = 8 if (iorder.ne.2 .and. iorder.ne.4) return c c check intl c ierror = 9 if (intl.ne.0 .and. intl.ne.1) return c c check that equation is elliptic c dlx = (b-a)/float(m) dly = (d-c)/float(n) do 30 i=2,m xi = a+float(i-1)*dlx call cofx (xi,ai,bi,ci) do 20 j=2,n yj = c+float(j-1)*dly call cofy (yj,dj,ej,fj) if (ai*dj .gt. 0.0) go to 10 ierror = 10 return 10 continue 20 continue 30 continue c c no error found c ierror = 0 return end subroutine chksn4(mbdcnd,nbdcnd,alpha,beta,cofx,singlr) c c this subroutine checks if the pde sepx4 c must solve is a singular operator c common /spl4/ kswx ,kswy ,k ,l , & ait ,bit ,cit ,dit , & mit ,nit ,is ,ms , & js ,ns ,dlx ,dly , & tdlx3 ,tdly3 ,dlx4 ,dly4 logical singlr external cofx singlr = .false. c c check if the boundary conditions are c entirely periodic and/or mixed c if ((mbdcnd.ne.0 .and. mbdcnd.ne.3) .or. & (nbdcnd.ne.0 .and. nbdcnd.ne.3)) return c c check that mixed conditions are pure neuman c if (mbdcnd .ne. 3) go to 10 if (alpha.ne.0.0 .or. beta.ne.0.0) return 10 continue c c check that non-derivative coefficient functions c are zero c do 30 i=is,ms xi = ait+float(i-1)*dlx call cofx (xi,ai,bi,ci) if (ci .ne. 0.0) return 30 continue c c the operator must be singular if this point is reached c singlr = .true. return end subroutine chksng (mbdcnd,nbdcnd,alpha,beta,gama,xnu,cofx,cofy, & singlr) c c this subroutine checks if the pde sepeli c must solve is a singular operator c external cofx external cofy c common /splp/ kswx ,kswy ,k ,l , & ait ,bit ,cit ,dit , & mit ,nit ,is ,ms , & js ,ns ,dlx ,dly , & tdlx3 ,tdly3 ,dlx4 ,dly4 logical singlr singlr = .false. c c check if the boundary conditions are c entirely periodic and/or mixed c if ((mbdcnd.ne.0 .and. mbdcnd.ne.3) .or. & (nbdcnd.ne.0 .and. nbdcnd.ne.3)) return c c check that mixed conditions are pure neuman c if (mbdcnd .ne. 3) go to 10 if (alpha.ne.0.0 .or. beta.ne.0.0) return 10 if (nbdcnd .ne. 3) go to 20 if (gama.ne.0.0 .or. xnu.ne.0.0) return 20 continue c c check that non-derivative coefficient functions c are zero c do 30 i=is,ms xi = ait+float(i-1)*dlx call cofx (xi,ai,bi,ci) if (ci .ne. 0.0) return 30 continue do 40 j=js,ns yj = cit+float(j-1)*dly call cofy (yj,dj,ej,fj) if (fj .ne. 0.0) return 40 continue c c the operator must be singular if this point is reached c singlr = .true. return end subroutine cmgnbn (nperod,n,mperod,m,a,b,c,idimy,y,ierror,w) c c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c * * c * f i s h p a k * c * * c * * c * a package of fortran subprograms for the solution of * c * * c * separable elliptic partial differential equations * c * * c * (version 3.1 , october 1980) * c * * c * by * c * * c * john adams, paul swarztrauber and roland sweet * c * * c * of * c * * c * the national center for atmospheric research * c * * c * boulder, colorado (80307) u.s.a. * c * * c * which is sponsored by * c * * c * the national science foundation * c * * c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c c c * * * * * * * * * purpose * * * * * * * * * * * * * * * * * * c c c subroutine cmgnbn solves the complex linear system of equations c c a(i)*x(i-1,j) + b(i)*x(i,j) + c(i)*x(i+1,j) c c + x(i,j-1) - 2.*x(i,j) + x(i,j+1) = y(i,j) c c for i = 1,2,...,m and j = 1,2,...,n. c c the indices i+1 and i-1 are evaluated modulo m, i.e., c x(0,j) = x(m,j) and x(m+1,j) = x(1,j), and x(i,0) may be equal to c 0, x(i,2), or x(i,n) and x(i,n+1) may be equal to 0, x(i,n-1), or c x(i,1) depending on an input parameter. c c c * * * * * * * * parameter description * * * * * * * * * * c c * * * * * * on input * * * * * * c c nperod c indicates the values that x(i,0) and x(i,n+1) are assumed to c have. c c = 0 if x(i,0) = x(i,n) and x(i,n+1) = x(i,1). c = 1 if x(i,0) = x(i,n+1) = 0 . c = 2 if x(i,0) = 0 and x(i,n+1) = x(i,n-1). c = 3 if x(i,0) = x(i,2) and x(i,n+1) = x(i,n-1). c = 4 if x(i,0) = x(i,2) and x(i,n+1) = 0. c c n c the number of unknowns in the j-direction. n must be greater c than 2. c c mperod c = 0 if a(1) and c(m) are not zero c = 1 if a(1) = c(m) = 0 c c m c the number of unknowns in the i-direction. n must be greater c than 2. c c a,b,c c one-dimensional complex arrays of length m that specify the c coefficients in the linear equations given above. if mperod = 0 c the array elements must not depend upon the index i, but must be c constant. specifically, the subroutine checks the following c condition c c a(i) = c(1) c c(i) = c(1) c b(i) = b(1) c c for i=1,2,...,m. c c idimy c the row (or first) dimension of the two-dimensional array y as c it appears in the program calling cmgnbn. this parameter is c used to specify the variable dimension of y. idimy must be at c least m. c c y c a two-dimensional complex array that specifies the values of the c right side of the linear system of equations given above. y c must be dimensioned at least m*n. c c w c a one-dimensional complex array that must be provided by the c user for work space. w may require up to 4*n + c (10 + int(log2(n)))*m locations. the actual number of locations c used is computed by cmgnbn and is returned in location w(1). c c c * * * * * * on output * * * * * * c c y c contains the solution x. c c ierror c an error flag which indicates invalid input parameters except c for number zero, a solution is not attempted. c c = 0 no error. c = 1 m .le. 2 . c = 2 n .le. 2 c = 3 idimy .lt. m c = 4 nperod .lt. 0 or nperod .gt. 4 c = 5 mperod .lt. 0 or mperod .gt. 1 c = 6 a(i) .ne. c(1) or c(i) .ne. c(1) or b(i) .ne. b(1) for c some i=1,2,...,m. c = 7 a(1) .ne. 0 or c(m) .ne. 0 and mperod = 1 c c w c w(1) contains the required length of w. c c * * * * * * * program specifications * * * * * * * * * * * * c c dimension of a(m),b(m),c(m),y(idimy,n),w(see parameter list) c arguments c c latest june 1979 c revision c c subprograms cmgnbn,cmposd,cmposn,cmposp,cmpcsg,cmpmrg, c required cmptrx,cmptr3,pimach c c special none c conditions c c common none c blocks c c i/o none c c precision single c c specialist roland sweet c c language fortran c c history written by roland sweet at ncar in june, 1977 c c algorithm the linear system is solved by a cyclic reduction c algorithm described in the reference. c c space 4944(decimal) = 11520(octal) locations on the ncar c required control data 7600 c c timing and the execution time t on the ncar control data c accuracy 7600 for subroutine cmgnbn is roughly proportional c to m*n*log2(n), but also depends on the input c parameter nperod. some typical values are listed c in the table below. c to measure the accuracy of the algorithm a c uniform random number generator was used to create c a solution array x for the system given in the c 'purpose' with c c a(i) = c(i) = -0.5*b(i) = 1, i=1,2,...,m c c and, when mperod = 1 c c a(1) = c(m) = 0 c a(m) = c(1) = 2. c c the solution x was substituted into the given sys- c tem and a right side y was computed. using this c array y subroutine cmgnbn was called to produce an c approximate solution z. then the relative error, c defined as c c e = max(cabs(z(i,j)-x(i,j)))/max(cabs(x(i,j))) c c where the two maxima are taken over all i=1,2,...,m c and j=1,2,...,n, was computed. the value of e is c given in the table below for some typical values of c m and n. c c c m (=n) mperod nperod t(msecs) e c ------ ------ ------ -------- ------ c c 31 0 0 77 1.e-12 c 31 1 1 45 4.e-13 c 31 1 3 91 2.e-12 c 32 0 0 59 7.e-14 c 32 1 1 65 5.e-13 c 32 1 3 97 2.e-13 c 33 0 0 80 6.e-13 c 33 1 1 67 5.e-13 c 33 1 3 76 3.e-12 c 63 0 0 350 5.e-12 c 63 1 1 215 6.e-13 c 63 1 3 412 1.e-11 c 64 0 0 264 1.e-13 c 64 1 1 287 3.e-12 c 64 1 3 421 3.e-13 c 65 0 0 338 2.e-12 c 65 1 1 292 5.e-13 c 65 1 3 329 1.e-11 c c portability american national standards institue fortran. c all machine dependent constants are located in the c the machine dependent constant pi is defined in c function pimach. c c required cos c resident c routines c c reference sweet, r., 'a cyclic reduction algorithm for c solving block tridiagonal systems of arbitrary c dimensions,' siam j. on numer. anal., c 14(sept., 1977), pp. 706-720. c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c c complex a ,b ,c ,y , & w ,a1 dimension y(idimy,1) dimension w(1) ,b(1) ,a(1) ,c(1) ierror = 0 if (m .le. 2) ierror = 1 if (n .le. 2) ierror = 2 if (idimy .lt. m) ierror = 3 if (nperod.lt.0 .or. nperod.gt.4) ierror = 4 if (mperod.lt.0 .or. mperod.gt.1) ierror = 5 if (mperod .eq. 1) go to 102 do 101 i=2,m if (cabs(a(i)-c(1)) .ne. 0.) go to 103 if (cabs(c(i)-c(1)) .ne. 0.) go to 103 if (cabs(b(i)-b(1)) .ne. 0.) go to 103 101 continue go to 104 102 if (cabs(a(1)).ne.0. .and. cabs(c(m)).ne.0.) ierror = 7 go to 104 103 ierror = 6 104 if (ierror .ne. 0) return iwba = m+1 iwbb = iwba+m iwbc = iwbb+m iwb2 = iwbc+m iwb3 = iwb2+m iww1 = iwb3+m iww2 = iww1+m iww3 = iww2+m iwd = iww3+m iwtcos = iwd+m iwp = iwtcos+4*n do 106 i=1,m k = iwba+i-1 w(k) = -a(i) k = iwbc+i-1 w(k) = -c(i) k = iwbb+i-1 w(k) = 2.-b(i) do 105 j=1,n y(i,j) = -y(i,j) 105 continue 106 continue mp = mperod+1 np = nperod+1 go to (114,107),mp 107 go to (108,109,110,111,123),np 108 call cmposp (m,n,w(iwba),w(iwbb),w(iwbc),y,idimy,w,w(iwb2), & w(iwb3),w(iww1),w(iww2),w(iww3),w(iwd),w(iwtcos), & w(iwp)) go to 112 109 call cmposd (m,n,1,w(iwba),w(iwbb),w(iwbc),y,idimy,w,w(iww1), & w(iwd),w(iwtcos),w(iwp)) go to 112 110 call cmposn (m,n,1,2,w(iwba),w(iwbb),w(iwbc),y,idimy,w,w(iwb2), & w(iwb3),w(iww1),w(iww2),w(iww3),w(iwd),w(iwtcos), & w(iwp)) go to 112 111 call cmposn (m,n,1,1,w(iwba),w(iwbb),w(iwbc),y,idimy,w,w(iwb2), & w(iwb3),w(iww1),w(iww2),w(iww3),w(iwd),w(iwtcos), & w(iwp)) 112 ipstor = real(w(iww1)) irev = 2 if (nperod .eq. 4) go to 124 113 go to (127,133),mp 114 continue c c reorder unknowns when mp =0 c mh = (m+1)/2 mhm1 = mh-1 modd = 1 if (mh*2 .eq. m) modd = 2 do 119 j=1,n do 115 i=1,mhm1 mhpi = mh+i mhmi = mh-i w(i) = y(mhmi,j)-y(mhpi,j) w(mhpi) = y(mhmi,j)+y(mhpi,j) 115 continue w(mh) = 2.*y(mh,j) go to (117,116),modd 116 w(m) = 2.*y(m,j) 117 continue do 118 i=1,m y(i,j) = w(i) 118 continue 119 continue k = iwbc+mhm1-1 i = iwba+mhm1 w(k) = (0.,0.) w(i) = (0.,0.) w(k+1) = 2.*w(k+1) go to (120,121),modd 120 continue k = iwbb+mhm1-1 w(k) = w(k)-w(i-1) w(iwbc-1) = w(iwbc-1)+w(iwbb-1) go to 122 121 w(iwbb-1) = w(k+1) 122 continue go to 107 c c reverse columns when nperod = 4 c 123 irev = 1 nby2 = n/2 124 do 126 j=1,nby2 mskip = n+1-j do 125 i=1,m a1 = y(i,j) y(i,j) = y(i,mskip) y(i,mskip) = a1 125 continue 126 continue go to (110,113),irev 127 continue do 132 j=1,n do 128 i=1,mhm1 mhmi = mh-i mhpi = mh+i w(mhmi) = .5*(y(mhpi,j)+y(i,j)) w(mhpi) = .5*(y(mhpi,j)-y(i,j)) 128 continue w(mh) = .5*y(mh,j) go to (130,129),modd 129 w(m) = .5*y(m,j) 130 continue do 131 i=1,m y(i,j) = w(i) 131 continue 132 continue 133 continue c c return storage requirements for w array. c w(1) = cmplx(float(ipstor+iwp-1),0.) return end subroutine cmpcsg (n,ijump,fnum,fden,a) complex a dimension a(1) c c c this subroutine computes required cosine values in ascending c order. when ijump .gt. 1 the routine computes values c c 2*cos(j*pi/l) , j=1,2,...,l and j .ne. 0(mod n/ijump+1) c c where l = ijump*(n/ijump+1). c c c when ijump = 1 it computes c c 2*cos((j-fnum)*pi/(n+fden)) , j=1, 2, ... ,n c c where c fnum = 0.5, fden = 0.0, for regular reduction values c fnum = 0.0, fden = 1.0, for b-r and c-r when istag = 1 c fnum = 0.0, fden = 0.5, for b-r and c-r when istag = 2 c fnum = 0.5, fden = 0.5, for b-r and c-r when istag = 2 c in cmposn only. c c pi = pimach() if (n .eq. 0) go to 105 if (ijump .eq. 1) go to 103 k3 = n/ijump+1 k4 = k3-1 pibyn = pi/float(n+ijump) do 102 k=1,ijump k1 = (k-1)*k3 k5 = (k-1)*k4 do 101 i=1,k4 x = k1+i k2 = k5+i a(k2) = cmplx(-2.*cos(x*pibyn),0.) 101 continue 102 continue go to 105 103 continue np1 = n+1 y = pi/(float(n)+fden) do 104 i=1,n x = float(np1-i)-fnum a(i) = cmplx(2.*cos(x*y),0.) 104 continue 105 continue return end subroutine cmpmrg (tcos,i1,m1,i2,m2,i3) complex tcos ,x ,y dimension tcos(1) c c c this subroutine merges two ascending strings of numbers in the c array tcos. the first string is of length m1 and starts at c tcos(i1+1). the second string is of length m2 and starts at c tcos(i2+1). the merged string goes into tcos(i3+1). c c j1 = 1 j2 = 1 j = i3 if (m1 .eq. 0) go to 107 if (m2 .eq. 0) go to 104 101 j = j+1 l = j1+i1 x = tcos(l) l = j2+i2 y = tcos(l) if (real(x-y)) 102,102,103 102 tcos(j) = x j1 = j1+1 if (j1 .gt. m1) go to 106 go to 101 103 tcos(j) = y j2 = j2+1 if (j2 .le. m2) go to 101 if (j1 .gt. m1) go to 109 104 k = j-j1+1 do 105 j=j1,m1 m = k+j l = j+i1 tcos(m) = tcos(l) 105 continue go to 109 106 continue if (j2 .gt. m2) go to 109 107 k = j-j2+1 do 108 j=j2,m2 m = k+j l = j+i2 tcos(m) = tcos(l) 108 continue 109 continue return end subroutine cmposd (mr,nr,istag,ba,bb,bc,q,idimq,b,w,d,tcos,p) c c subroutine to solve poisson's equation for dirichlet boundary c conditions. c c istag = 1 if the last diagonal block is the matrix a. c istag = 2 if the last diagonal block is the matrix a+i. c complex ba ,bb ,bc ,q , & b ,w ,d ,tcos , & p ,t dimension q(idimq,1) ,ba(1) ,bb(1) ,bc(1) , & tcos(1) ,b(1) ,d(1) ,w(1) , & p(1) m = mr n = nr fi = 1./float(istag) ip = -m ipstor = 0 jsh = 0 go to (101,102),istag 101 kr = 0 irreg = 1 if (n .gt. 1) go to 106 tcos(1) = (0.,0.) go to 103 102 kr = 1 jstsav = 1 irreg = 2 if (n .gt. 1) go to 106 tcos(1) = cmplx(-1.,0.) 103 do 104 i=1,m b(i) = q(i,1) 104 continue call cmptrx (1,0,m,ba,bb,bc,b,tcos,d,w) do 105 i=1,m q(i,1) = b(i) 105 continue go to 183 106 lr = 0 do 107 i=1,m p(i) = cmplx(0.,0.) 107 continue nun = n jst = 1 jsp = n c c irreg = 1 when no irregularities have occurred, otherwise it is 2. c 108 l = 2*jst nodd = 2-2*((nun+1)/2)+nun c c nodd = 1 when nun is odd, otherwise it is 2. c go to (110,109),nodd 109 jsp = jsp-l go to 111 110 jsp = jsp-jst if (irreg .ne. 1) jsp = jsp-l 111 continue c c regular reduction c call cmpcsg (jst,1,0.5,0.0,tcos) if (l .gt. jsp) go to 118 do 117 j=l,jsp,l jm1 = j-jsh jp1 = j+jsh jm2 = j-jst jp2 = j+jst jm3 = jm2-jsh jp3 = jp2+jsh if (jst .ne. 1) go to 113 do 112 i=1,m b(i) = 2.*q(i,j) q(i,j) = q(i,jm2)+q(i,jp2) 112 continue go to 115 113 do 114 i=1,m t = q(i,j)-q(i,jm1)-q(i,jp1)+q(i,jm2)+q(i,jp2) b(i) = t+q(i,j)-q(i,jm3)-q(i,jp3) q(i,j) = t 114 continue 115 continue call cmptrx (jst,0,m,ba,bb,bc,b,tcos,d,w) do 116 i=1,m q(i,j) = q(i,j)+b(i) 116 continue 117 continue c c reduction for last unknown c 118 go to (119,136),nodd 119 go to (152,120),irreg c c odd number of unknowns c 120 jsp = jsp+l j = jsp jm1 = j-jsh jp1 = j+jsh jm2 = j-jst jp2 = j+jst jm3 = jm2-jsh go to (123,121),istag 121 continue if (jst .ne. 1) go to 123 do 122 i=1,m b(i) = q(i,j) q(i,j) = cmplx(0.,0.) 122 continue go to 130 123 go to (124,126),noddpr 124 do 125 i=1,m ip1 = ip+i b(i) = .5*(q(i,jm2)-q(i,jm1)-q(i,jm3))+p(ip1)+q(i,j) 125 continue go to 128 126 do 127 i=1,m b(i) = .5*(q(i,jm2)-q(i,jm1)-q(i,jm3))+q(i,jp2)-q(i,jp1)+q(i,j) 127 continue 128 do 129 i=1,m q(i,j) = .5*(q(i,j)-q(i,jm1)-q(i,jp1)) 129 continue 130 call cmptrx (jst,0,m,ba,bb,bc,b,tcos,d,w) ip = ip+m ipstor = max0(ipstor,ip+m) do 131 i=1,m ip1 = ip+i p(ip1) = q(i,j)+b(i) b(i) = q(i,jp2)+p(ip1) 131 continue if (lr .ne. 0) go to 133 do 132 i=1,jst krpi = kr+i tcos(krpi) = tcos(i) 132 continue go to 134 133 continue call cmpcsg (lr,jstsav,0.,fi,tcos(jst+1)) call cmpmrg (tcos,0,jst,jst,lr,kr) 134 continue call cmpcsg (kr,jstsav,0.0,fi,tcos) call cmptrx (kr,kr,m,ba,bb,bc,b,tcos,d,w) do 135 i=1,m ip1 = ip+i q(i,j) = q(i,jm2)+b(i)+p(ip1) 135 continue lr = kr kr = kr+l go to 152 c c even number of unknowns c 136 jsp = jsp+l j = jsp jm1 = j-jsh jp1 = j+jsh jm2 = j-jst jp2 = j+jst jm3 = jm2-jsh go to (137,138),irreg 137 continue jstsav = jst ideg = jst kr = l go to 139 138 call cmpcsg (kr,jstsav,0.0,fi,tcos) call cmpcsg (lr,jstsav,0.0,fi,tcos(kr+1)) ideg = kr kr = kr+jst 139 if (jst .ne. 1) go to 141 irreg = 2 do 140 i=1,m b(i) = q(i,j) q(i,j) = q(i,jm2) 140 continue go to 150 141 do 142 i=1,m b(i) = q(i,j)+.5*(q(i,jm2)-q(i,jm1)-q(i,jm3)) 142 continue go to (143,145),irreg 143 do 144 i=1,m q(i,j) = q(i,jm2)+.5*(q(i,j)-q(i,jm1)-q(i,jp1)) 144 continue irreg = 2 go to 150 145 continue go to (146,148),noddpr 146 do 147 i=1,m ip1 = ip+i q(i,j) = q(i,jm2)+p(ip1) 147 continue ip = ip-m go to 150 148 do 149 i=1,m q(i,j) = q(i,jm2)+q(i,j)-q(i,jm1) 149 continue 150 call cmptrx (ideg,lr,m,ba,bb,bc,b,tcos,d,w) do 151 i=1,m q(i,j) = q(i,j)+b(i) 151 continue 152 nun = nun/2 noddpr = nodd jsh = jst jst = 2*jst if (nun .ge. 2) go to 108 c c start solution. c j = jsp do 153 i=1,m b(i) = q(i,j) 153 continue go to (154,155),irreg 154 continue call cmpcsg (jst,1,0.5,0.0,tcos) ideg = jst go to 156 155 kr = lr+jst call cmpcsg (kr,jstsav,0.0,fi,tcos) call cmpcsg (lr,jstsav,0.0,fi,tcos(kr+1)) ideg = kr 156 continue call cmptrx (ideg,lr,m,ba,bb,bc,b,tcos,d,w) jm1 = j-jsh jp1 = j+jsh go to (157,159),irreg 157 do 158 i=1,m q(i,j) = .5*(q(i,j)-q(i,jm1)-q(i,jp1))+b(i) 158 continue go to 164 159 go to (160,162),noddpr 160 do 161 i=1,m ip1 = ip+i q(i,j) = p(ip1)+b(i) 161 continue ip = ip-m go to 164 162 do 163 i=1,m q(i,j) = q(i,j)-q(i,jm1)+b(i) 163 continue 164 continue c c start back substitution. c jst = jst/2 jsh = jst/2 nun = 2*nun if (nun .gt. n) go to 183 do 182 j=jst,n,l jm1 = j-jsh jp1 = j+jsh jm2 = j-jst jp2 = j+jst if (j .gt. jst) go to 166 do 165 i=1,m b(i) = q(i,j)+q(i,jp2) 165 continue go to 170 166 if (jp2 .le. n) go to 168 do 167 i=1,m b(i) = q(i,j)+q(i,jm2) 167 continue if (jst .lt. jstsav) irreg = 1 go to (170,171),irreg 168 do 169 i=1,m b(i) = q(i,j)+q(i,jm2)+q(i,jp2) 169 continue 170 continue call cmpcsg (jst,1,0.5,0.0,tcos) ideg = jst jdeg = 0 go to 172 171 if (j+l .gt. n) lr = lr-jst kr = jst+lr call cmpcsg (kr,jstsav,0.0,fi,tcos) call cmpcsg (lr,jstsav,0.0,fi,tcos(kr+1)) ideg = kr jdeg = lr 172 continue call cmptrx (ideg,jdeg,m,ba,bb,bc,b,tcos,d,w) if (jst .gt. 1) go to 174 do 173 i=1,m q(i,j) = b(i) 173 continue go to 182 174 if (jp2 .gt. n) go to 177 175 do 176 i=1,m q(i,j) = .5*(q(i,j)-q(i,jm1)-q(i,jp1))+b(i) 176 continue go to 182 177 go to (175,178),irreg 178 if (j+jsh .gt. n) go to 180 do 179 i=1,m ip1 = ip+i q(i,j) = b(i)+p(ip1) 179 continue ip = ip-m go to 182 180 do 181 i=1,m q(i,j) = b(i)+q(i,j)-q(i,jm1) 181 continue 182 continue l = l/2 go to 164 183 continue c c return storage requirements for p vectors. c w(1) = cmplx(float(ipstor),0.) return end subroutine cmposn (m,n,istag,mixbnd,a,bb,c,q,idimq,b,b2,b3,w,w2, & w3,d,tcos,p) c c subroutine to solve poisson's equation with neumann boundary c conditions. c c istag = 1 if the last diagonal block is a. c istag = 2 if the last diagonal block is a-i. c mixbnd = 1 if have neumann boundary conditions at both boundaries. c mixbnd = 2 if have neumann boundary conditions at bottom and c dirichlet condition at top. (for this case, must have istag = 1.) c complex a ,bb ,c ,q , & b ,b2 ,b3 ,w , & w2 ,w3 ,d ,tcos , & p ,fi ,t dimension a(1) ,bb(1) ,c(1) ,q(idimq,1) , & b(1) ,b2(1) ,b3(1) ,w(1) , & w2(1) ,w3(1) ,d(1) ,tcos(1) , & k(4) ,p(1) equivalence (k(1),k1) ,(k(2),k2) ,(k(3),k3) ,(k(4),k4) fistag = 3-istag fnum = 1./float(istag) fden = 0.5*float(istag-1) mr = m ip = -mr ipstor = 0 i2r = 1 jr = 2 nr = n nlast = n kr = 1 lr = 0 go to (101,103),istag 101 continue do 102 i=1,mr q(i,n) = .5*q(i,n) 102 continue go to (103,104),mixbnd 103 if (n .le. 3) go to 155 104 continue jr = 2*i2r nrod = 1 if ((nr/2)*2 .eq. nr) nrod = 0 go to (105,106),mixbnd 105 jstart = 1 go to 107 106 jstart = jr nrod = 1-nrod 107 continue jstop = nlast-jr if (nrod .eq. 0) jstop = jstop-i2r call cmpcsg (i2r,1,0.5,0.0,tcos) i2rby2 = i2r/2 if (jstop .ge. jstart) go to 108 j = jr go to 116 108 continue c c regular reduction. c do 115 j=jstart,jstop,jr jp1 = j+i2rby2 jp2 = j+i2r jp3 = jp2+i2rby2 jm1 = j-i2rby2 jm2 = j-i2r jm3 = jm2-i2rby2 if (j .ne. 1) go to 109 jm1 = jp1 jm2 = jp2 jm3 = jp3 109 continue if (i2r .ne. 1) go to 111 if (j .eq. 1) jm2 = jp2 do 110 i=1,mr b(i) = 2.*q(i,j) q(i,j) = q(i,jm2)+q(i,jp2) 110 continue go to 113 111 continue do 112 i=1,mr fi = q(i,j) q(i,j) = q(i,j)-q(i,jm1)-q(i,jp1)+q(i,jm2)+q(i,jp2) b(i) = fi+q(i,j)-q(i,jm3)-q(i,jp3) 112 continue 113 continue call cmptrx (i2r,0,mr,a,bb,c,b,tcos,d,w) do 114 i=1,mr q(i,j) = q(i,j)+b(i) 114 continue c c end of reduction for regular unknowns. c 115 continue c c begin special reduction for last unknown. c j = jstop+jr 116 nlast = j jm1 = j-i2rby2 jm2 = j-i2r jm3 = jm2-i2rby2 if (nrod .eq. 0) go to 128 c c odd number of unknowns c if (i2r .ne. 1) go to 118 do 117 i=1,mr b(i) = fistag*q(i,j) q(i,j) = q(i,jm2) 117 continue go to 126 118 do 119 i=1,mr b(i) = q(i,j)+.5*(q(i,jm2)-q(i,jm1)-q(i,jm3)) 119 continue if (nrodpr .ne. 0) go to 121 do 120 i=1,mr ii = ip+i q(i,j) = q(i,jm2)+p(ii) 120 continue ip = ip-mr go to 123 121 continue do 122 i=1,mr q(i,j) = q(i,j)-q(i,jm1)+q(i,jm2) 122 continue 123 if (lr .eq. 0) go to 124 call cmpcsg (lr,1,0.5,fden,tcos(kr+1)) go to 126 124 continue do 125 i=1,mr b(i) = fistag*b(i) 125 continue 126 continue call cmpcsg (kr,1,0.5,fden,tcos) call cmptrx (kr,lr,mr,a,bb,c,b,tcos,d,w) do 127 i=1,mr q(i,j) = q(i,j)+b(i) 127 continue kr = kr+i2r go to 151 128 continue c c even number of unknowns c jp1 = j+i2rby2 jp2 = j+i2r if (i2r .ne. 1) go to 135 do 129 i=1,mr b(i) = q(i,j) 129 continue call cmptrx (1,0,mr,a,bb,c,b,tcos,d,w) ip = 0 ipstor = mr go to (133,130),istag 130 do 131 i=1,mr p(i) = b(i) b(i) = b(i)+q(i,n) 131 continue tcos(1) = cmplx(1.,0.) tcos(2) = cmplx(0.,0.) call cmptrx (1,1,mr,a,bb,c,b,tcos,d,w) do 132 i=1,mr q(i,j) = q(i,jm2)+p(i)+b(i) 132 continue go to 150 133 continue do 134 i=1,mr p(i) = b(i) q(i,j) = q(i,jm2)+2.*q(i,jp2)+3.*b(i) 134 continue go to 150 135 continue do 136 i=1,mr b(i) = q(i,j)+.5*(q(i,jm2)-q(i,jm1)-q(i,jm3)) 136 continue if (nrodpr .ne. 0) go to 138 do 137 i=1,mr ii = ip+i b(i) = b(i)+p(ii) 137 continue go to 140 138 continue do 139 i=1,mr b(i) = b(i)+q(i,jp2)-q(i,jp1) 139 continue 140 continue call cmptrx (i2r,0,mr,a,bb,c,b,tcos,d,w) ip = ip+mr ipstor = max0(ipstor,ip+mr) do 141 i=1,mr ii = ip+i p(ii) = b(i)+.5*(q(i,j)-q(i,jm1)-q(i,jp1)) b(i) = p(ii)+q(i,jp2) 141 continue if (lr .eq. 0) go to 142 call cmpcsg (lr,1,0.5,fden,tcos(i2r+1)) call cmpmrg (tcos,0,i2r,i2r,lr,kr) go to 144 142 do 143 i=1,i2r ii = kr+i tcos(ii) = tcos(i) 143 continue 144 call cmpcsg (kr,1,0.5,fden,tcos) if (lr .ne. 0) go to 145 go to (146,145),istag 145 continue call cmptrx (kr,kr,mr,a,bb,c,b,tcos,d,w) go to 148 146 continue do 147 i=1,mr b(i) = fistag*b(i) 147 continue 148 continue do 149 i=1,mr ii = ip+i q(i,j) = q(i,jm2)+p(ii)+b(i) 149 continue 150 continue lr = kr kr = kr+jr 151 continue go to (152,153),mixbnd 152 nr = (nlast-1)/jr+1 if (nr .le. 3) go to 155 go to 154 153 nr = nlast/jr if (nr .le. 1) go to 192 154 i2r = jr nrodpr = nrod go to 104 155 continue c c begin solution c j = 1+jr jm1 = j-i2r jp1 = j+i2r jm2 = nlast-i2r if (nr .eq. 2) go to 184 if (lr .ne. 0) go to 170 if (n .ne. 3) go to 161 c c case n = 3. c go to (156,168),istag 156 continue do 157 i=1,mr b(i) = q(i,2) 157 continue tcos(1) = cmplx(0.,0.) call cmptrx (1,0,mr,a,bb,c,b,tcos,d,w) do 158 i=1,mr q(i,2) = b(i) b(i) = 4.*b(i)+q(i,1)+2.*q(i,3) 158 continue tcos(1) = cmplx(-2.,0.) tcos(2) = cmplx(2.,0.) i1 = 2 i2 = 0 call cmptrx (i1,i2,mr,a,bb,c,b,tcos,d,w) do 159 i=1,mr q(i,2) = q(i,2)+b(i) b(i) = q(i,1)+2.*q(i,2) 159 continue tcos(1) = (0.,0.) call cmptrx (1,0,mr,a,bb,c,b,tcos,d,w) do 160 i=1,mr q(i,1) = b(i) 160 continue jr = 1 i2r = 0 go to 194 c c case n = 2**p+1 c 161 continue go to (162,170),istag 162 continue do 163 i=1,mr b(i) = q(i,j)+.5*q(i,1)-q(i,jm1)+q(i,nlast)-q(i,jm2) 163 continue call cmpcsg (jr,1,0.5,0.0,tcos) call cmptrx (jr,0,mr,a,bb,c,b,tcos,d,w) do 164 i=1,mr q(i,j) = .5*(q(i,j)-q(i,jm1)-q(i,jp1))+b(i) b(i) = q(i,1)+2.*q(i,nlast)+4.*q(i,j) 164 continue jr2 = 2*jr call cmpcsg (jr,1,0.0,0.0,tcos) do 165 i=1,jr i1 = jr+i i2 = jr+1-i tcos(i1) = -tcos(i2) 165 continue call cmptrx (jr2,0,mr,a,bb,c,b,tcos,d,w) do 166 i=1,mr q(i,j) = q(i,j)+b(i) b(i) = q(i,1)+2.*q(i,j) 166 continue call cmpcsg (jr,1,0.5,0.0,tcos) call cmptrx (jr,0,mr,a,bb,c,b,tcos,d,w) do 167 i=1,mr q(i,1) = .5*q(i,1)-q(i,jm1)+b(i) 167 continue go to 194 c c case of general n with nr = 3 . c 168 do 169 i=1,mr b(i) = q(i,2) q(i,2) = (0.,0.) b2(i) = q(i,3) b3(i) = q(i,1) 169 continue jr = 1 i2r = 0 j = 2 go to 177 170 continue do 171 i=1,mr b(i) = .5*q(i,1)-q(i,jm1)+q(i,j) 171 continue if (nrod .ne. 0) go to 173 do 172 i=1,mr ii = ip+i b(i) = b(i)+p(ii) 172 continue go to 175 173 do 174 i=1,mr b(i) = b(i)+q(i,nlast)-q(i,jm2) 174 continue 175 continue do 176 i=1,mr t = .5*(q(i,j)-q(i,jm1)-q(i,jp1)) q(i,j) = t b2(i) = q(i,nlast)+t b3(i) = q(i,1)+2.*t 176 continue 177 continue k1 = kr+2*jr-1 k2 = kr+jr tcos(k1+1) = (-2.,0.) k4 = k1+3-istag call cmpcsg (k2+istag-2,1,0.0,fnum,tcos(k4)) k4 = k1+k2+1 call cmpcsg (jr-1,1,0.0,1.0,tcos(k4)) call cmpmrg (tcos,k1,k2,k1+k2,jr-1,0) k3 = k1+k2+lr call cmpcsg (jr,1,0.5,0.0,tcos(k3+1)) k4 = k3+jr+1 call cmpcsg (kr,1,0.5,fden,tcos(k4)) call cmpmrg (tcos,k3,jr,k3+jr,kr,k1) if (lr .eq. 0) go to 178 call cmpcsg (lr,1,0.5,fden,tcos(k4)) call cmpmrg (tcos,k3,jr,k3+jr,lr,k3-lr) call cmpcsg (kr,1,0.5,fden,tcos(k4)) 178 k3 = kr k4 = kr call cmptr3 (mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3) do 179 i=1,mr b(i) = b(i)+b2(i)+b3(i) 179 continue tcos(1) = (2.,0.) call cmptrx (1,0,mr,a,bb,c,b,tcos,d,w) do 180 i=1,mr q(i,j) = q(i,j)+b(i) b(i) = q(i,1)+2.*q(i,j) 180 continue call cmpcsg (jr,1,0.5,0.0,tcos) call cmptrx (jr,0,mr,a,bb,c,b,tcos,d,w) if (jr .ne. 1) go to 182 do 181 i=1,mr q(i,1) = b(i) 181 continue go to 194 182 continue do 183 i=1,mr q(i,1) = .5*q(i,1)-q(i,jm1)+b(i) 183 continue go to 194 184 continue if (n .ne. 2) go to 188 c c case n = 2 c do 185 i=1,mr b(i) = q(i,1) 185 continue tcos(1) = (0.,0.) call cmptrx (1,0,mr,a,bb,c,b,tcos,d,w) do 186 i=1,mr q(i,1) = b(i) b(i) = 2.*(q(i,2)+b(i))*fistag 186 continue tcos(1) = cmplx(-fistag,0.) tcos(2) = cmplx(2.,0.) call cmptrx (2,0,mr,a,bb,c,b,tcos,d,w) do 187 i=1,mr q(i,1) = q(i,1)+b(i) 187 continue jr = 1 i2r = 0 go to 194 188 continue c c case of general n and nr = 2 . c do 189 i=1,mr ii = ip+i b3(i) = (0.,0.) b(i) = q(i,1)+2.*p(ii) q(i,1) = .5*q(i,1)-q(i,jm1) b2(i) = 2.*(q(i,1)+q(i,nlast)) 189 continue k1 = kr+jr-1 tcos(k1+1) = (-2.,0.) k4 = k1+3-istag call cmpcsg (kr+istag-2,1,0.0,fnum,tcos(k4)) k4 = k1+kr+1 call cmpcsg (jr-1,1,0.0,1.0,tcos(k4)) call cmpmrg (tcos,k1,kr,k1+kr,jr-1,0) call cmpcsg (kr,1,0.5,fden,tcos(k1+1)) k2 = kr k4 = k1+k2+1 call cmpcsg (lr,1,0.5,fden,tcos(k4)) k3 = lr k4 = 0 call cmptr3 (mr,a,bb,c,k,b,b2,b3,tcos,d,w,w2,w3) do 190 i=1,mr b(i) = b(i)+b2(i) 190 continue tcos(1) = (2.,0.) call cmptrx (1,0,mr,a,bb,c,b,tcos,d,w) do 191 i=1,mr q(i,1) = q(i,1)+b(i) 191 continue go to 194 192 do 193 i=1,mr b(i) = q(i,nlast) 193 continue go to 196 194 continue c c start back substitution. c j = nlast-jr do 195 i=1,mr b(i) = q(i,nlast)+q(i,j) 195 continue 196 jm2 = nlast-i2r if (jr .ne. 1) go to 198 do 197 i=1,mr q(i,nlast) = (0.,0.) 197 continue go to 202 198 continue if (nrod .ne. 0) go to 200 do 199 i=1,mr ii = ip+i q(i,nlast) = p(ii) 199 continue ip = ip-mr go to 202 200 do 201 i=1,mr q(i,nlast) = q(i,nlast)-q(i,jm2) 201 continue 202 continue call cmpcsg (kr,1,0.5,fden,tcos) call cmpcsg (lr,1,0.5,fden,tcos(kr+1)) if (lr .ne. 0) go to 204 do 203 i=1,mr b(i) = fistag*b(i) 203 continue 204 continue call cmptrx (kr,lr,mr,a,bb,c,b,tcos,d,w) do 205 i=1,mr q(i,nlast) = q(i,nlast)+b(i) 205 continue nlastp = nlast 206 continue jstep = jr jr = i2r i2r = i2r/2 if (jr .eq. 0) go to 222 go to (207,208),mixbnd 207 jstart = 1+jr go to 209 208 jstart = jr 209 continue kr = kr-jr if (nlast+jr .gt. n) go to 210 kr = kr-jr nlast = nlast+jr jstop = nlast-jstep go to 211 210 continue jstop = nlast-jr 211 continue lr = kr-jr call cmpcsg (jr,1,0.5,0.0,tcos) do 221 j=jstart,jstop,jstep jm2 = j-jr jp2 = j+jr if (j .ne. jr) go to 213 do 212 i=1,mr b(i) = q(i,j)+q(i,jp2) 212 continue go to 215 213 continue do 214 i=1,mr b(i) = q(i,j)+q(i,jm2)+q(i,jp2) 214 continue 215 continue if (jr .ne. 1) go to 217 do 216 i=1,mr q(i,j) = (0.,0.) 216 continue go to 219 217 continue jm1 = j-i2r jp1 = j+i2r do 218 i=1,mr q(i,j) = .5*(q(i,j)-q(i,jm1)-q(i,jp1)) 218 continue 219 continue call cmptrx (jr,0,mr,a,bb,c,b,tcos,d,w) do 220 i=1,mr q(i,j) = q(i,j)+b(i) 220 continue 221 continue nrod = 1 if (nlast+i2r .le. n) nrod = 0 if (nlastp .ne. nlast) go to 194 go to 206 222 continue c c return storage requirements for p vectors. c w(1) = cmplx(float(ipstor),0.) return end subroutine cmposp (m,n,a,bb,c,q,idimq,b,b2,b3,w,w2,w3,d,tcos,p) c c subroutine to solve poisson equation with periodic boundary c conditions. c complex a ,bb ,c ,q , & b ,b2 ,b3 ,w , & w2 ,w3 ,d ,tcos , & p ,s ,t dimension a(1) ,bb(1) ,c(1) ,q(idimq,1) , & b(1) ,b2(1) ,b3(1) ,w(1) , & w2(1) ,w3(1) ,d(1) ,tcos(1) , & p(1) mr = m nr = (n+1)/2 nrm1 = nr-1 if (2*nr .ne. n) go to 107 c c even number of unknowns c do 102 j=1,nrm1 nrmj = nr-j nrpj = nr+j do 101 i=1,mr s = q(i,nrmj)-q(i,nrpj) t = q(i,nrmj)+q(i,nrpj) q(i,nrmj) = s q(i,nrpj) = t 101 continue 102 continue do 103 i=1,mr q(i,nr) = 2.*q(i,nr) q(i,n) = 2.*q(i,n) 103 continue call cmposd (mr,nrm1,1,a,bb,c,q,idimq,b,w,d,tcos,p) ipstor = real(w(1)) call cmposn (mr,nr+1,1,1,a,bb,c,q(1,nr),idimq,b,b2,b3,w,w2,w3,d, & tcos,p) ipstor = max0(ipstor,int(real(w(1)))) do 105 j=1,nrm1 nrmj = nr-j nrpj = nr+j do 104 i=1,mr s = .5*(q(i,nrpj)+q(i,nrmj)) t = .5*(q(i,nrpj)-q(i,nrmj)) q(i,nrmj) = s q(i,nrpj) = t 104 continue 105 continue do 106 i=1,mr q(i,nr) = .5*q(i,nr) q(i,n) = .5*q(i,n) 106 continue go to 118 107 continue c c odd number of unknowns c do 109 j=1,nrm1 nrpj = n+1-j do 108 i=1,mr s = q(i,j)-q(i,nrpj) t = q(i,j)+q(i,nrpj) q(i,j) = s q(i,nrpj) = t 108 continue 109 continue do 110 i=1,mr q(i,nr) = 2.*q(i,nr) 110 continue lh = nrm1/2 do 112 j=1,lh nrmj = nr-j do 111 i=1,mr s = q(i,j) q(i,j) = q(i,nrmj) q(i,nrmj) = s 111 continue 112 continue call cmposd (mr,nrm1,2,a,bb,c,q,idimq,b,w,d,tcos,p) ipstor = real(w(1)) call cmposn (mr,nr,2,1,a,bb,c,q(1,nr),idimq,b,b2,b3,w,w2,w3,d, & tcos,p) ipstor = max0(ipstor,int(real(w(1)))) do 114 j=1,nrm1 nrpj = nr+j do 113 i=1,mr s = .5*(q(i,nrpj)+q(i,j)) t = .5*(q(i,nrpj)-q(i,j)) q(i,nrpj) = t q(i,j) = s 113 continue 114 continue do 115 i=1,mr q(i,nr) = .5*q(i,nr) 115 continue do 117 j=1,lh nrmj = nr-j do 116 i=1,mr s = q(i,j) q(i,j) = q(i,nrmj) q(i,nrmj) = s 116 continue 117 continue 118 continue c c return storage requirements for p vectors. c w(1) = cmplx(float(ipstor),0.) return end subroutine cmptr3 (m,a,b,c,k,y1,y2,y3,tcos,d,w1,w2,w3) complex a ,b ,c ,y1 , & y2 ,y3 ,tcos ,d , & w1 ,w2 ,w3 ,x , & xx ,z dimension a(1) ,b(1) ,c(1) ,k(4) , & tcos(1) ,y1(1) ,y2(1) ,y3(1) , & d(1) ,w1(1) ,w2(1) ,w3(1) c c subroutine to solve tridiagonal systems c mm1 = m-1 k1 = k(1) k2 = k(2) k3 = k(3) k4 = k(4) f1 = k1+1 f2 = k2+1 f3 = k3+1 f4 = k4+1 k2k3k4 = k2+k3+k4 if (k2k3k4 .eq. 0) go to 101 l1 = f1/f2 l2 = f1/f3 l3 = f1/f4 lint1 = 1 lint2 = 1 lint3 = 1 kint1 = k1 kint2 = kint1+k2 kint3 = kint2+k3 101 continue do 115 n=1,k1 x = tcos(n) if (k2k3k4 .eq. 0) go to 107 if (n .ne. l1) go to 103 do 102 i=1,m w1(i) = y1(i) 102 continue 103 if (n .ne. l2) go to 105 do 104 i=1,m w2(i) = y2(i) 104 continue 105 if (n .ne. l3) go to 107 do 106 i=1,m w3(i) = y3(i) 106 continue 107 continue z = 1./(b(1)-x) d(1) = c(1)*z y1(1) = y1(1)*z y2(1) = y2(1)*z y3(1) = y3(1)*z do 108 i=2,m z = 1./(b(i)-x-a(i)*d(i-1)) d(i) = c(i)*z y1(i) = (y1(i)-a(i)*y1(i-1))*z y2(i) = (y2(i)-a(i)*y2(i-1))*z y3(i) = (y3(i)-a(i)*y3(i-1))*z 108 continue do 109 ip=1,mm1 i = m-ip y1(i) = y1(i)-d(i)*y1(i+1) y2(i) = y2(i)-d(i)*y2(i+1) y3(i) = y3(i)-d(i)*y3(i+1) 109 continue if (k2k3k4 .eq. 0) go to 115 if (n .ne. l1) go to 111 i = lint1+kint1 xx = x-tcos(i) do 110 i=1,m y1(i) = xx*y1(i)+w1(i) 110 continue lint1 = lint1+1 l1 = (float(lint1)*f1)/f2 111 if (n .ne. l2) go to 113 i = lint2+kint2 xx = x-tcos(i) do 112 i=1,m y2(i) = xx*y2(i)+w2(i) 112 continue lint2 = lint2+1 l2 = (float(lint2)*f1)/f3 113 if (n .ne. l3) go to 115 i = lint3+kint3 xx = x-tcos(i) do 114 i=1,m y3(i) = xx*y3(i)+w3(i) 114 continue lint3 = lint3+1 l3 = (float(lint3)*f1)/f4 115 continue return end subroutine cmptrx (idegbr,idegcr,m,a,b,c,y,tcos,d,w) c c subroutine to solve a system of linear equations where the c coefficient matrix is a rational function in the matrix given by c tridiagonal ( . . . , a(i), b(i), c(i), . . . ). c complex a ,b ,c ,y , & tcos ,d ,w ,x , & xx ,z dimension a(1) ,b(1) ,c(1) ,y(1) , & tcos(1) ,d(1) ,w(1) mm1 = m-1 fb = idegbr+1 fc = idegcr+1 l = fb/fc lint = 1 do 108 k=1,idegbr x = tcos(k) if (k .ne. l) go to 102 i = idegbr+lint xx = x-tcos(i) do 101 i=1,m w(i) = y(i) y(i) = xx*y(i) 101 continue 102 continue z = 1./(b(1)-x) d(1) = c(1)*z y(1) = y(1)*z do 103 i=2,mm1 z = 1./(b(i)-x-a(i)*d(i-1)) d(i) = c(i)*z y(i) = (y(i)-a(i)*y(i-1))*z 103 continue z = b(m)-x-a(m)*d(mm1) if (cabs(z) .ne. 0.) go to 104 y(m) = (0.,0.) go to 105 104 y(m) = (y(m)-a(m)*y(mm1))/z 105 continue do 106 ip=1,mm1 i = m-ip y(i) = y(i)-d(i)*y(i+1) 106 continue if (k .ne. l) go to 108 do 107 i=1,m y(i) = y(i)+w(i) 107 continue lint = lint+1 l = (float(lint)*fb)/fc 108 continue return end subroutine cofx (x,af,bf,cf) c c set coefficients in the x-direction. c af = (x+1.)**2 bf = 2.0*(x+1.) cf = -x return end subroutine cofx4(x,af,bf,cf) c c set coefficients in the x-direction. c af = (x+1.)**2 bf = 2.0*(x+1.) cf = -x return end subroutine cofy (y,df,ef,ff) c c set coefficients in y direction c df = exp(y) ef = 0.0 ff = -y return end subroutine compb (n,ierror,an,bn,cn,b,ah,bh) c c compb computes the roots of the b polynomials using subroutine c tevls which is a modification the eispack program tqlrat. c ierror is set to 4 if either tevls fails or if a(j+1)*c(j) is c less than zero for some j. ah,bh are temporary work arrays. c dimension an(1) ,bn(1) ,cn(1) ,b(1) , & ah(1) ,bh(1) common /cblkt/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik eps = epmach() bnorm = abs(bn(1)) do 102 j=2,nm bnorm = amax1(bnorm,abs(bn(j))) arg = an(j)*cn(j-1) if (arg) 119,101,101 101 b(j) = sign(sqrt(arg),an(j)) 102 continue cnv = eps*bnorm if = 2**k kdo = k-1 do 108 l=1,kdo ir = l-1 i2 = 2**ir i4 = i2+i2 ipl = i4-1 ifd = if-i4 do 107 i=i4,ifd,i4 call indxb (i,l,ib,nb) if (nb) 108,108,103 103 js = i-ipl jf = js+nb-1 ls = 0 do 104 j=js,jf ls = ls+1 bh(ls) = bn(j) ah(ls) = b(j) 104 continue call tevls (nb,bh,ah,ierror) if (ierror) 118,105,118 105 lh = ib-1 do 106 j=1,nb lh = lh+1 b(lh) = -bh(j) 106 continue 107 continue 108 continue do 109 j=1,nm b(j) = -bn(j) 109 continue if (npp) 117,110,117 110 nmp = nm+1 nb = nm+nmp do 112 j=1,nb l1 = mod(j-1,nmp)+1 l2 = mod(j+nm-1,nmp)+1 arg = an(l1)*cn(l2) if (arg) 119,111,111 111 bh(j) = sign(sqrt(arg),-an(l1)) ah(j) = -bn(l1) 112 continue call tevls (nb,ah,bh,ierror) if (ierror) 118,113,118 113 call indxb (if,k-1,j2,lh) call indxb (if/2,k-1,j1,lh) j2 = j2+1 lh = j2 n2m2 = j2+nm+nm-2 114 d1 = abs(b(j1)-b(j2-1)) d2 = abs(b(j1)-b(j2)) d3 = abs(b(j1)-b(j2+1)) if ((d2 .lt. d1) .and. (d2 .lt. d3)) go to 115 b(lh) = b(j2) j2 = j2+1 lh = lh+1 if (j2-n2m2) 114,114,116 115 j2 = j2+1 j1 = j1+1 if (j2-n2m2) 114,114,116 116 b(lh) = b(n2m2+1) call indxb (if,k-1,j1,j2) j2 = j1+nmp+nmp call ppadd (nm+1,ierror,an,cn,b(j1),b(j1),b(j2)) 117 return 118 ierror = 4 return 119 ierror = 5 return end subroutine cosgen (n,ijump,fnum,fden,a) dimension a(1) c c c this subroutine computes required cosine values in ascending c order. when ijump .gt. 1 the routine computes values c c 2*cos(j*pi/l) , j=1,2,...,l and j .ne. 0(mod n/ijump+1) c c where l = ijump*(n/ijump+1). c c c when ijump = 1 it computes c c 2*cos((j-fnum)*pi/(n+fden)) , j=1, 2, ... ,n c c where c fnum = 0.5, fden = 0.0, for regular reduction values c fnum = 0.0, fden = 1.0, for b-r and c-r when istag = 1 c fnum = 0.0, fden = 0.5, for b-r and c-r when istag = 2 c fnum = 0.5, fden = 0.5, for b-r and c-r when istag = 2 c in poisn2 only. c c pi = pimach() if (n .eq. 0) go to 105 if (ijump .eq. 1) go to 103 k3 = n/ijump+1 k4 = k3-1 pibyn = pi/float(n+ijump) do 102 k=1,ijump k1 = (k-1)*k3 k5 = (k-1)*k4 do 101 i=1,k4 x = k1+i k2 = k5+i a(k2) = -2.*cos(x*pibyn) 101 continue 102 continue go to 105 103 continue np1 = n+1 y = pi/(float(n)+fden) do 104 i=1,n x = float(np1-i)-fnum a(i) = 2.*cos(x*y) 104 continue 105 continue return end subroutine cpadd (n,ierror,a,c,cbp,bp,bh) c c cpadd computes the eigenvalues of the periodic tridiagonal matrix c with coefficients an,bn,cn c c n is the order of the bh and bp polynomials c on output bp contians the eigenvalues c cbp is the same as bp except type complex c bh is used to temporarily store the roots of the b hat polynomial c which enters through bp c complex cf ,cx ,fsg ,hsg , & dd ,f ,fp ,fpp , & cdis ,r1 ,r2 ,r3 , & cbp dimension a(1) ,c(1) ,bp(1) ,bh(1) , & cbp(1) common /ccblk/ npp ,k ,eps ,cnv , & nm ,ncmplx ,ik external pgsf ,pppsf ,ppgsf scnv = sqrt(cnv) iz = n izm = iz-1 izm2 = iz-2 if (bp(n)-bp(1)) 101,142,103 101 do 102 j=1,n nt = n-j bh(j) = bp(nt+1) 102 continue go to 105 103 do 104 j=1,n bh(j) = bp(j) 104 continue 105 ncmplx = 0 modiz = mod(iz,2) is = 1 if (modiz) 106,107,106 106 if (a(1)) 110,142,107 107 xl = bh(1) db = bh(3)-bh(1) 108 xl = xl-db if (pgsf(xl,iz,c,a,bh)) 108,108,109 109 sgn = -1. cbp(1) = cmplx(bcrh(xl,bh(1),iz,c,a,bh,pgsf,sgn),0.) is = 2 110 if = iz-1 if (modiz) 111,112,111 111 if (a(1)) 112,142,115 112 xr = bh(iz) db = bh(iz)-bh(iz-2) 113 xr = xr+db if (pgsf(xr,iz,c,a,bh)) 113,114,114 114 sgn = 1. cbp(iz) = cmplx(bcrh(bh(iz),xr,iz,c,a,bh,pgsf,sgn),0.) if = iz-2 115 do 136 ig=is,if,2 xl = bh(ig) xr = bh(ig+1) sgn = -1. xm = bcrh(xl,xr,iz,c,a,bh,pppsf,sgn) psg = pgsf(xm,iz,c,a,bh) if (abs(psg)-eps) 118,118,116 116 if (psg*ppgsf(xm,iz,c,a,bh)) 117,118,119 c c case of a real zero c 117 sgn = 1. cbp(ig) = cmplx(bcrh(bh(ig),xm,iz,c,a,bh,pgsf,sgn),0.) sgn = -1. cbp(ig+1) = cmplx(bcrh(xm,bh(ig+1),iz,c,a,bh,pgsf,sgn),0.) go to 136 c c case of a multiple zero c 118 cbp(ig) = cmplx(xm,0.) cbp(ig+1) = cmplx(xm,0.) go to 136 c c case of a complex zero c 119 it = 0 icv = 0 cx = cmplx(xm,0.) 120 fsg = (1.,0.) hsg = (1.,0.) fp = (0.,0.) fpp = (0.,0.) do 121 j=1,iz dd = 1./(cx-bh(j)) fsg = fsg*a(j)*dd hsg = hsg*c(j)*dd fp = fp+dd fpp = fpp-dd*dd 121 continue if (modiz) 123,122,123 122 f = (1.,0.)-fsg-hsg go to 124 123 f = (1.,0.)+fsg+hsg 124 i3 = 0 if (cabs(fp)) 126,126,125 125 i3 = 1 r3 = -f/fp 126 i2 = 0 if (cabs(fpp)) 132,132,127 127 i2 = 1 cdis = csqrt(fp**2-2.*f*fpp) r1 = cdis-fp r2 = -fp-cdis if (cabs(r1)-cabs(r2)) 129,129,128 128 r1 = r1/fpp go to 130 129 r1 = r2/fpp 130 r2 = 2.*f/fpp/r1 if (cabs(r2) .lt. cabs(r1)) r1 = r2 if (i3) 133,133,131 131 if (cabs(r3) .lt. cabs(r1)) r1 = r3 go to 133 132 r1 = r3 133 cx = cx+r1 it = it+1 if (it .gt. 50) go to 142 if (cabs(r1) .gt. scnv) go to 120 if (icv) 134,134,135 134 icv = 1 go to 120 135 cbp(ig) = cx cbp(ig+1) = conjg(cx) 136 continue if (cabs(cbp(n))-cabs(cbp(1))) 137,142,139 137 nhalf = n/2 do 138 j=1,nhalf nt = n-j cx = cbp(j) cbp(j) = cbp(nt+1) cbp(nt+1) = cx 138 continue 139 ncmplx = 1 do 140 j=2,iz if (aimag(cbp(j))) 143,140,143 140 continue ncmplx = 0 do 141 j=2,iz bp(j) = real(cbp(j)) 141 continue go to 143 142 ierror = 4 143 continue return end subroutine cproc (nd,bd,nm1,bm1,nm2,bm2,na,aa,x,y,m,a,b,c,d,w,yy) c c proc applies a sequence of matrix operations to the vector x and c stores the result in y c aa array containing scalar multipliers of the vector x c nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively c bd,bm1,bm2 are arrays containing roots of certian b polynomials c na is the length of the array aa c x,y the matrix operations are applied to x and the result is y c a,b,c are arrays which contain the tridiagonal matrix c m is the order of the matrix c d,w are work arrays c isgn determines whether or not a change in sign is made c complex y ,d ,w ,bd , & crt ,den ,y1 ,y2 , & x ,a ,b ,c dimension a(1) ,b(1) ,c(1) ,x(1) , & y(1) ,d(1) ,w(1) ,bd(1) , & bm1(1) ,bm2(1) ,aa(1) ,yy(1) do 101 j=1,m y(j) = x(j) 101 continue mm = m-1 id = nd m1 = nm1 m2 = nm2 ia = na 102 iflg = 0 if (id) 109,109,103 103 crt = bd(id) id = id-1 c c begin solution to system c d(m) = a(m)/(b(m)-crt) w(m) = y(m)/(b(m)-crt) do 104 j=2,mm k = m-j den = b(k+1)-crt-c(k+1)*d(k+2) d(k+1) = a(k+1)/den w(k+1) = (y(k+1)-c(k+1)*w(k+2))/den 104 continue den = b(1)-crt-c(1)*d(2) if (cabs(den)) 105,106,105 105 y(1) = (y(1)-c(1)*w(2))/den go to 107 106 y(1) = (1.,0.) 107 do 108 j=2,m y(j) = w(j)-d(j)*y(j-1) 108 continue 109 if (m1) 110,110,112 110 if (m2) 121,121,111 111 rt = bm2(m2) m2 = m2-1 go to 117 112 if (m2) 113,113,114 113 rt = bm1(m1) m1 = m1-1 go to 117 114 if (abs(bm1(m1))-abs(bm2(m2))) 116,116,115 115 rt = bm1(m1) m1 = m1-1 go to 117 116 rt = bm2(m2) m2 = m2-1 117 y1 = (b(1)-rt)*y(1)+c(1)*y(2) if (mm-2) 120,118,118 c c matrix multiplication c 118 do 119 j=2,mm y2 = a(j)*y(j-1)+(b(j)-rt)*y(j)+c(j)*y(j+1) y(j-1) = y1 y1 = y2 119 continue 120 y(m) = a(m)*y(m-1)+(b(m)-rt)*y(m) y(m-1) = y1 iflg = 1 go to 102 121 if (ia) 124,124,122 122 rt = aa(ia) ia = ia-1 iflg = 1 c c scalar multiplication c do 123 j=1,m y(j) = rt*y(j) 123 continue 124 if (iflg) 125,125,102 125 return end subroutine cprocp (nd,bd,nm1,bm1,nm2,bm2,na,aa,x,y,m,a,b,c,d,u,yy) c c cprocp applies a sequence of matrix operations to the vector x and c stores the result in y c c bd,bm1,bm2 are arrays containing roots of certian b polynomials c nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively c aa array containing scalar multipliers of the vector x c na is the length of the array aa c x,y the matrix operations are applied to x and the result is y c a,b,c are arrays which contain the tridiagonal matrix c m is the order of the matrix c d,u are work arrays c isgn determines whether or not a change in sign is made c complex y ,d ,u ,v , & den ,bh ,ym ,am , & y1 ,y2 ,yh ,bd , & crt ,x ,a ,b ,c dimension a(1) ,b(1) ,c(1) ,x(1) , & y(1) ,d(1) ,u(1) ,bd(1) , & bm1(1) ,bm2(1) ,aa(1) ,yy(1) do 101 j=1,m y(j) = x(j) 101 continue mm = m-1 mm2 = m-2 id = nd m1 = nm1 m2 = nm2 ia = na 102 iflg = 0 if (id) 111,111,103 103 crt = bd(id) id = id-1 iflg = 1 c c begin solution to system c bh = b(m)-crt ym = y(m) den = b(1)-crt d(1) = c(1)/den u(1) = a(1)/den y(1) = y(1)/den v = c(m) if (mm2-2) 106,104,104 104 do 105 j=2,mm2 den = b(j)-crt-a(j)*d(j-1) d(j) = c(j)/den u(j) = -a(j)*u(j-1)/den y(j) = (y(j)-a(j)*y(j-1))/den bh = bh-v*u(j-1) ym = ym-v*y(j-1) v = -v*d(j-1) 105 continue 106 den = b(m-1)-crt-a(m-1)*d(m-2) d(m-1) = (c(m-1)-a(m-1)*u(m-2))/den y(m-1) = (y(m-1)-a(m-1)*y(m-2))/den am = a(m)-v*d(m-2) bh = bh-v*u(m-2) ym = ym-v*y(m-2) den = bh-am*d(m-1) if (cabs(den)) 107,108,107 107 y(m) = (ym-am*y(m-1))/den go to 109 108 y(m) = (1.,0.) 109 y(m-1) = y(m-1)-d(m-1)*y(m) do 110 j=2,mm k = m-j y(k) = y(k)-d(k)*y(k+1)-u(k)*y(m) 110 continue 111 if (m1) 112,112,114 112 if (m2) 123,123,113 113 rt = bm2(m2) m2 = m2-1 go to 119 114 if (m2) 115,115,116 115 rt = bm1(m1) m1 = m1-1 go to 119 116 if (abs(bm1(m1))-abs(bm2(m2))) 118,118,117 117 rt = bm1(m1) m1 = m1-1 go to 119 118 rt = bm2(m2) m2 = m2-1 c c matrix multiplication c 119 yh = y(1) y1 = (b(1)-rt)*y(1)+c(1)*y(2)+a(1)*y(m) if (mm-2) 122,120,120 120 do 121 j=2,mm y2 = a(j)*y(j-1)+(b(j)-rt)*y(j)+c(j)*y(j+1) y(j-1) = y1 y1 = y2 121 continue 122 y(m) = a(m)*y(m-1)+(b(m)-rt)*y(m)+c(m)*yh y(m-1) = y1 iflg = 1 go to 102 123 if (ia) 126,126,124 124 rt = aa(ia) ia = ia-1 iflg = 1 c c scalar multiplication c do 125 j=1,m y(j) = rt*y(j) 125 continue 126 if (iflg) 127,127,102 127 return end subroutine cprod (nd,bd,nm1,bm1,nm2,bm2,na,aa,x,yy,m,a,b,c,d,w,y) c c prod applies a sequence of matrix operations to the vector x and c stores the result in yy (complex case) c aa array containing scalar multipliers of the vector x c nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively c bd,bm1,bm2 are arrays containing roots of certian b polynomials c na is the length of the array aa c x,yy the matrix operations are applied to x and the result is yy c a,b,c are arrays which contain the tridiagonal matrix c m is the order of the matrix c d,w,y are working arrays c isgn determines whether or not a change in sign is made c integer na integer nd integer nm1 integer nm2 c real aa(na) complex bd(nd) real x(m) complex y(m) complex d ,w , & crt ,den ,y1 ,y2 dimension a(1) ,b(1) ,c(1) , & d(1) ,w(1) , & bm1(1) ,bm2(1) ,yy(1) c do 101 j=1,m y(j) = cmplx(x(j),0.) 101 continue mm = m-1 id = nd m1 = nm1 m2 = nm2 ia = na 102 iflg = 0 if (id) 109,109,103 103 crt = bd(id) id = id-1 c c begin solution to system c d(m) = a(m)/(b(m)-crt) w(m) = y(m)/(b(m)-crt) do 104 j=2,mm k = m-j den = b(k+1)-crt-c(k+1)*d(k+2) d(k+1) = a(k+1)/den w(k+1) = (y(k+1)-c(k+1)*w(k+2))/den 104 continue den = b(1)-crt-c(1)*d(2) if (cabs(den)) 105,106,105 105 y(1) = (y(1)-c(1)*w(2))/den go to 107 106 y(1) = (1.,0.) 107 do 108 j=2,m y(j) = w(j)-d(j)*y(j-1) 108 continue 109 if (m1) 110,110,112 110 if (m2) 121,121,111 111 rt = bm2(m2) m2 = m2-1 go to 117 112 if (m2) 113,113,114 113 rt = bm1(m1) m1 = m1-1 go to 117 114 if (abs(bm1(m1))-abs(bm2(m2))) 116,116,115 115 rt = bm1(m1) m1 = m1-1 go to 117 116 rt = bm2(m2) m2 = m2-1 117 y1 = (b(1)-rt)*y(1)+c(1)*y(2) if (mm-2) 120,118,118 c c matrix multiplication c 118 do 119 j=2,mm y2 = a(j)*y(j-1)+(b(j)-rt)*y(j)+c(j)*y(j+1) y(j-1) = y1 y1 = y2 119 continue 120 y(m) = a(m)*y(m-1)+(b(m)-rt)*y(m) y(m-1) = y1 iflg = 1 go to 102 121 if (ia) 124,124,122 122 rt = aa(ia) ia = ia-1 iflg = 1 c c scalar multiplication c do 123 j=1,m y(j) = rt*y(j) 123 continue 124 if (iflg) 125,125,102 125 do 126 j=1,m yy(j) = real(y(j)) 126 continue return end subroutine cprodp (nd,bd,nm1,bm1,nm2,bm2,na,aa,x,yy,m,a,b,c,d,u,y) c c prodp applies a sequence of matrix operations to the vector x and c stores the result in yy periodic boundary conditions c and complex case c c bd,bm1,bm2 are arrays containing roots of certian b polynomials c nd,nm1,nm2 are the lengths of the arrays bd,bm1,bm2 respectively c aa array containing scalar multipliers of the vector x c na is the length of the array aa c x,yy the matrix operations are applied to x and the result is yy c a,b,c are arrays which contain the tridiagonal matrix c m is the order of the matrix c d,u,y are working arrays c isgn determines whether or not a change in sign is made c complex y ,d ,u ,v , & den ,bh ,ym ,am , & y1 ,y2 ,yh ,bd , & crt dimension a(1) ,b(1) ,c(1) ,x(1) , & y(1) ,d(1) ,u(1) ,bd(1) , & bm1(1) ,bm2(1) ,aa(1) ,yy(1) do 101 j=1,m y(j) = cmplx(x(j),0.) 101 continue mm = m-1 mm2 = m-2 id = nd m1 = nm1 m2 = nm2 ia = na 102 iflg = 0 if (id) 111,111,103 103 crt = bd(id) id = id-1 iflg = 1 c c begin solution to system c bh = b(m)-crt ym = y(m) den = b(1)-crt d(1) = c(1)/den u(1) = a(1)/den y(1) = y(1)/den v = cmplx(c(m),0.) if (mm2-2) 106,104,104 104 do 105 j=2,mm2 den = b(j)-crt-a(j)*d(j-1) d(j) = c(j)/den u(j) = -a(j)*u(j-1)/den y(j) = (y(j)-a(j)*y(j-1))/den bh = bh-v*u(j-1) ym = ym-v*y(j-1) v = -v*d(j-1) 105 continue 106 den = b(m-1)-crt-a(m-1)*d(m-2) d(m-1) = (c(m-1)-a(m-1)*u(m-2))/den y(m-1) = (y(m-1)-a(m-1)*y(m-2))/den am = a(m)-v*d(m-2) bh = bh-v*u(m-2) ym = ym-v*y(m-2) den = bh-am*d(m-1) if (cabs(den)) 107,108,107 107 y(m) = (ym-am*y(m-1))/den go to 109 108 y(m) = (1.,0.) 109 y(m-1) = y(m-1)-d(m-1)*y(m) do 110 j=2,mm k = m-j y(k) = y(k)-d(k)*y(k+1)-u(k)*y(m) 110 continue 111 if (m1) 112,112,114 112 if (m2) 123,123,113 113 rt = bm2(m2) m2 = m2-1 go to 119 114 if (m2) 115,115,116 115 rt = bm1(m1) m1 = m1-1 go to 119 116 if (abs(bm1(m1))-abs(bm2(m2))) 118,118,117 117 rt = bm1(m1) m1 = m1-1 go to 119 118 rt = bm2(m2) m2 = m2-1 c c matrix multiplication c 119 yh = y(1) y1 = (b(1)-rt)*y(1)+c(1)*y(2)+a(1)*y(m) if (mm-2) 122,120,120 120 do 121 j=2,mm y2 = a(j)*y(j-1)+(b(j)-rt)*y(j)+c(j)*y(j+1) y(j-1) = y1 y1 = y2 121 continue 122 y(m) = a(m)*y(m-1)+(b(m)-rt)*y(m)+c(m)*yh y(m-1) = y1 iflg = 1 go to 102 123 if (ia) 126,126,124 124 rt = aa(ia) ia = ia-1 iflg = 1 c c scalar multiplication c do 125 j=1,m y(j) = rt*y(j) 125 continue 126 if (iflg) 127,127,102 127 do 128 j=1,m yy(j) = real(y(j)) 128 continue return end subroutine defe4(cofx,idmn,usol,grhs) c c this subroutine first approximates the truncation error given by c trun1(x,y)=dlx**2*tx+dly**2*ty where c tx=afun(x)*uxxxx/12.0+bfun(x)*uxxx/6.0 on the interior and c at the boundaries if periodic(here uxxx,uxxxx are the third c and fourth partial derivatives of u with respect to x). c tx is of the form afun(x)/3.0*(uxxxx/4.0+uxxx/dlx) c at x=a or x=b if the boundary condition there is mixed. c tx=0.0 along specified boundaries. ty has symmetric form c in y with x,afun(x),bfun(x) replaced by y,dfun(y),efun(y). c the second order solution in usol is used to approximate c (via second order finite differencing) the trun1ation error c and the result is added to the right hand side in grhs c and then transferred to usol to be used as a new right c hand side when calling blktri for a fourth order solution. c common /spl4/ kswx ,kswy ,k ,l , & ait ,bit ,cit ,dit , & mit ,nit ,is ,ms , & js ,ns ,dlx ,dly , & tdlx3 ,tdly3 ,dlx4 ,dly4 dimension grhs(idmn,1) ,usol(idmn,1) external cofx c c c compute truncation error approximation over the entire mesh c do 30 i=is,ms xi = ait+float(i-1)*dlx call cofx (xi,ai,bi,ci) do 30 j=js,ns c c compute partial derivative approximations at (xi,yj) c call dx4(usol,idmn,i,j,uxxx,uxxxx) call dy4(usol,idmn,i,j,uyyy,uyyyy) tx = ai*uxxxx/12.0+bi*uxxx/6.0 ty=uyyyy/12.0 c c reset form of trun1ation if at boundary which is non-periodic c if (kswx.eq.1 .or. (i.gt.1 .and. i.lt.k)) go to 10 tx = ai/3.0*(uxxxx/4.0+uxxx/dlx) 10 if (kswy.eq.1 .or. (j.gt.1 .and. j.lt.l)) go to 20 ty = (uyyyy/4.0+uyyy/dly)/3.0 20 grhs(i,j)=grhs(i,j)+dly**2*(dlx**2*tx+dly**2*ty) 30 continue c c reset the right hand side in usol c do 60 i=is,ms do 50 j=js,ns usol(i,j) = grhs(i,j) 50 continue 60 continue return end subroutine defer (cofx,cofy,idmn,usol,grhs) c c this subroutine first approximates the truncation error given by c trun1(x,y)=dlx**2*tx+dly**2*ty where c tx=afun(x)*uxxxx/12.0+bfun(x)*uxxx/6.0 on the interior and c at the boundaries if periodic(here uxxx,uxxxx are the third c and fourth partial derivatives of u with respect to x). c tx is of the form afun(x)/3.0*(uxxxx/4.0+uxxx/dlx) c at x=a or x=b if the boundary condition there is mixed. c tx=0.0 along specified boundaries. ty has symmetric form c in y with x,afun(x),bfun(x) replaced by y,dfun(y),efun(y). c the second order solution in usol is used to approximate c (via second order finite differencing) the trun1ation error c and the result is added to the right hand side in grhs c and then transferred to usol to be used as a new right c hand side when calling blktri for a fourth order solution. c common /splp/ kswx ,kswy ,k ,l , & ait ,bit ,cit ,dit , & mit ,nit ,is ,ms , & js ,ns ,dlx ,dly , & tdlx3 ,tdly3 ,dlx4 ,dly4 dimension grhs(idmn,1) ,usol(idmn,1) external cofx ,cofy c c compute trun1ation error approximation over the entire mesh c do 40 j=js,ns yj = cit+float(j-1)*dly call cofy (yj,dj,ej,fj) do 30 i=is,ms xi = ait+float(i-1)*dlx call cofx (xi,ai,bi,ci) c c compute partial derivative approximations at (xi,yj) c call dx (usol,idmn,i,j,uxxx,uxxxx) call dy (usol,idmn,i,j,uyyy,uyyyy) tx = ai*uxxxx/12.0+bi*uxxx/6.0 ty = dj*uyyyy/12.0+ej*uyyy/6.0 c c reset form of trun1ation if at boundary which is non-periodic c if (kswx.eq.1 .or. (i.gt.1 .and. i.lt.k)) go to 10 tx = ai/3.0*(uxxxx/4.0+uxxx/dlx) 10 if (kswy.eq.1 .or. (j.gt.1 .and. j.lt.l)) go to 20 ty = dj/3.0*(uyyyy/4.0+uyyy/dly) 20 grhs(i,j) = grhs(i,j)+dlx**2*tx+dly**2*ty 30 continue 40 continue c c reset the right hand side in usol c do 60 i=is,ms do 50 j=js,ns usol(i,j) = grhs(i,j) 50 continue 60 continue return end subroutine dx (u,idmn,i,j,uxxx,uxxxx) c c this program computes second order finite difference c approximations to the third and fourth x c partial derivatives of u at the (i,j) mesh point c common /splp/ kswx ,kswy ,k ,l , & ait ,bit ,cit ,dit , & mit ,nit ,is ,ms , & js ,ns ,dlx ,dly , & tdlx3 ,tdly3 ,dlx4 ,dly4 dimension u(idmn,1) if (i.gt.2 .and. i.lt.(k-1)) go to 50 if (i .eq. 1) go to 10 if (i .eq. 2) go to 30 if (i .eq. k-1) go to 60 if (i .eq. k) go to 80 c c compute partial derivative approximations at x=a c 10 if (kswx .eq. 1) go to 20 uxxx = (-5.0*u(1,j)+18.0*u(2,j)-24.0*u(3,j)+14.0*u(4,j)- & 3.0*u(5,j))/(tdlx3) uxxxx = (3.0*u(1,j)-14.0*u(2,j)+26.0*u(3,j)-24.0*u(4,j)+ & 11.0*u(5,j)-2.0*u(6,j))/dlx4 return c c periodic at x=a c 20 uxxx = (-u(k-2,j)+2.0*u(k-1,j)-2.0*u(2,j)+u(3,j))/(tdlx3) uxxxx = (u(k-2,j)-4.0*u(k-1,j)+6.0*u(1,j)-4.0*u(2,j)+u(3,j))/dlx4 return c c compute partial derivative approximations at x=a+dlx c 30 if (kswx .eq. 1) go to 40 uxxx = (-3.0*u(1,j)+10.0*u(2,j)-12.0*u(3,j)+6.0*u(4,j)-u(5,j))/ & tdlx3 uxxxx = (2.0*u(1,j)-9.0*u(2,j)+16.0*u(3,j)-14.0*u(4,j)+6.0*u(5,j)- & u(6,j))/dlx4 return c c periodic at x=a+dlx c 40 uxxx = (-u(k-1,j)+2.0*u(1,j)-2.0*u(3,j)+u(4,j))/(tdlx3) uxxxx = (u(k-1,j)-4.0*u(1,j)+6.0*u(2,j)-4.0*u(3,j)+u(4,j))/dlx4 return c c compute partial derivative approximations on the interior c 50 continue uxxx = (-u(i-2,j)+2.0*u(i-1,j)-2.0*u(i+1,j)+u(i+2,j))/tdlx3 uxxxx = (u(i-2,j)-4.0*u(i-1,j)+6.0*u(i,j)-4.0*u(i+1,j)+u(i+2,j))/ & dlx4 return c c compute partial derivative approximations at x=b-dlx c 60 if (kswx .eq. 1) go to 70 uxxx = (u(k-4,j)-6.0*u(k-3,j)+12.0*u(k-2,j)-10.0*u(k-1,j)+ & 3.0*u(k,j))/tdlx3 uxxxx = (-u(k-5,j)+6.0*u(k-4,j)-14.0*u(k-3,j)+16.0*u(k-2,j)- & 9.0*u(k-1,j)+2.0*u(k,j))/dlx4 return c c periodic at x=b-dlx c 70 uxxx = (-u(k-3,j)+2.0*u(k-2,j)-2.0*u(1,j)+u(2,j))/tdlx3 uxxxx = (u(k-3,j)-4.0*u(k-2,j)+6.0*u(k-1,j)-4.0*u(1,j)+u(2,j))/ & dlx4 return c c compute partial derivative approximations at x=b c 80 uxxx = -(3.0*u(k-4,j)-14.0*u(k-3,j)+24.0*u(k-2,j)-18.0*u(k-1,j)+ & 5.0*u(k,j))/tdlx3 uxxxx = (-2.0*u(k-5,j)+11.0*u(k-4,j)-24.0*u(k-3,j)+26.0*u(k-2,j)- & 14.0*u(k-1,j)+3.0*u(k,j))/dlx4 return end subroutine dx4(u,idmn,i,j,uxxx,uxxxx) c c this program computes second order finite difference c approximations to the third and fourth x c partial derivatives of u at the (i,j) mesh point c common /spl4/ kswx ,kswy ,k ,l , & ait ,bit ,cit ,dit , & mit ,nit ,is ,ms , & js ,ns ,dlx ,dly , & tdlx3 ,tdly3 ,dlx4 ,dly4 dimension u(idmn,1) if (i.gt.2 .and. i.lt.(k-1)) go to 50 if (i .eq. 1) go to 10 if (i .eq. 2) go to 30 if (i .eq. k-1) go to 60 if (i .eq. k) go to 80 c c compute partial derivative approximations at x=a c 10 if (kswx .eq. 1) go to 20 uxxx = (-5.0*u(1,j)+18.0*u(2,j)-24.0*u(3,j)+14.0*u(4,j)- & 3.0*u(5,j))/(tdlx3) uxxxx = (3.0*u(1,j)-14.0*u(2,j)+26.0*u(3,j)-24.0*u(4,j)+ & 11.0*u(5,j)-2.0*u(6,j))/dlx4 return c c periodic at x=a c 20 uxxx = (-u(k-2,j)+2.0*u(k-1,j)-2.0*u(2,j)+u(3,j))/(tdlx3) uxxxx = (u(k-2,j)-4.0*u(k-1,j)+6.0*u(1,j)-4.0*u(2,j)+u(3,j))/dlx4 return c c compute partial derivative approximations at x=a+dlx c 30 if (kswx .eq. 1) go to 40 uxxx = (-3.0*u(1,j)+10.0*u(2,j)-12.0*u(3,j)+6.0*u(4,j)-u(5,j))/ & tdlx3 uxxxx = (2.0*u(1,j)-9.0*u(2,j)+16.0*u(3,j)-14.0*u(4,j)+6.0*u(5,j)- & u(6,j))/dlx4 return c c periodic at x=a+dlx c 40 uxxx = (-u(k-1,j)+2.0*u(1,j)-2.0*u(3,j)+u(4,j))/(tdlx3) uxxxx = (u(k-1,j)-4.0*u(1,j)+6.0*u(2,j)-4.0*u(3,j)+u(4,j))/dlx4 return c c compute partial derivative approximations on the interior c 50 continue uxxx = (-u(i-2,j)+2.0*u(i-1,j)-2.0*u(i+1,j)+u(i+2,j))/tdlx3 uxxxx = (u(i-2,j)-4.0*u(i-1,j)+6.0*u(i,j)-4.0*u(i+1,j)+u(i+2,j))/ & dlx4 return c c compute partial derivative approximations at x=b-dlx c 60 if (kswx .eq. 1) go to 70 uxxx = (u(k-4,j)-6.0*u(k-3,j)+12.0*u(k-2,j)-10.0*u(k-1,j)+ & 3.0*u(k,j))/tdlx3 uxxxx = (-u(k-5,j)+6.0*u(k-4,j)-14.0*u(k-3,j)+16.0*u(k-2,j)- & 9.0*u(k-1,j)+2.0*u(k,j))/dlx4 return c c periodic at x=b-dlx c 70 uxxx = (-u(k-3,j)+2.0*u(k-2,j)-2.0*u(1,j)+u(2,j))/tdlx3 uxxxx = (u(k-3,j)-4.0*u(k-2,j)+6.0*u(k-1,j)-4.0*u(1,j)+u(2,j))/ & dlx4 return c c compute partial derivative approximations at x=b c 80 uxxx = -(3.0*u(k-4,j)-14.0*u(k-3,j)+24.0*u(k-2,j)-18.0*u(k-1,j)+ & 5.0*u(k,j))/tdlx3 uxxxx = (-2.0*u(k-5,j)+11.0*u(k-4,j)-24.0*u(k-3,j)+26.0*u(k-2,j)- & 14.0*u(k-1,j)+3.0*u(k,j))/dlx4 return end subroutine dy (u,idmn,i,j,uyyy,uyyyy) c c this program computes second order finite difference c approximations to the third and fourth y c partial derivatives of u at the (i,j) mesh point c common /splp/ kswx ,kswy ,k ,l , & ait ,bit ,cit ,dit , & mit ,nit ,is ,ms , & js ,ns ,dlx ,dly , & tdlx3 ,tdly3 ,dlx4 ,dly4 dimension u(idmn,1) if (j.gt.2 .and. j.lt.(l-1)) go to 50 if (j .eq. 1) go to 10 if (j .eq. 2) go to 30 if (j .eq. l-1) go to 60 if (j .eq. l) go to 80 c c compute partial derivative approximations at y=c c 10 if (kswy .eq. 1) go to 20 uyyy = (-5.0*u(i,1)+18.0*u(i,2)-24.0*u(i,3)+14.0*u(i,4)- & 3.0*u(i,5))/tdly3 uyyyy = (3.0*u(i,1)-14.0*u(i,2)+26.0*u(i,3)-24.0*u(i,4)+ & 11.0*u(i,5)-2.0*u(i,6))/dly4 return c c periodic at x=a c 20 uyyy = (-u(i,l-2)+2.0*u(i,l-1)-2.0*u(i,2)+u(i,3))/tdly3 uyyyy = (u(i,l-2)-4.0*u(i,l-1)+6.0*u(i,1)-4.0*u(i,2)+u(i,3))/dly4 return c c compute partial derivative approximations at y=c+dly c 30 if (kswy .eq. 1) go to 40 uyyy = (-3.0*u(i,1)+10.0*u(i,2)-12.0*u(i,3)+6.0*u(i,4)-u(i,5))/ & tdly3 uyyyy = (2.0*u(i,1)-9.0*u(i,2)+16.0*u(i,3)-14.0*u(i,4)+6.0*u(i,5)- & u(i,6))/dly4 return c c periodic at y=c+dly c 40 uyyy = (-u(i,l-1)+2.0*u(i,1)-2.0*u(i,3)+u(i,4))/tdly3 uyyyy = (u(i,l-1)-4.0*u(i,1)+6.0*u(i,2)-4.0*u(i,3)+u(i,4))/dly4 return c c compute partial derivative approximations on the interior c 50 continue uyyy = (-u(i,j-2)+2.0*u(i,j-1)-2.0*u(i,j+1)+u(i,j+2))/tdly3 uyyyy = (u(i,j-2)-4.0*u(i,j-1)+6.0*u(i,j)-4.0*u(i,j+1)+u(i,j+2))/ & dly4 return c c compute partial derivative approximations at y=d-dly c 60 if (kswy .eq. 1) go to 70 uyyy = (u(i,l-4)-6.0*u(i,l-3)+12.0*u(i,l-2)-10.0*u(i,l-1)+ & 3.0*u(i,l))/tdly3 uyyyy = (-u(i,l-5)+6.0*u(i,l-4)-14.0*u(i,l-3)+16.0*u(i,l-2)- & 9.0*u(i,l-1)+2.0*u(i,l))/dly4 return c c periodic at y=d-dly c 70 continue uyyy = (-u(i,l-3)+2.0*u(i,l-2)-2.0*u(i,1)+u(i,2))/tdly3 uyyyy = (u(i,l-3)-4.0*u(i,l-2)+6.0*u(i,l-1)-4.0*u(i,1)+u(i,2))/ & dly4 return c c compute partial derivative approximations at y=d c 80 uyyy = -(3.0*u(i,l-4)-14.0*u(i,l-3)+24.0*u(i,l-2)-18.0*u(i,l-1)+ & 5.0*u(i,l))/tdly3 uyyyy = (-2.0*u(i,l-5)+11.0*u(i,l-4)-24.0*u(i,l-3)+26.0*u(i,l-2)- & 14.0*u(i,l-1)+3.0*u(i,l))/dly4 return end subroutine dy4(u,idmn,i,j,uyyy,uyyyy) c c this program computes second order finite difference c approximations to the third and fourth y c partial derivatives of u at the (i,j) mesh point c common /spl4/ kswx ,kswy ,k ,l , & ait ,bit ,cit ,dit , & mit ,nit ,is ,ms , & js ,ns ,dlx ,dly , & tdlx3 ,tdly3 ,dlx4 ,dly4 dimension u(idmn,1) if (j.gt.2 .and. j.lt.(l-1)) go to 50 if (j .eq. 1) go to 10 if (j .eq. 2) go to 30 if (j .eq. l-1) go to 60 if (j .eq. l) go to 80 c c compute partial derivative approximations at y=c c 10 if (kswy .eq. 1) go to 20 uyyy = (-5.0*u(i,1)+18.0*u(i,2)-24.0*u(i,3)+14.0*u(i,4)- & 3.0*u(i,5))/tdly3 uyyyy = (3.0*u(i,1)-14.0*u(i,2)+26.0*u(i,3)-24.0*u(i,4)+ & 11.0*u(i,5)-2.0*u(i,6))/dly4 return c c periodic at x=a c 20 uyyy = (-u(i,l-2)+2.0*u(i,l-1)-2.0*u(i,2)+u(i,3))/tdly3 uyyyy = (u(i,l-2)-4.0*u(i,l-1)+6.0*u(i,1)-4.0*u(i,2)+u(i,3))/dly4 return c c compute partial derivative approximations at y=c+dly c 30 if (kswy .eq. 1) go to 40 uyyy = (-3.0*u(i,1)+10.0*u(i,2)-12.0*u(i,3)+6.0*u(i,4)-u(i,5))/ & tdly3 uyyyy = (2.0*u(i,1)-9.0*u(i,2)+16.0*u(i,3)-14.0*u(i,4)+6.0*u(i,5)- & u(i,6))/dly4 return c c periodic at y=c+dly c 40 uyyy = (-u(i,l-1)+2.0*u(i,1)-2.0*u(i,3)+u(i,4))/tdly3 uyyyy = (u(i,l-1)-4.0*u(i,1)+6.0*u(i,2)-4.0*u(i,3)+u(i,4))/dly4 return c c compute partial derivative approximations on the interior c 50 continue uyyy = (-u(i,j-2)+2.0*u(i,j-1)-2.0*u(i,j+1)+u(i,j+2))/tdly3 uyyyy = (u(i,j-2)-4.0*u(i,j-1)+6.0*u(i,j)-4.0*u(i,j+1)+u(i,j+2))/ & dly4 return c c compute partial derivative approximations at y=d-dly c 60 if (kswy .eq. 1) go to 70 uyyy = (u(i,l-4)-6.0*u(i,l-3)+12.0*u(i,l-2)-10.0*u(i,l-1)+ & 3.0*u(i,l))/tdly3 uyyyy = (-u(i,l-5)+6.0*u(i,l-4)-14.0*u(i,l-3)+16.0*u(i,l-2)- & 9.0*u(i,l-1)+2.0*u(i,l))/dly4 return c c periodic at y=d-dly c 70 continue uyyy = (-u(i,l-3)+2.0*u(i,l-2)-2.0*u(i,1)+u(i,2))/tdly3 uyyyy = (u(i,l-3)-4.0*u(i,l-2)+6.0*u(i,l-1)-4.0*u(i,1)+u(i,2))/ & dly4 return c c compute partial derivative approximations at y=d c 80 uyyy = -(3.0*u(i,l-4)-14.0*u(i,l-3)+24.0*u(i,l-2)-18.0*u(i,l-1)+ & 5.0*u(i,l))/tdly3 uyyyy = (-2.0*u(i,l-5)+11.0*u(i,l-4)-24.0*u(i,l-3)+26.0*u(i,l-2)- & 14.0*u(i,l-1)+3.0*u(i,l))/dly4 return end function epmach () c c this program computes an approximate machiine epsilon (accuracy) c common /value/ v eps = 1. 101 eps = eps/10. call store (eps+1.) if (v-1.) 102,102,101 102 epmach = 100.*eps return end subroutine genbun (nperod,n,mperod,m,a,b,c,idimy,y,ierror,w) c c c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c * * c * f i s h p a k * c * * c * * c * a package of fortran subprograms for the solution of * c * * c * separable elliptic partial differential equations * c * * c * (version 3.1 , october 1980) * c * * c * by * c * * c * john adams, paul swarztrauber and roland sweet * c * * c * of * c * * c * the national center for atmospheric research * c * * c * boulder, colorado (80307) u.s.a. * c * * c * which is sponsored by * c * * c * the national science foundation * c * * c * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * c c c * * * * * * * * * purpose * * * * * * * * * * * * * * * * * * c c c subroutine genbun solves the linear system of equations c c a(i)*x(i-1,j) + b(i)*x(i,j) + c(i)*x(i+1,j) c c + x(i,j-1) - 2.*x(i,j) + x(i,j+1) = y(i,j) c c for i = 1,2,...,m and j = 1,2,...,n. c c the indices i+1 and i-1 are evaluated modulo m, i.e., c x(0,j) = x(m,j) and x(m+1,j) = x(1,j), and x(i,0) may be equal to c 0, x(i,2), or x(i,n) and x(i,n+1) may be equal to 0, x(i,n-1), or c x(i,1) depending on an input parameter. c c c * * * * * * * * parameter description * * * * * * * * * * c c * * * * * * on input * * * * * * c c nperod c indicates the values that x(i,0) and x(i,n+1) are assumed to c have. c c = 0 if x(i,0) = x(i,n) and x(i,n+1) = x(i,1). c = 1 if x(i,0) = x(i,n+1) = 0 . c = 2 if x(i,0) = 0 and x(i,n+1) = x(i,n-1). c = 3 if x(i,0) = x(i,2) and x(i,n+1) = x(i,n-1). c = 4 if x(i,0) = x(i,2) and x(i,n+1) = 0. c c n c the number of unknowns in the j-direction. n must be greater c than 2. c c mperod c = 0 if a(1) and c(m) are not zero c = 1 if a(1) = c(m) = 0 c c m c the number of unknowns in the i-direction. m must be greater c than 2. c c a,b,c c one-dimensional arrays of length m that specify the c coefficients in the linear equations given above. if mperod = 0 c the array elements must not depend upon the index i, but must be c constant. specifically, the subroutine checks the following c condition c c a(i) = c(1) c c(i) = c(1) c b(i) = b(1) c c for i=1,2,...,m. c c idimy c the row (or first) dimension of the two-dimensional array y as c it appears in the program calling genbun. this parameter is c used to specify the variable dime