#include c c purpose mexsepeli solves for the 2nd order finite c difference approximation to the separable c elliptic equation arising from conformal c mapping and grid generation. c c c c c The calling sequence from MATLAB should be c c >> [x, y, perturb, ierror] = mexsepeli ( x, y, ... c l2, m2, ... c seta, sxi, ... c nwrk, nx2 ); c c c Output Parameters c x,y: c represent "x" and "y" grid coordinates c perturb: c Optional. Don't think this is needed for our c application. c ierror: c Optional. An error flag that indicates invalid c input parameters or failure to find a solution. c = 0 no error c = 4 if attempt to find a solution fails. c (the linear system generated is not c diagonally dominant.) c c Input Parameters c x,y: c represent "x" and "y" grid coordinates c l2: c x independent variable ranges from 0 to l2 c m2: c y independent variable ranges from 0 to m2 c seta; c digitized points on boundaries 1,3 c sxi: c digitized points on boundaries 2,4 c c c So c 1) nlhs = 2 c 2) plhs(1) = x c plhs(2) = y c 3) nrhs = 6 c 4) prhs(1) = x c prhs(2) = y c prhs(3) = l2 c prhs(4) = m2 c prhs(5) = seta c prhs(6) = sxi c c c $Id: mexsepeli.f,v 1.2 1997/10/03 19:43:04 jevans Exp jevans $ c Currently locked by $Locker: jevans $ (not locked if blank) c (Time in GMT, EST=GMT-5:00 c $Log: mexsepeli.f,v $ c Revision 1.2 1997/10/03 19:43:04 jevans c Had to change compile time options to include "-align dcommons" c and "-r8". c c the "mexCopyReal8ToPtr" and related functions suffer from the c fact that they assume all arrays to be a single column. Since c my x and y are dimensioned to be larger than their logical c size (i.e. they are 800x800, but might logically be only 100x80 c or so), those arrays had to be columnified before passing them c back to MATLAB, and "uncolumnified" upon passing them into c MATLAB. c c Revision 1.1 1997/10/03 17:50:43 jevans c Initial revision c c c subroutine mexFunction ( nlhs, plhs, nrhs, prhs ) POINTER plhs(*), prhs(*) integer nlhs, nrhs include 'cgridgen.inc' c c These are the number of rows and columns (m by n) of the c matrices we work with. integer size_m, size_n, size_mn integer i, j, k c c Need some temporary space to store the double-precision c matlab arguments. The correct number of rows and columns c is not known at first. real*8 tempx(0:nx2,0:ny2), tempy(0:nx2,0:ny2) c c input ranges for independent variables. c 0 < x < l2, 0 < y < y2 integer l2, m2 c c Pointers to input arguments POINTER p real*8 pa(1) c c Check for proper number of arguments. c There must be at least 1 argument on the left hand size and c no more than 3. if ( nlhs .lt. 1 ) then call mexErrMsgTxt ( 'Not enough input args for mexsepeli.' ) elseif ( nlhs .gt. 3 ) then call mexErrMsgTxt ( 'Too many input args for mexsepeli.' ) endif if ( nrhs .ne. 6 ) then call mexErrMsgTxt('mexsepeli requires 6 input arguments') endif c c Check that all input arguments were numeric. if ( mxIsNumeric(prhs(1)) .ne. 1 ) then call mexErrMsgTxt + ( 'x input array should be numeric' ) endif if ( mxIsNumeric(prhs(2)) .ne. 1 ) then call mexErrMsgTxt + ( 'y input array should be numeric' ) endif if ( mxIsNumeric(prhs(3)) .ne. 1 ) then call mexErrMsgTxt + ( 'l2 input should be numeric' ) endif if ( mxIsNumeric(prhs(4)) .ne. 1 ) then call mexErrMsgTxt + ( 'm2 input should be numeric' ) endif if ( mxIsNumeric(prhs(5)) .ne. 1 ) then call mexErrMsgTxt + ( 'seta input array should be numeric' ) endif if ( mxIsNumeric(prhs(6)) .ne. 1 ) then call mexErrMsgTxt + ( 'sxi input array should be numeric' ) endif c c Check to see that the l2 and m2 arguments were scalars size_m = mxGetM ( prhs(3) ) size_n = mxGetN ( prhs(3) ) if ( size_n .ne. 1 .or. size_m .ne. 1 ) then call mexErrMsgTxt + ('l2 argument should be a scalar' ) endif size_m = mxGetM ( prhs(4) ) size_n = mxGetN ( prhs(4) ) if ( size_n .ne. 1 .or. size_m .ne. 1 ) then call mexErrMsgTxt + ('m2 argument should be a scalar' ) endif c c Load the data into Fortran matrices. c c x p = mxGetPr(prhs(1)) size_m = mxGetM ( prhs(1) ) size_n = mxGetN ( prhs(1) ) size_mn = size_m*size_n call mxCopyPtrToReal8(p,tempx,size_mn) k = 0 do 10 j = 0, size_n - 1 do 5 i = 0, size_m - 1 x(i,j) = tempx(k,0) 5 k = k + 1 10 continue continue c c y p = mxGetPr(prhs(2)) size_m = mxGetM ( prhs(2) ) size_n = mxGetN ( prhs(2) ) size_mn = size_m*size_n call mxCopyPtrToReal8(p,tempy,size_mn) k = 0 do 20 j = 0, size_n - 1 do 15 i = 0, size_m - 1 y(i,j) = tempy(k,0) 15 k = k + 1 20 continue continue c c l2 p = mxGetPr(prhs(3)) call mxCopyPtrToReal8(p,pa,1) l2 = pa(1) c c m2 p = mxGetPr(prhs(4)) call mxCopyPtrToReal8(p,pa,1) m2 = pa(1) c c seta p = mxGetPr(prhs(5)) size_m = mxGetM ( prhs(5) ) size_n = mxGetN ( prhs(5) ) neta = size_m*size_n call mxCopyPtrToReal8(p,seta,neta) c c sxi p = mxGetPr(prhs(6)) size_m = mxGetM ( prhs(6) ) size_n = mxGetN ( prhs(6) ) nxi = size_m*size_n call mxCopyPtrToReal8(p,sxi,nxi) c c Construct right hand side do 30 j = 0, neta - 1 do 25 i = 0, nxi - 1 rhs(i,j) = 0. 25 continue 30 continue c c Call the computational subroutine. ewrk(1) = nwrk call sepeli(0,2,0.,float(l2),l2,1,wrk,wrk,wrk,wrk,0.,float(m2),m2, * 1,wrk,wrk,wrk,wrk,rhs,x,nx2+1,ewrk,pertrb,ierr) if ( ierr .eq. 1 ) then call mexErrMsgTxt + ( 'range of independent variables is out of whack' ) else if ( ierr .eq. 2 ) then call mexErrMsgTxt + ( 'boundary condition mbdcnd wrongly specified' ) else if ( ierr .eq. 3 ) then call mexErrMsgTxt + ( 'boundary condition nbdcnd wrongly specified' ) else if ( ierr .eq. 4 ) then call mexErrMsgTxt + ( 'linear system generated is not diagonally dominant' ) else if ( ierr .eq. 5 ) then call mexErrMsgTxt + ( 'idmn is too small' ) else if ( ierr .eq. 6 ) then call mexErrMsgTxt + ( 'm is too small or too large' ) else if ( ierr .eq. 7 ) then call mexErrMsgTxt + ( 'n is too small or too large' ) else if ( ierr .eq. 8 ) then call mexErrMsgTxt + ( 'iorder is not 2 or 4' ) else if ( ierr .eq. 9 ) then call mexErrMsgTxt + ( 'intl is not 0 or 1' ) else if ( ierr .eq. 10 ) then call mexErrMsgTxt + ( 'afun*dfun less than or equal to 0' ) else if ( ierr .eq. 11 ) then call mexErrMsgTxt + ( 'work space length input in w(1) is not right' ) end if ewrk(1) = nwrk call sepeli(0,2,0.,float(l2),l2,1,wrk,wrk,wrk,wrk,0.,float(m2),m2, * 1,wrk,wrk,wrk,wrk,rhs,y,nx2+1,ewrk,pertrb,ierr) if ( ierr .eq. 1 ) then call mexErrMsgTxt + ( 'range of independent variables is out of whack' ) else if ( ierr .eq. 2 ) then call mexErrMsgTxt + ( 'boundary condition mbdcnd wrongly specified' ) else if ( ierr .eq. 3 ) then call mexErrMsgTxt + ( 'boundary condition nbdcnd wrongly specified' ) else if ( ierr .eq. 4 ) then call mexErrMsgTxt + ( 'linear system generated is not diagonally dominant' ) else if ( ierr .eq. 5 ) then call mexErrMsgTxt + ( 'idmn is too small' ) else if ( ierr .eq. 6 ) then call mexErrMsgTxt + ( 'm is too small or too large' ) else if ( ierr .eq. 7 ) then call mexErrMsgTxt + ( 'n is too small or too large' ) else if ( ierr .eq. 8 ) then call mexErrMsgTxt + ( 'iorder is not 2 or 4' ) else if ( ierr .eq. 9 ) then call mexErrMsgTxt + ( 'intl is not 0 or 1' ) else if ( ierr .eq. 10 ) then call mexErrMsgTxt + ( 'afun*dfun less than or equal to 0' ) else if ( ierr .eq. 11 ) then call mexErrMsgTxt + ( 'work space length input in w(1) is not right' ) end if c c Create matrices for output arguments c Since the mxCopyReal8ToPtr routine just assumes one long c column matrix, we have to "columnify" both x and y. c c x size_m = mxGetM ( prhs(1) ) size_n = mxGetN ( prhs(1) ) size_mn = size_m*size_n plhs(1) = mxCreateFull(size_m, size_n, 0) p = mxGetPr(plhs(1)) k = 0 do 40 j = 0, size_n - 1 do 35 i = 0, size_m - 1 tempx(k,0) = x(i,j) k = k + 1 35 continue 40 continue call mxCopyReal8ToPtr ( tempx, p, size_mn ) c c y size_m = mxGetM ( prhs(2) ) size_n = mxGetN ( prhs(2) ) size_mn = size_m*size_n plhs(2) = mxCreateFull(size_m, size_n, 0) p = mxGetPr(plhs(2)) k = 0 do 50 j = 0, size_n - 1 do 45 i = 0, size_m - 1 tempy(k,0) = y(i,j) k = k + 1 45 continue 50 continue call mxCopyReal8ToPtr ( tempy, p, size_mn ) return end c Subroutine SEPELI(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D, c * N,NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,GRHS,USOL,IDMN,W,PERTRB, c * IERROR) Subroutine SEPELI(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D, * N,NBDCND,BDC,GAMA,BDD,XNU,GRHS,USOL,IDMN,W,PERTRB, * IERROR) c c dimension of bda(n+1), bdb(n+1), bdc(m+1), bdd(m+1), c arguments usol(idmn,n+1),grhs(idmn,n+1), c w (see argument list) c c latest revision march 1985 c c purpose sepeli solves for either the second-order c finite difference approximation or a c fourth-order approximation to a separable c elliptic equation c c 2 2 c af(x)*d u/dx + bf(x)*du/dx + cf(x)*u + c 2 2 c df(y)*d u/dy + ef(y)*du/dy + ff(y)*u c c = g(x,y) c c on a rectangle (x greater than or equal to a c and less than or equal to b; y greater than c or equal to c and less than or equal to d). c any combination of periodic or mixed boundary c conditions is allowed. c c the possible boundary conditions are: c in the x-direction: c (0) periodic, u(x+b-a,y)=u(x,y) for all c y,x (1) u(a,y), u(b,y) are specified for c all y c (2) u(a,y), du(b,y)/dx+beta*u(b,y) are c specified for all y c (3) du(a,y)/dx+alpha*u(a,y),du(b,y)/dx+ c beta*u(b,y) are specified for all y c (4) du(a,y)/dx+alpha*u(a,y),u(b,y) are c specified for all y c c in the y-direction: c (0) periodic, u(x,y+d-c)=u(x,y) for all x,y c (1) u(x,c),u(x,d) are specified for all x c (2) u(x,c),du(x,d)/dy+xnu*u(x,d) are c specified for all x c (3) du(x,c)/dy+gama*u(x,c),du(x,d)/dy+ c xnu*u(x,d) are specified for all x c (4) du(x,c)/dy+gama*u(x,c),u(x,d) are c specified for all x c c usage call sepeli (intl,iorder,a,b,m,mbdcnd,bda, c alpha,bdb,beta,c,d,n,nbdcnd,bdc, c gama,bdd,xnu,cofx,cofy,grhs,usol, c idmn,w,pertrb,ierror) c c arguments c on input intl c = 0 on initial entry to sepeli or if any c of the arguments c,d, n, nbdcnd, cofy c are changed from a previous call c = 1 if c, d, n, nbdcnd, cofy are unchanged c from the previous call. c c iorder c = 2 if a second-order approximation c is sought c = 4 if a fourth-order approximation c is sought c c a,b c the range of the x-independent variable, c i.e., x is greater than or equal to a c and less than or equal to b. a must be c less than b. c c m c the number of panels into which the c interval [a,b] is subdivided. hence, c there will be m+1 grid points in the x- c direction given by xi=a+(i-1)*dlx c for i=1,2,...,m+1 where dlx=(b-a)/m is c the panel width. m must be less than c idmn and greater than 5. c c mbdcnd c indicates the type of boundary condition c at x=a and x=b c c = 0 if the solution is periodic in x, i.e., c u(x+b-a,y)=u(x,y) for all y,x c = 1 if the solution is specified at x=a c and x=b, i.e., u(a,y) and u(b,y) are c specified for all y c = 2 if the solution is specified at x=a and c the boundary condition is mixed at x=b, c i.e., u(a,y) and du(b,y)/dx+beta*u(b,y) c are specified for all y c = 3 if the boundary conditions at x=a and c x=b are mixed, i.e., c du(a,y)/dx+alpha*u(a,y) and c du(b,y)/dx+beta*u(b,y) are specified c for all y c = 4 if the boundary condition at x=a is c mixed and the solution is specified c at x=b, i.e., du(a,y)/dx+alpha*u(a,y) c and u(b,y) are specified for all y c c bda c a one-dimensional array of length n+1 c that specifies the values of c du(a,y)/dx+ alpha*u(a,y) at x=a, when c mbdcnd=3 or 4. c bda(j) = du(a,yj)/dx+alpha*u(a,yj), c j=1,2,...,n+1. when mbdcnd has any other c other value, bda is a dummy parameter. c c alpha c the scalar multiplying the solution in c case of a mixed boundary condition at x=a c (see argument bda). if mbdcnd is not c equal to 3 or 4 then alpha is a dummy c parameter. c c bdb c a one-dimensional array of length n+1 c that specifies the values of c du(b,y)/dx+ beta*u(b,y) at x=b. c when mbdcnd=2 or 3 c bdb(j) = du(b,yj)/dx+beta*u(b,yj), c j=1,2,...,n+1. when mbdcnd has any other c other value, bdb is a dummy parameter. c c beta c the scalar multiplying the solution in c case of a mixed boundary condition at c x=b (see argument bdb). if mbdcnd is c not equal to 2 or 3 then beta is a dummy c parameter. c c c,d c the range of the y-independent variable, c i.e., y is greater than or equal to c c and less than or equal to d. c must be c less than d. c c n c the number of panels into which the c interval [c,d] is subdivided. c hence, there will be n+1 grid points c in the y-direction given by c yj=c+(j-1)*dly for j=1,2,...,n+1 where c dly=(d-c)/n is the panel width. c in addition, n must be greater than 4. c c nbdcnd c indicates the types of boundary conditions c at y=c and y=d c c = 0 if the solution is periodic in y, c i.e., u(x,y+d-c)=u(x,y) for all x,y c = 1 if the solution is specified at y=c c and y = d, i.e., u(x,c) and u(x,d) c are specified for all x c = 2 if the solution is specified at y=c c and the boundary condition is mixed c at y=d, i.e., u(x,c) and c du(x,d)/dy+xnu*u(x,d) are specified c for all x c = 3 if the boundary conditions are mixed c at y=c and y=d, i.e., c du(x,d)/dy+gama*u(x,c) and c du(x,d)/dy+xnu*u(x,d) are specified c for all x c = 4 if the boundary condition is mixed c at y=c and the solution is specified c at y=d, i.e. du(x,c)/dy+gama*u(x,c) c and u(x,d) are specified for all x c c bdc c a one-dimensional array of length m+1 c that specifies the value of c du(x,c)/dy+gama*u(x,c) at y=c. c when nbdcnd=3 or 4 bdc(i) = du(xi,c)/dy + c gama*u(xi,c), i=1,2,...,m+1. c when nbdcnd has any other value, bdc c is a dummy parameter. c c gama c the scalar multiplying the solution in c case of a mixed boundary condition at c y=c (see argument bdc). if nbdcnd is c not equal to 3 or 4 then gama is a dummy c parameter. c c bdd c a one-dimensional array of length m+1 c that specifies the value of c du(x,d)/dy + xnu*u(x,d) at y=c. c when nbdcnd=2 or 3 bdd(i) = du(xi,d)/dy + c xnu*u(xi,d), i=1,2,...,m+1. c when nbdcnd has any other value, bdd c is a dummy parameter. c c xnu c the scalar multiplying the solution in c case of a mixed boundary condition at c y=d (see argument bdd). if nbdcnd is c not equal to 2 or 3 then xnu is a c dummy parameter. c c cofx c a user-supplied subprogram with c parameters x, afun, bfun, cfun which c returns the values of the x-dependent c coefficients af(x), bf(x), cf(x) in the c elliptic equation at x. c c cofy c a user-supplied subprogram with parameters c y, dfun, efun, ffun which returns the c values of the y-dependent coefficients c df(y), ef(y), ff(y) in the elliptic c equation at y. c c note: cofx and cofy must be declared c external in the calling routine. c the values returned in afun and dfun c must satisfy afun*dfun greater than 0 c for a less than x less than b, c less c than y less than d (see ierror=10). c the coefficients provided may lead to a c matrix equation which is not diagonally c dominant in which case solution may fail c (see ierror=4). c c grhs c a two-dimensional array that specifies the c values of the right-hand side of the c elliptic equation, i.e., c grhs(i,j)=g(xi,yi), for i=2,...,m, c j=2,...,n. at the boundaries, grhs is c defined by c c mbdcnd grhs(1,j) grhs(m+1,j) c ------ --------- ----------- c 0 g(a,yj) g(b,yj) c 1 * * c 2 * g(b,yj) j=1,2,...,n+1 c 3 g(a,yj) g(b,yj) c 4 g(a,yj) * c c nbdcnd grhs(i,1) grhs(i,n+1) c ------ --------- ----------- c 0 g(xi,c) g(xi,d) c 1 * * c 2 * g(xi,d) i=1,2,...,m+1 c 3 g(xi,c) g(xi,d) c 4 g(xi,c) * c c where * means these quantities are not used. c grhs should be dimensioned idmn by at least c n+1 in the calling routine. c c usol c a two-dimensional array that specifies the c values of the solution along the boundaries. c at the boundaries, usol is defined by c c mbdcnd usol(1,j) usol(m+1,j) c ------ --------- ----------- c 0 * * c 1 u(a,yj) u(b,yj) c 2 u(a,yj) * j=1,2,...,n+1 c 3 * * c 4 * u(b,yj) c c nbdcnd usol(i,1) usol(i,n+1) c ------ --------- ----------- c 0 * * c 1 u(xi,c) u(xi,d) c 2 u(xi,c) * i=1,2,...,m+1 c 3 * * c 4 * u(xi,d) c c where * means the quantities are not used c in the solution. c c if iorder=2, the user may equivalence grhs c and usol to save space. note that in this c case the tables specifying the boundaries c of the grhs and usol arrays determine the c boundaries uniquely except at the corners. c if the tables call for both g(x,y) and c u(x,y) at a corner then the solution must c be chosen. for example, if mbdcnd=2 and c nbdcnd=4, then u(a,c), u(a,d), u(b,d) must c be chosen at the corners in addition c to g(b,c). c c if iorder=4, then the two arrays, usol and c grhs, must be distinct. c c usol should be dimensioned idmn by at least c n+1 in the calling routine. c c idmn c the row (or first) dimension of the arrays c grhs and usol as it appears in the program c calling sepeli. this parameter is used c to specify the variable dimension of grhs c and usol. idmn must be at least 7 and c greater than or equal to m+1. c c w c a one-dimensional array that must be c provided by the user for work space. c let k=int(log2(n+1))+1 and set l=2**(k+1). c then (k-2)*l+k+10*n+12*m+27 will suffice c as a length of w. the actual length of w c in the calling routine must be set in w(1) c (see ierror=11). c c on output usol c contains the approximate solution to the c elliptic equation. c usol(i,j) is the approximation to u(xi,yj) c for i=1,2...,m+1 and j=1,2,...,n+1. c the approximation has error c o(dlx**2+dly**2) if called with iorder=2 c and o(dlx**4+dly**4) if called with c iorder=4. c c w c contains intermediate values that must not c be destroyed if sepeli is called again with c intl=1. in addition w(1) contains the c exact minimal length (in floating point) c required for the work space (see ierror=11). c c pertrb c if a combination of periodic or derivative c boundary conditions c (i.e., alpha=beta=0 if mbdcnd=3; c gama=xnu=0 if nbdcnd=3) is specified c and if the coefficients of u(x,y) in the c separable elliptic equation are zero c (i.e., cf(x)=0 for x greater than or equal c to a and less than or equal to b; c ff(y)=0 for y greater than or equal to c c and less than or equal to d) then a c solution may not exist. pertrb is a c constant calculated and subtracted from c the right-hand side of the matrix equations c generated by sepeli which insures that a c solution exists. sepeli then computes this c solution which is a weighted minimal least c squares solution to the original problem. c c ierror c an error flag that indicates invalid input c parameters or failure to find a solution c = 0 no error c = 1 if a greater than b or c greater than d c = 2 if mbdcnd less than 0 or mbdcnd greater c than 4 c = 3 if nbdcnd less than 0 or nbdcnd greater c than 4 c = 4 if attempt to find a solution fails. c (the linear system generated is not c diagonally dominant.) c = 5 if idmn is too small c (see discussion of idmn) c = 6 if m is too small or too large c (see discussion of m) c = 7 if n is too small (see discussion of n) c = 8 if iorder is not 2 or 4 c = 9 if intl is not 0 or 1 c = 10 if afun*dfun less than or equal to 0 c for some interior mesh point (xi,yj) c = 11 if the work space length input in w(1) c is less than the exact minimal work c space length required output in w(1). c c note (concerning ierror=4): for the c coefficients input through cofx, cofy, c the discretization may lead to a block c tridiagonal linear system which is not c diagonally dominant (for example, this c happens if cfun=0 and bfun/(2.*dlx) greater c than afun/dlx**2). in this case solution c may fail. this cannot happen in the limit c as dlx, dly approach zero. hence, the c condition may be remedied by taking larger c values for m or n. c c special conditions see cofx, cofy argument descriptions above. c c i/o none c c precision single c c required library blktri, comf, q8qst4, and sepaux, which are c files loaded by default on ncar's cray machines. c c language fortran c c history developed at ncar during 1975-76 by c john c. adams of the scientific computing c division. released on ncar's public software c libraries in january 1980. c c portability fortran 66 c c algorithm sepeli automatically discretizes the c separable elliptic equation which is then c solved by a generalized cyclic reduction c algorithm in the subroutine, blktri. the c fourth-order solution is obtained using c 'deferred corrections' which is described c and referenced in sections, references and c method. c c timing the operational count is proportional to c m*n*log2(n). c c accuracy the following accuracy results were obtained c on a cdc 7600. note that the fourth-order c accuracy is not realized until the mesh is c sufficiently refined. c c second-order fourth-order c m n error error c c 6 6 6.8e-1 1.2e0 c 14 14 1.4e-1 1.8e-1 c 30 30 3.2e-2 9.7e-3 c 62 62 7.5e-3 3.0e-4 c 126 126 1.8e-3 3.5e-6 c c portability fortran 66 c c references keller, h.b., numerical methods for two-point c boundary-value problems, blaisdel (1968), c waltham, mass. c c swarztrauber, p., and r. sweet (1975): c efficient fortran subprograms for the c solution of elliptic partial differential c equations. ncar technical note c ncar-tn/ia-109, pp. 135-137. c*********************************************************************** Dimension GRHS(IDMN,1), USOL(IDMN,1) Dimension BDA(1), BDB(1), BDC(1), BDD(1), W(1) c External COFX, COFY Logical Q8Q4 Save Q8Q4 Data Q8Q4 /.TRUE./ c If (Q8Q4) Then c call q8qst4('loclib','sepeli','sepeli','version 01') Q8Q4 = .FALSE. End If c c check input parameters c c Call CHKPRM(INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,COFX,COFY,IDMN, c * IERROR) Call CHKPRM(INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,IDMN, * IERROR) If (IERROR.NE.0) Return c c compute minimum work space and check work space length input c L = N + 1 If (NBDCND.EQ.0) L = N LOGB2N = INT(ALOG(FLOAT(L)+0.5)/ALOG(2.0)) + 1 LL = 2 ** (LOGB2N+1) K = M + 1 L = N + 1 LENGTH = (LOGB2N-2) * LL + LOGB2N + MAX0(2*L,6*K) + 5 If (NBDCND.EQ.0) LENGTH = LENGTH + 2 * L IERROR = 11 LINPUT = INT(W(1)+0.5) LOUTPT = LENGTH + 6 * (K+L) + 1 W(1) = FLOAT(LOUTPT) If (LOUTPT.GT.LINPUT) Return IERROR = 0 c c set work space indices c I1 = LENGTH + 2 I2 = I1 + L I3 = I2 + L I4 = I3 + L I5 = I4 + L I6 = I5 + L I7 = I6 + L I8 = I7 + K I9 = I8 + K I10 = I9 + K I11 = I10 + K I12 = I11 + K I13 = 2 c Call SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,N, c * NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,W(I1),W(I2),W(I3),W(I4), c * W(I5),W(I6),W(I7),W(I8),W(I9),W(I10),W(I11),W(I12),GRHS,USOL, c * IDMN,W(I13),PERTRB,IERROR) Call SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D,N, * NBDCND,BDC,GAMA,BDD,XNU,W(I1),W(I2),W(I3),W(I4), * W(I5),W(I6),W(I7),W(I8),W(I9),W(I10),W(I11),W(I12),GRHS,USOL, * IDMN,W(I13),PERTRB,IERROR) Return End c Subroutine SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D, c * N,NBDCND,BDC,GAMA,BDD,XNU,COFX,COFY,AN,BN,CN,DN,UN,ZN,AM,BM, c * CM,DM,UM,ZM,GRHS,USOL,IDMN,W,PERTRB,IERROR) Subroutine SPELIP(INTL,IORDER,A,B,M,MBDCND,BDA,ALPHA,BDB,BETA,C,D, * N,NBDCND,BDC,GAMA,BDD,XNU,AN,BN,CN,DN,UN,ZN,AM,BM, * CM,DM,UM,ZM,GRHS,USOL,IDMN,W,PERTRB,IERROR) c c spelip sets up vectors and arrays for input to blktri c and computes a second order solution in usol. a return jump to c sepeli occurrs if iorder=2. if iorder=4 a fourth order c solution is generated in usol. c Dimension BDA(1), BDB(1), BDC(1), BDD(1), W(1) Dimension GRHS(IDMN,1), USOL(IDMN,1) Dimension AN(1), BN(1), CN(1), DN(1), UN(1), ZN(1) Dimension AM(1), BM(1), CM(1), DM(1), UM(1), ZM(1) Common /SPLP/ KSWX, KSWY, K, L, AIT, BIT, CIT, DIT, MIT, NIT, IS, * MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4 Logical SINGLR c External COFX, COFY c c set parameters internally c KSWX = MBDCND + 1 KSWY = NBDCND + 1 K = M + 1 L = N + 1 AIT = A BIT = B CIT = C DIT = D c c set right hand side values from grhs in usol on the interior c and non-specified boundaries. c Do 20 I = 2, M Do 10 J = 2, N USOL(I,J) = GRHS(I,J) 10 Continue 20 Continue If (KSWX.NE.2.AND.KSWX.NE.3) Then Do 30 J = 2, N USOL(1,J) = GRHS(1,J) 30 Continue End If If (KSWX.NE.2.AND.KSWX.NE.5) Then Do 40 J = 2, N USOL(K,J) = GRHS(K,J) 40 Continue End If If (KSWY.NE.2.AND.KSWY.NE.3) Then Do 50 I = 2, M USOL(I,1) = GRHS(I,1) 50 Continue End If If (KSWY.NE.2.AND.KSWY.NE.5) Then Do 60 I = 2, M USOL(I,L) = GRHS(I,L) 60 Continue End If If (KSWX.NE.2.AND.KSWX.NE.3.AND.KSWY.NE.2.AND.KSWY.NE.3) * USOL(1,1) = GRHS(1,1) If (KSWX.NE.2.AND.KSWX.NE.5.AND.KSWY.NE.2.AND.KSWY.NE.3) * USOL(K,1) = GRHS(K,1) If (KSWX.NE.2.AND.KSWX.NE.3.AND.KSWY.NE.2.AND.KSWY.NE.5) * USOL(1,L) = GRHS(1,L) If (KSWX.NE.2.AND.KSWX.NE.5.AND.KSWY.NE.2.AND.KSWY.NE.5) * USOL(K,L) = GRHS(K,L) I1 = 1 c c set switches for periodic or non-periodic boundaries c MP = 1 NP = 1 If (KSWX.EQ.1) MP = 0 If (KSWY.EQ.1) NP = 0 c c set dlx,dly and size of block tri-diagonal system generated c in nint,mint c DLX = (BIT-AIT) / FLOAT(M) MIT = K - 1 If (KSWX.EQ.2) MIT = K - 2 If (KSWX.EQ.4) MIT = K DLY = (DIT-CIT) / FLOAT(N) NIT = L - 1 If (KSWY.EQ.2) NIT = L - 2 If (KSWY.EQ.4) NIT = L TDLX3 = 2.0 * DLX ** 3 DLX4 = DLX ** 4 TDLY3 = 2.0 * DLY ** 3 DLY4 = DLY ** 4 c c set subscript limits for portion of array to input to blktri c IS = 1 JS = 1 If (KSWX.EQ.2.OR.KSWX.EQ.3) IS = 2 If (KSWY.EQ.2.OR.KSWY.EQ.3) JS = 2 NS = NIT + JS - 1 MS = MIT + IS - 1 c c set x - direction c Do 70 I = 1, MIT XI = AIT + FLOAT(IS+I-2) * DLX Call COFX(XI,AI,BI,CI) AXI = (AI/DLX-0.5*BI) / DLX BXI = -2. * AI / DLX ** 2 + CI CXI = (AI/DLX+0.5*BI) / DLX AM(I) = AXI BM(I) = BXI CM(I) = CXI 70 Continue c c set y direction c Do 80 J = 1, NIT YJ = CIT + FLOAT(JS+J-2) * DLY Call COFY(YJ,DJ,EJ,FJ) DYJ = (DJ/DLY-0.5*EJ) / DLY EYJ = (-2.*DJ/DLY**2+FJ) FYJ = (DJ/DLY+0.5*EJ) / DLY AN(J) = DYJ BN(J) = EYJ CN(J) = FYJ 80 Continue c c adjust edges in x direction unless periodic c AX1 = AM(1) CXM = CM(MIT) Go To (130,90,110,120,100), KSWX c c dirichlet-dirichlet in x direction c 90 AM(1) = 0.0 CM(MIT) = 0.0 Go To 130 c c mixed-dirichlet in x direction c 100 AM(1) = 0.0 BM(1) = BM(1) + 2. * ALPHA * DLX * AX1 CM(1) = CM(1) + AX1 CM(MIT) = 0.0 Go To 130 c c dirichlet-mixed in x direction c 110 AM(1) = 0.0 AM(MIT) = AM(MIT) + CXM BM(MIT) = BM(MIT) - 2. * BETA * DLX * CXM CM(MIT) = 0.0 Go To 130 c c mixed - mixed in x direction c 120 Continue AM(1) = 0.0 BM(1) = BM(1) + 2. * DLX * ALPHA * AX1 CM(1) = CM(1) + AX1 AM(MIT) = AM(MIT) + CXM BM(MIT) = BM(MIT) - 2. * DLX * BETA * CXM CM(MIT) = 0.0 130 Continue c c adjust in y direction unless periodic c DY1 = AN(1) FYN = CN(NIT) Go To (180,140,160,170,150), KSWY c c dirichlet-dirichlet in y direction c 140 Continue AN(1) = 0.0 CN(NIT) = 0.0 Go To 180 c c mixed-dirichlet in y direction c 150 Continue AN(1) = 0.0 BN(1) = BN(1) + 2. * DLY * GAMA * DY1 CN(1) = CN(1) + DY1 CN(NIT) = 0.0 Go To 180 c c dirichlet-mixed in y direction c 160 AN(1) = 0.0 AN(NIT) = AN(NIT) + FYN BN(NIT) = BN(NIT) - 2. * DLY * XNU * FYN CN(NIT) = 0.0 Go To 180 c c mixed - mixed direction in y direction c 170 Continue AN(1) = 0.0 BN(1) = BN(1) + 2. * DLY * GAMA * DY1 CN(1) = CN(1) + DY1 AN(NIT) = AN(NIT) + FYN BN(NIT) = BN(NIT) - 2.0 * DLY * XNU * FYN CN(NIT) = 0.0 180 If (KSWX.NE.1) Then c c adjust usol along x edge c Do 190 J = JS, NS If (KSWX.EQ.2.OR.KSWX.EQ.3) Then USOL(IS,J) = USOL(IS,J) - AX1 * USOL(1,J) Else USOL(IS,J) = USOL(IS,J) + 2.0 * DLX * AX1 * BDA(J) End If If (KSWX.EQ.2.OR.KSWX.EQ.5) Then USOL(MS,J) = USOL(MS,J) - CXM * USOL(K,J) Else USOL(MS,J) = USOL(MS,J) - 2.0 * DLX * CXM * BDB(J) End If 190 Continue End If If (KSWY.NE.1) Then c c adjust usol along y edge c Do 200 I = IS, MS If (KSWY.EQ.2.OR.KSWY.EQ.3) Then USOL(I,JS) = USOL(I,JS) - DY1 * USOL(I,1) Else USOL(I,JS) = USOL(I,JS) + 2.0 * DLY * DY1 * BDC(I) End If If (KSWY.EQ.2.OR.KSWY.EQ.5) Then USOL(I,NS) = USOL(I,NS) - FYN * USOL(I,L) Else USOL(I,NS) = USOL(I,NS) - 2.0 * DLY * FYN * BDD(I) End If 200 Continue End If c c save adjusted edges in grhs if iorder=4 c If (IORDER.EQ.4) Then Do 210 J = JS, NS GRHS(IS,J) = USOL(IS,J) GRHS(MS,J) = USOL(MS,J) 210 Continue Do 220 I = IS, MS GRHS(I,JS) = USOL(I,JS) GRHS(I,NS) = USOL(I,NS) 220 Continue End If IORD = IORDER PERTRB = 0.0 c c check if operator is singular c c Call CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,COFX,COFY,SINGLR) Call CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,SINGLR) c c compute non-zero eigenvector in null space of transpose c if singular c If (SINGLR) Call SEPTRI(MIT,AM,BM,CM,DM,UM,ZM) If (SINGLR) Call SEPTRI(NIT,AN,BN,CN,DN,UN,ZN) c c make initialization call to blktri c If (INTL.EQ.0) Call BLKTRI(INTL,NP,NIT,AN,BN,CN,MP,MIT,AM,BM,CM, * IDMN,USOL(IS,JS),IERROR,W) If (IERROR.NE.0) Return c c adjust right hand side if necessary c 230 Continue If (SINGLR) Call SEPORT(USOL,IDMN,ZN,ZM,PERTRB) c c compute solution c Call BLKTRI(I1,NP,NIT,AN,BN,CN,MP,MIT,AM,BM,CM,IDMN,USOL(IS,JS), * IERROR,W) If (IERROR.NE.0) Return c c set periodic boundaries if necessary c If (KSWX.EQ.1) Then Do 240 J = 1, L USOL(K,J) = USOL(1,J) 240 Continue End If If (KSWY.EQ.1) Then Do 250 I = 1, K USOL(I,L) = USOL(I,1) 250 Continue End If c c minimize solution with respect to weighted least squares c norm if operator is singular c If (SINGLR) Call SEPMIN(USOL,IDMN,ZN,ZM,PRTRB) c c return if deferred corrections and a fourth order solution are c not flagged c If (IORD.EQ.2) Return IORD = 2 c c compute new right hand side for fourth order solution c c Call DEFER(COFX,COFY,IDMN,USOL,GRHS) Call DEFER(IDMN,USOL,GRHS) Go To 230 End c Subroutine CHKPRM(INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND,COFX,COFY, c * IDMN,IERROR) Subroutine CHKPRM(INTL,IORDER,A,B,M,MBDCND,C,D,N,NBDCND, * IDMN,IERROR) c c this program checks the input parameters for errors c 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 20 I = 2, M XI = A + FLOAT(I-1) * DLX Call COFX(XI,AI,BI,CI) Do 10 J = 2, N YJ = C + FLOAT(J-1) * DLY Call COFY(YJ,DJ,EJ,FJ) If (AI*DJ.LE.0.0) Then IERROR = 10 Return End If 10 Continue 20 Continue c c no error found c IERROR = 0 Return End c Subroutine CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU,COFX,COFY, c * SINGLR) Subroutine CHKSNG(MBDCND,NBDCND,ALPHA,BETA,GAMA,XNU, * SINGLR) c c this subroutine checks if the pde sepeli c must solve is a singular operator 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.EQ.3) Then If (ALPHA.NE.0.0.OR.BETA.NE.0.0) Return End If If (NBDCND.EQ.3) Then If (GAMA.NE.0.0.OR.XNU.NE.0.0) Return End If c c check that non-derivative coefficient functions c are zero c Do 10 I = IS, MS XI = AIT + FLOAT(I-1) * DLX Call COFX(XI,AI,BI,CI) If (CI.NE.0.0) Return 10 Continue Do 20 J = JS, NS YJ = CIT + FLOAT(J-1) * DLY Call COFY(YJ,DJ,EJ,FJ) If (FJ.NE.0.0) Return 20 Continue c c the operator must be singular if this point is reached c SINGLR = .TRUE. Return End c Subroutine DEFER(COFX,COFY,IDMN,USOL,GRHS) Subroutine DEFER(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 truncation 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) c External COFX, COFY c c compute truncation error approximation over the entire mesh c Do 20 J = JS, NS YJ = CIT + FLOAT(J-1) * DLY Call COFY(YJ,DJ,EJ,FJ) Do 10 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 SEPDX(USOL,IDMN,I,J,UXXX,UXXXX) Call SEPDY(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 truncation if at boundary which is non-periodic c If (.NOT.(KSWX.EQ.1.OR.(I.GT.1.AND.I.LT.K))) Then TX = AI / 3.0 * (UXXXX/4.0+UXXX/DLX) End If If (.NOT.(KSWY.EQ.1.OR.(J.GT.1.AND.J.LT.L))) Then TY = DJ / 3.0 * (UYYYY/4.0+UYYY/DLY) End If GRHS(I,J) = GRHS(I,J) + DLX ** 2 * TX + DLY ** 2 * TY 10 Continue 20 Continue c c reset the right hand side in usol c Do 40 I = IS, MS Do 30 J = JS, NS USOL(I,J) = GRHS(I,J) 30 Continue 40 Continue Return c c revision history--- c c december 1979 first added to nssl c----------------------------------------------------------------------- End Subroutine BLKTRI(IFLG,NP,N,AN,BN,CN,MP,M,AM,BM,CM,IDIMY,Y,IERROR, * W) 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 revision january 1985 c c usage call blktri (iflg,np,n,an,bn,cn,mp,m,am,bm, c cm,idimy,y,ierror,w) c c purpose blktri solves a system of linear equations c of the form c c an(j)*x(i,j-1) + am(i)*x(i-1,j) + c (bn(j)+bm(i))*x(i,j) + cn(j)*x(i,j+1) + c 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 c 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 c discretization of separable elliptic c equations. boundary conditions may be c dirichlet, neumann, or periodic. c c arguments c c on input iflg c c = 0 initialization only. c certain quantities that depend on np, c n, an, bn, and cn are computed and c stored in the work array w. c c = 1 the quantities that were computed c in the initialization are used c to obtain the solution x(i,j). c c note: c a call with iflg=0 takes c approximately one half the time c as a call with iflg = 1. c however, the initialization does c not have to be repeated unless np, c n, an, bn, or cn change. c c np c = 0 if an(1) and cn(n) are not zero, c which corresponds to periodic c bounary conditions. c c = 1 if an(1) and cn(n) are zero. c c n c the number of unknowns in the j-direction. c n must be greater than 4. c the operation count is proportional to c mnlog2(n), hence n should be selected c less than or equal to m. c c an,bn,cn c one-dimensional arrays of length n c that specify the coefficients in the c linear equations given above. c c mp c = 0 if am(1) and cm(m) are not zero, c which corresponds to periodic c boundary conditions. c c = 1 if am(1) = cm(m) = 0 . c c m c the number of unknowns in the i-direction. c m must be greater than 4. c c am,bm,cm c one-dimensional arrays of length m that c specify the coefficients in the linear c equations given above. c c idimy c the row (or first) dimension of the c two-dimensional array y as it appears c in the program calling blktri. c this parameter is used to specify the c variable dimension of y. c idimy must be at least m. c c y c a two-dimensional array that specifies c the values of the right side of the linear c system of equations given above. c y must be dimensioned at least m*n. c c w c a one-dimensional array that must be c provided by the user for work space. c if np=1 define k=int(log2(n))+1 and c set l=2**(k+1) then w must have dimension c (k-2)*l+k+5+max(2n,6m) c c if np=0 define k=int(log2(n-1))+1 and c set l=2**(k+1) then w must have dimension c (k-2)*l+k+5+2n+max(2n,6m) c c **important** c for purposes of checking, the required c dimension of w is computed by blktri and c stored in w(1) in floating point format. c c arguments c c on output y c contains the solution x. c c ierror c an error flag that indicates invalid c input parameters. except for number zer0, c 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 c that depend on the coefficient arrays c an, bn, cn. check these arrays. c = 5 an(j)*cn(j-1) is less than 0 for some j. c c possible reasons for this condition are c 1. the arrays an and cn are not correct c 2. too large a grid spacing was used c in the discretization of the elliptic c equation. c 3. the linear equations resulted from a c partial differential equation which c was not elliptic. c c w c contains intermediate values that must c not be destroyed if blktri will be called c again with iflg=1. w(1) contains the c number of locations required by w in c floating point format. c c c special conditions the algorithm may fail if abs(bm(i)+bn(j)) c is less than abs(am(i))+abs(an(j))+ c abs(cm(i))+abs(cn(j)) c for some i and j. the algorithm will also c fail if an(j)*cn(j-1) is less than zero for c some j. c see the description of the output parameter c ierror. c c i/o none c c precision single c c required library comf and q8qst4, which are automatically loaded c files ncar's cray machines. c c language fortran c c history written by paul swarztrauber at ncar in the c early 1970's. rewritten and released in c january, 1980. c c algorithm generalized cyclic reduction c c portability fortran 66. approximate machine accuracy c is computed in function epmach. c c references swarztrauber,p. and r. sweet, 'efficient c fortran subprograms for the solution of c elliptic equations' c ncar tn/ia-109, july, 1975, 138 pp. c c swarztrauber p. n.,a direct method for c the discrete solution of separable c elliptic equations, s.i.a.m. c j. numer. anal.,11(1974) pp. 1136-1150. c*********************************************************************** Dimension AN(1), BN(1), CN(1), AM(1), BM(1), CM(1), Y(IDIMY,1), * W(1) External PROD, PRODP, CPROD, CPRODP Common /CBLKT/ NPP, K, EPS, CNV, NM, NCMPLX, IK c c the following call is for monitoring library use at ncar c Logical Q8Q4 Save Q8Q4 Data Q8Q4 /.TRUE./ If (Q8Q4) Then c call q8qst4('loclib','blktri','blktri','version 01') Q8Q4 = .FALSE. End If c c test m and n for the proper form c NM = N IERROR = 0 If (M.LT.5) Then IERROR = 1 Else If (NM.LT.3) Then IERROR = 2 Else If (IDIMY.LT.M) Then IERROR = 3 Else NH = N NPP = NP If (NPP.NE.0) Then NH = NH + 1 End If IK = 2 K = 1 10 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 If (NPP.NE.0) Then c c divide w into working sub arrays c IW1 = IWAH IWBH = IW1 + NM W(1) = FLOAT(IW1-1+MAX0(2*NM,6*M)) Else IWBH = IWAH + NM + NM IW1 = IWBH W(1) = FLOAT(IW1-1+MAX0(2*NM,6*M)) NM = NM - 1 End If c c subroutine comp b computes the roots of the b polynomials c If (IERROR.EQ.0) Then 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)) Else If (MP.NE.0) Then c c subroutine blktr1 solves the linear system c 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 End If End If End If End If End If Return End Subroutine BLKTR1(N,AN,BN,CN,M,AM,BM,CM,IDIMY,Y,B,W1,W2,W3,WD,WW, * WU,PRDCT,CPRDCT) c c blktr1 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 prodp or prod depending on whether the boundary c conditions in the m direction are periodic or not c cprdct is either cprodp or cprod which are the complex versions c of prodp and prod. these are called in the event that some c of the roots of the b sub p polynomial 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) Common /CBLKT/ NPP, K, EPS, CNV, NM, NCMPLX, IK c c begin reduction phase c KDO = K - 1 Do 40 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 30 I = I4, IF, I4 If (I.LE.NM) Then IPI1 = I + I1 IPI2 = I + I2 IPI3 = I + I3 Call INDXC(I,IR,IDXC,NC) If (I.LT.IF) Then 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.GT.NM) Then Do 10 J = 1, M W3(J) = 0. W2(J) = 0. 10 Continue Else 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) End If Do 20 J = 1, M Y(J,I) = W1(J) + W2(J) + Y(J,I) 20 Continue End If End If 30 Continue 40 Continue If (NPP.EQ.0) Then c c the periodic case is treated using the capacitance matrix method c 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 50 J = 1, M W2(J) = W1(J) 50 Continue Do 70 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 60 J = 1, M W1(J) = Y(J,I) + W1(J) 60 Continue Call PRDCT(NZ,B(IZ),NM1,B(IM1),NP1,B(IP1),0,DUM,W1,W1,M,AM,BM, * CM,WD,WW,WU) 70 Continue Do 100 LL = 2, K L = K - LL + 1 IR = L - 1 I2 = 2 ** IR I1 = I2 / 2 I4 = I2 + I2 IFD = IF - I2 Do 90 I = I2, IFD, I4 If (I-I2-IZR.EQ.0) Then If (I.GT.NM) Go To 100 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 80 J = 1, M W2(J) = Y(J,I) + W2(J) 80 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.EQ.NM) Go To 110 End If 90 Continue 100 Continue 110 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 INDXB(IF/2,K-1,IM1,NM1) Call INDXB(IF,K-1,IP,NP) If (NCMPLX.NE.0) Then 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) Else 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) End If Do 130 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) 130 Continue Do 150 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 140 J = 1, M Y(J,I) = Y(J,I) - W1(J) 140 Continue 150 Continue c IZR = NM Do 180 L = 1, KDO IR = L - 1 I2 = 2 ** IR I1 = I2 / 2 I3 = I2 + I1 I4 = I2 + I2 IRM1 = IR - 1 Do 170 I = I4, IF, I4 IPI1 = I + I1 IPI2 = I + I2 IPI3 = I + I3 If (IPI2.NE.IZR) Then If (I-IZR) 170, 180, 170 End If 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 160 J = 1, M Y(J,I) = Y(J,I) - W2(J) 160 Continue IZR = I Go To 180 170 Continue 180 Continue End If c c begin back substitution phase c Do 230 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 220 I = I2, IFD, I4 If (I.LE.NM) Then 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.LE.I2) Then Do 190 J = 1, M W1(J) = 0. 190 Continue Else Call PRDCT(NM1,B(IM1),0,DUM,0,DUM,NA,AN(IDXA),Y(1,IMI2), * W1,M,AM,BM,CM,WD,WW,WU) End If If (IPI2.GT.NM) Then Do 200 J = 1, M W2(J) = 0. 200 Continue Else Call PRDCT(NP1,B(IP1),0,DUM,0,DUM,NC,CN(IDXC),Y(1,IPI2), * W2,M,AM,BM,CM,WD,WW,WU) End If Do 210 J = 1, M W1(J) = Y(J,I) + W1(J) + W2(J) 210 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) End If 220 Continue 230 Continue Return End Subroutine INDXB(I,IR,IDX,IDP) c c b(idx) is the location of the first root of the b(i,ir) polynomial c Common /CBLKT/ NPP, K, EPS, CNV, NM, NCMPLX, IK IDP = 0 If (IR) 30, 10, 20 10 If (I.GT.NM) Go To 30 IDX = I IDP = 1 Return 20 IZH = 2 ** IR ID = I - IZH - IZH IDX = ID + ID + (IR-1) * IK + IR + (IK-I) / IZH + 4 IPL = IZH - 1 IDP = IZH + IZH - 1 If (I-IPL-NM.GT.0) Then IDP = 0 Return End If If (I+IPL-NM.GT.0) Then IDP = NM + IPL - I + 1 End If 30 Return End Subroutine INDXA(I,IR,IDXA,NA) Common /CBLKT/ NPP, K, EPS, CNV, NM, NCMPLX, IK NA = 2 ** IR IDXA = I - NA + 1 If (I.GT.NM) Then NA = 0 End If Return End Subroutine INDXC(I,IR,IDXC,NC) Common /CBLKT/ NPP, K, EPS, CNV, NM, NCMPLX, IK NC = 2 ** IR IDXC = I If (IDXC+NC-1-NM.GT.0) Then NC = 0 End If Return End Subroutine PROD(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,W,U) c c prod applies a sequence of matrix operations to the vector x and c stores the result in y 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,w,u are working arrays c is determines whether or not a change in sign is made c Dimension A(1), B(1), C(1), X(1), Y(1), D(1), W(1), BD(1), * BM1(1), BM2(1), AA(1), U(1) Do 10 J = 1, M W(J) = X(J) Y(J) = W(J) 10 Continue MM = M - 1 ID = ND IBR = 0 M1 = NM1 M2 = NM2 IA = NA 20 If (IA.GT.0) Then RT = AA(IA) If (ND.EQ.0) RT = -RT IA = IA - 1 c c scalar multiplication c Do 30 J = 1, M Y(J) = RT * W(J) 30 Continue End If If (ID.GT.0) Then RT = BD(ID) ID = ID - 1 If (ID.EQ.0) IBR = 1 c c begin solution to system c D(M) = A(M) / (B(M)-RT) W(M) = Y(M) / (B(M)-RT) Do 40 J = 2, MM K = M - J DEN = B(K+1) - RT - 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 40 Continue DEN = B(1) - RT - C(1) * D(2) W(1) = 1. If (DEN.NE.0) Then W(1) = (Y(1)-C(1)*W(2)) / DEN End If Do 50 J = 2, M W(J) = W(J) - D(J) * W(J-1) 50 Continue If (NA) 80, 80, 20 60 Do 70 J = 1, M Y(J) = W(J) 70 Continue IBR = 1 Go To 20 80 If (M1.LE.0) Then If (M2) 60, 60, 90 End If If (M2.GT.0) Then If (ABS(BM1(M1)).LE.ABS(BM2(M2))) Go To 90 End If If (IBR.LE.0) Then If (ABS(BM1(M1)-BD(ID)).LT.ABS(BM1(M1)-RT)) Go To 60 End If RT = RT - BM1(M1) M1 = M1 - 1 Go To 100 90 If (IBR.LE.0) Then If (ABS(BM2(M2)-BD(ID)).LT.ABS(BM2(M2)-RT)) Go To 60 End If RT = RT - BM2(M2) M2 = M2 - 1 100 Do 110 J = 1, M Y(J) = Y(J) + RT * W(J) 110 Continue Go To 20 End If Return End Subroutine PRODP(ND,BD,NM1,BM1,NM2,BM2,NA,AA,X,Y,M,A,B,C,D,U,W) c c prodp applies a sequence of matrix operations to the vector x and c stores the result in y periodic boundary conditions 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,w are working arrays c is determines whether or not a change in sign is made c Dimension A(1), B(1), C(1), X(1), Y(1), D(1), U(1), BD(1), * BM1(1), BM2(1), AA(1), W(1) Do 10 J = 1, M Y(J) = X(J) W(J) = Y(J) 10 Continue MM = M - 1 MM2 = M - 2 ID = ND IBR = 0 M1 = NM1 M2 = NM2 IA = NA 20 If (IA.GT.0) Then RT = AA(IA) If (ND.EQ.0) RT = -RT IA = IA - 1 Do 30 J = 1, M Y(J) = RT * W(J) 30 Continue End If If (ID.GT.0) Then RT = BD(ID) ID = ID - 1 If (ID.EQ.0) IBR = 1 c c begin solution to system c BH = B(M) - RT YM = Y(M) DEN = B(1) - RT D(1) = C(1) / DEN U(1) = A(1) / DEN W(1) = Y(1) / DEN V = C(M) If (MM2.GE.2) Then Do 40 J = 2, MM2 DEN = B(J) - RT - A(J) * D(J-1) D(J) = C(J) / DEN U(J) = -A(J) * U(J-1) / DEN W(J) = (Y(J)-A(J)*W(J-1)) / DEN BH = BH - V * U(J-1) YM = YM - V * W(J-1) V = -V * D(J-1) 40 Continue End If DEN = B(M-1) - RT - A(M-1) * D(M-2) D(M-1) = (C(M-1)-A(M-1)*U(M-2)) / DEN W(M-1) = (Y(M-1)-A(M-1)*W(M-2)) / DEN AM = A(M) - V * D(M-2) BH = BH - V * U(M-2) YM = YM - V * W(M-2) DEN = BH - AM * D(M-1) If (DEN.NE.0) Then W(M) = (YM-AM*W(M-1)) / DEN Else W(M) = 1. End If W(M-1) = W(M-1) - D(M-1) * W(M) Do 50 J = 2, MM K = M - J W(K) = W(K) - D(K) * W(K+1) - U(K) * W(M) 50 Continue If (NA) 80, 80, 20 60 Do 70 J = 1, M Y(J) = W(J) 70 Continue IBR = 1 Go To 20 80 If (M1.LE.0) Then If (M2) 60, 60, 90 End If If (M2.GT.0) Then If (ABS(BM1(M1)).LE.ABS(BM2(M2))) Go To 90 End If If (IBR.LE.0) Then If (ABS(BM1(M1)-BD(ID)).LT.ABS(BM1(M1)-RT)) Go To 60 End If RT = RT - BM1(M1) M1 = M1 - 1 Go To 100 90 If (IBR.LE.0) Then If (ABS(BM2(M2)-BD(ID)).LT.ABS(BM2(M2)-RT)) Go To 60 End If RT = RT - BM2(M2) M2 = M2 - 1 100 Do 110 J = 1, M Y(J) = Y(J) + RT * W(J) 110 Continue Go To 20 End If 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 Complex Y, D, W, BD, CRT, DEN, Y1, Y2 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 10 J = 1, M Y(J) = CMPLX(X(J),0.) 10 Continue MM = M - 1 ID = ND M1 = NM1 M2 = NM2 IA = NA 20 IFLG = 0 If (ID.GT.0) Then 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 30 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 30 Continue DEN = B(1) - CRT - C(1) * D(2) If (CABS(DEN).NE.0) Then Y(1) = (Y(1)-C(1)*W(2)) / DEN Else Y(1) = (1.,0.) End If Do 40 J = 2, M Y(J) = W(J) - D(J) * Y(J-1) 40 Continue End If If (M1.LE.0) Then If (M2.LE.0) Go To 60 RT = BM2(M2) M2 = M2 - 1 Else If (M2.LE.0) Then RT = BM1(M1) M1 = M1 - 1 Else If (ABS(BM1(M1)).GT.ABS(BM2(M2))) Then RT = BM1(M1) M1 = M1 - 1 Else RT = BM2(M2) M2 = M2 - 1 End If End If End If Y1 = (B(1)-RT) * Y(1) + C(1) * Y(2) If (MM.GE.2) Then c c matrix multiplication c Do 50 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 50 Continue End If Y(M) = A(M) * Y(M-1) + (B(M)-RT) * Y(M) Y(M-1) = Y1 IFLG = 1 Go To 20 60 If (IA.GT.0) Then RT = AA(IA) IA = IA - 1 IFLG = 1 c c scalar multiplication c Do 70 J = 1, M Y(J) = RT * Y(J) 70 Continue End If If (IFLG.GT.0) Go To 20 Do 80 J = 1, M YY(J) = REAL(Y(J)) 80 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 10 J = 1, M Y(J) = CMPLX(X(J),0.) 10 Continue MM = M - 1 MM2 = M - 2 ID = ND M1 = NM1 M2 = NM2 IA = NA 20 IFLG = 0 If (ID.GT.0) Then 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.GE.2) Then Do 30 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) 30 Continue End If 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).NE.0) Then Y(M) = (YM-AM*Y(M-1)) / DEN Else Y(M) = (1.,0.) End If Y(M-1) = Y(M-1) - D(M-1) * Y(M) Do 40 J = 2, MM K = M - J Y(K) = Y(K) - D(K) * Y(K+1) - U(K) * Y(M) 40 Continue End If If (M1.LE.0) Then If (M2.LE.0) Go To 60 RT = BM2(M2) M2 = M2 - 1 Else If (M2.LE.0) Then RT = BM1(M1) M1 = M1 - 1 Else If (ABS(BM1(M1)).GT.ABS(BM2(M2))) Then RT = BM1(M1) M1 = M1 - 1 Else RT = BM2(M2) M2 = M2 - 1 End If End If End If c c matrix multiplication c YH = Y(1) Y1 = (B(1)-RT) * Y(1) + C(1) * Y(2) + A(1) * Y(M) If (MM.GE.2) Then Do 50 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 50 Continue End If Y(M) = A(M) * Y(M-1) + (B(M)-RT) * Y(M) + C(M) * YH Y(M-1) = Y1 IFLG = 1 Go To 20 60 If (IA.GT.0) Then RT = AA(IA) IA = IA - 1 IFLG = 1 c c scalar multiplication c Do 70 J = 1, M Y(J) = RT * Y(J) 70 Continue End If If (IFLG.GT.0) Go To 20 Do 80 J = 1, M YY(J) = REAL(Y(J)) 80 Continue Return End Subroutine PPADD(N,IERROR,A,C,CBP,BP,BH) c c ppadd 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 /CBLKT/ NPP, K, EPS, CNV, NM, NCMPLX, IK External PSGF, PPSPF, PPSGF SCNV = SQRT(CNV) IZ = N IZM = IZ - 1 IZM2 = IZ - 2 If (BP(N)-BP(1)) 10, 240, 30 10 Do 20 J = 1, N NT = N - J BH(J) = BP(NT+1) 20 Continue Go To 50 30 Do 40 J = 1, N BH(J) = BP(J) 40 Continue 50 NCMPLX = 0 MODIZ = MOD(IZ,2) IS = 1 If (MODIZ.NE.0) Then If (A(1)) 80, 240, 60 End If 60 XL = BH(1) DB = BH(3) - BH(1) 70 XL = XL - DB If (PSGF(XL,IZ,C,A,BH).LE.0) Go To 70 SGN = -1. CBP(1) = CMPLX(BSRH(XL,BH(1),IZ,C,A,BH,PSGF,SGN),0.) IS = 2 80 IF = IZ - 1 If (MODIZ.NE.0) Then If (A(1)) 90, 240, 110 End If 90 XR = BH(IZ) DB = BH(IZ) - BH(IZ-2) 100 XR = XR + DB If (PSGF(XR,IZ,C,A,BH).LT.0) Go To 100 SGN = 1. CBP(IZ) = CMPLX(BSRH(BH(IZ),XR,IZ,C,A,BH,PSGF,SGN),0.) IF = IZ - 2 110 Do 180 IG = IS, IF, 2 XL = BH(IG) XR = BH(IG+1) SGN = -1. XM = BSRH(XL,XR,IZ,C,A,BH,PPSPF,SGN) PSG = PSGF(XM,IZ,C,A,BH) If (ABS(PSG).GT.EPS) Then If (PSG*PPSGF(XM,IZ,C,A,BH)) 120, 130, 140 c c case of a real zero c 120 SGN = 1. CBP(IG) = CMPLX(BSRH(BH(IG),XM,IZ,C,A,BH,PSGF,SGN),0.) SGN = -1. CBP(IG+1) = CMPLX(BSRH(XM,BH(IG+1),IZ,C,A,BH,PSGF,SGN),0.) Go To 180 End If c c case of a multiple zero c 130 CBP(IG) = CMPLX(XM,0.) CBP(IG+1) = CMPLX(XM,0.) Go To 180 c c case of a complex zero c 140 IT = 0 ICV = 0 CX = CMPLX(XM,0.) 150 FSG = (1.,0.) HSG = (1.,0.) FP = (0.,0.) FPP = (0.,0.) Do 160 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 160 Continue If (MODIZ.EQ.0) Then F = (1.,0.) - FSG - HSG Else F = (1.,0.) + FSG + HSG End If I3 = 0 If (CABS(FP).GT.0) Then I3 = 1 R3 = -F / FP End If I2 = 0 If (CABS(FPP).GT.0) Then I2 = 1 CDIS = CSQRT(FP**2-2.*F*FPP) R1 = CDIS - FP R2 = -FP - CDIS If (CABS(R1).GT.CABS(R2)) Then R1 = R1 / FPP Else R1 = R2 / FPP End If R2 = 2. * F / FPP / R1 If (CABS(R2).LT.CABS(R1)) R1 = R2 If (I3.LE.0) Go To 170 If (CABS(R3).LT.CABS(R1)) R1 = R3 Else R1 = R3 End If 170 CX = CX + R1 IT = IT + 1 If (IT.GT.50) Go To 240 If (CABS(R1).GT.SCNV) Go To 150 If (ICV.LE.0) Then ICV = 1 Go To 150 End If CBP(IG) = CX CBP(IG+1) = CONJG(CX) 180 Continue If (CABS(CBP(N))-CABS(CBP(1))) 190, 240, 210 190 NHALF = N / 2 Do 200 J = 1, NHALF NT = N - J CX = CBP(J) CBP(J) = CBP(NT+1) CBP(NT+1) = CX 200 Continue 210 NCMPLX = 1 Do 220 J = 2, IZ If (AIMAG(CBP(J)).NE.0) Go To 250 220 Continue NCMPLX = 0 Do 230 J = 2, IZ BP(J) = REAL(CBP(J)) 230 Continue Go To 250 240 IERROR = 4 250 Continue Return End Function PSGF(X,IZ,C,A,BH) Dimension A(1), C(1), BH(1) FSG = 1. HSG = 1. Do 10 J = 1, IZ DD = 1. / (X-BH(J)) FSG = FSG * A(J) * DD HSG = HSG * C(J) * DD 10 Continue If (MOD(IZ,2).EQ.0) Then PSGF = 1. - FSG - HSG Return End If PSGF = 1. + FSG + HSG Return End Function BSRH(XLL,XRR,IZ,C,A,BH,F,SGN) Dimension A(1), C(1), BH(1) Common /CBLKT/ NPP, K, EPS, CNV, NM, NCMPLX, IK XL = XLL XR = XRR DX = .5 * ABS(XR-XL) 10 X = .5 * (XL+XR) If (SGN*F(X,IZ,C,A,BH)) 30, 50, 20 20 XR = X Go To 40 30 XL = X 40 DX = .5 * DX If (DX.GT.CNV) Go To 10 50 BSRH = .5 * (XL+XR) Return End Function PPSGF(X,IZ,C,A,BH) Dimension A(1), C(1), BH(1) SUM = 0. Do 10 J = 1, IZ SUM = SUM - 1. / (X-BH(J)) ** 2 10 Continue PPSGF = SUM Return End Function PPSPF(X,IZ,C,A,BH) Dimension A(1), C(1), BH(1) SUM = 0. Do 10 J = 1, IZ SUM = SUM + 1. / (X-BH(J)) 10 Continue PPSPF = SUM 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(DUM) BNORM = ABS(BN(1)) Do 10 J = 2, NM BNORM = AMAX1(BNORM,ABS(BN(J))) ARG = AN(J) * CN(J-1) If (ARG.LT.0) Go To 110 B(J) = SIGN(SQRT(ARG),AN(J)) 10 Continue CNV = EPS * BNORM IF = 2 ** K KDO = K - 1 Do 50 L = 1, KDO IR = L - 1 I2 = 2 ** IR I4 = I2 + I2 IPL = I4 - 1 IFD = IF - I4 Do 40 I = I4, IFD, I4 Call INDXB(I,L,IB,NB) If (NB.LE.0) Go To 50 JS = I - IPL JF = JS + NB - 1 LS = 0 Do 20 J = JS, JF LS = LS + 1 BH(LS) = BN(J) AH(LS) = B(J) 20 Continue Call TEVLS(NB,BH,AH,IERROR) If (IERROR.NE.0) Go To 100 LH = IB - 1 Do 30 J = 1, NB LH = LH + 1 B(LH) = -BH(J) 30 Continue 40 Continue 50 Continue Do 60 J = 1, NM B(J) = -BN(J) 60 Continue If (NPP.EQ.0) Then NMP = NM + 1 NB = NM + NMP Do 70 J = 1, NB L1 = MOD(J-1,NMP) + 1 L2 = MOD(J+NM-1,NMP) + 1 ARG = AN(L1) * CN(L2) If (ARG.LT.0) Go To 110 BH(J) = SIGN(SQRT(ARG),-AN(L1)) AH(J) = -BN(L1) 70 Continue Call TEVLS(NB,AH,BH,IERROR) If (IERROR.NE.0) Go To 100 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 80 D1 = ABS(B(J1)-B(J2-1)) D2 = ABS(B(J1)-B(J2)) D3 = ABS(B(J1)-B(J2+1)) If (.NOT.((D2.LT.D1).AND.(D2.LT.D3))) Then B(LH) = B(J2) J2 = J2 + 1 LH = LH + 1 If (J2-N2M2) 80, 80, 90 End If J2 = J2 + 1 J1 = J1 + 1 If (J2.LE.N2M2) Go To 80 90 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)) End If Return 100 IERROR = 4 Return 110 IERROR = 5 Return End Subroutine TEVLS(N,D,E2,IERR) c Integer I, J, L, M, N, II, L1, MML, IERR Real D(N), E2(N) Real B, C, F, G, H, P, R, S, MACHEP c c real sqrt,abs,sign c Common /CBLKT/ NPP, K, MACHEP, CNV, NM, NCMPLX, IK c c this subroutine is a modification of the eispack subroutine tqlrat c algorithm 464, comm. acm 16, 689(1973) by reinsch. c c this subroutine finds the eigenvalues of a symmetric c tridiagonal matrix by the rational ql method. c c on input- c c n is the order of the matrix, c c d contains the diagonal elements of the input matrix, c c e2 contains the subdiagonal elements of the c input matrix in its last n-1 positions. e2(1) is arbitrary. c c on output- c c d contains the eigenvalues in ascending order. if an c error exit is made, the eigenvalues are correct and c ordered for indices 1,2,...ierr-1, but may not be c the smallest eigenvalues, c c e2 has been destroyed, c c ierr is set to c zero for normal return, c j if the j-th eigenvalue has not been c determined after 30 iterations. c c questions and comments should be directed to b. s. garbow, c applied mathematics division, argonne national laboratory c c c ********** machep is a machine dependent parameter specifying c the relative precision of floating point arithmetic. c c ********** c IERR = 0 If (N.NE.1) Then c Do 10 I = 2, N E2(I-1) = E2(I) * E2(I) 10 Continue c F = 0.0 B = 0.0 E2(N) = 0.0 c Do 90 L = 1, N J = 0 H = MACHEP * (ABS(D(L))+SQRT(E2(L))) If (B.LE.H) Then B = H C = B * B End If c c ********** look for small squared sub-diagonal element ********** c Do 20 M = L, N If (E2(M).LE.C) Go To 30 c c ********** e2(n) is always zero, so there is no exit c through the bottom of the loop ********** c 20 Continue c 30 If (M.NE.L) Then 40 If (J.EQ.30) Go To 110 J = J + 1 c c ********** form shift ********** c L1 = L + 1 S = SQRT(E2(L)) G = D(L) P = (D(L1)-G) / (2.0*S) R = SQRT(P*P+1.0) D(L) = S / (P+SIGN(R,P)) H = G - D(L) c Do 50 I = L1, N D(I) = D(I) - H 50 Continue c F = F + H c c ********** rational ql transformation ********** c G = D(M) If (G.EQ.0.0) G = B H = G S = 0.0 MML = M - L c c ********** for i=m-1 step -1 until l do -- ********** c Do 60 II = 1, MML I = M - II P = G * H R = P + E2(I) E2(I+1) = S * R S = E2(I) / R D(I+1) = H + S * (H+D(I)) G = D(I) - E2(I) / G If (G.EQ.0.0) G = B H = G * P / R 60 Continue c E2(L) = S * G D(L) = H c c ********** guard against underflowed h ********** c If (H.NE.0.0) Then If (ABS(E2(L)).GT.ABS(C/H)) Then E2(L) = H * E2(L) If (E2(L).NE.0.0) Go To 40 End If End If End If P = D(L) + F c c ********** order eigenvalues ********** c If (L.NE.1) Then c c ********** for i=l step -1 until 2 do -- ********** c Do 70 II = 2, L I = L + 2 - II If (P.GE.D(I-1)) Go To 80 D(I) = D(I-1) 70 Continue End If c I = 1 80 D(I) = P 90 Continue c If (ABS(D(N)).GE.ABS(D(1))) Go To 120 NHALF = N / 2 Do 100 I = 1, NHALF NTOP = N - I DHOLD = D(I) D(I) = D(NTOP+1) D(NTOP+1) = DHOLD 100 Continue Go To 120 c c ********** set error -- no convergence to an c eigenvalue after 30 iterations ********** c 110 IERR = L End If 120 Return c c ********** last card of tqlrat ********** c c c revision history--- c c december 1979 first added to nssl c----------------------------------------------------------------------- End c package comf the entries in this package are lowlevel c entries, supporting ulib packages blktri c and cblktri. that is, these routines are c not called directly by users, but rather c by entries within blktri and cblktri. c description of entries epmach and pimach c follow below. c c latest revision january 1985 c c special conditions none c c i/o none c c precision single c c required library none c files c c language fortran c ******************************************************************** c c function epmach (dum) c c purpose to compute an approximate machine accuracy c epsilon according to the following definition: c epsilon is the smallest number such that c (1.+epsilon).gt.1.) c c usage eps = epmach (dum) c c arguments c on input dum c dummy value c c arguments c on output none c c history the original version, written when the c blktri package was converted from the c cdc 7600 to run on the cray-1, calculated c machine accuracy by successive divisions c by 10. use of this constant caused blktri c to compute solutions on the cray-1 with four c fewer places of accuracy than the version c on the 7600. it was found that computing c machine accuracy by successive divisions c of 2 produced a machine accuracy 29% less c than the value generated by successive c divisions by 10, and that use of this c machine constant in the blktri package c recovered the accuracy that appeared to c be lost on conversion. c c algorithm computes machine accuracy by successive c divisions of two. c c portability this code will execute on machines other c than the cray1, but the returned value may c be unsatisfactory. see history above. c ******************************************************************** c c function pimach (dum) c c purpose to supply the value of the constant pi c correct to machine precision where c pi=3.141592653589793238462643383279502884197 c 1693993751058209749446 c c usage pi = pimach (dum) c c arguments c on input dum c dummy value c c arguments c on output none c c algorithm the value of pi is set in a constant. c c portability this entry is portable, but users should c check to see whether greater accuracy is c required. c c*********************************************************************** Function EPMACH(DUM) Common /VALUE/ V EPS = 1. 10 EPS = EPS / 2. Call STORE(EPS+1.) If (V.GT.1.) Go To 10 EPMACH = 100. * EPS Return End Subroutine STORE(X) Common /VALUE/ V V = X Return End Function PIMACH(DUM) c pi=3.1415926535897932384626433832795028841971693993751058209749446 c PIMACH = 3.14159265358979 Return End c package sepaux contains no user entry points. c c latest revision march 1985 c c purpose this package contains auxiliary routines for c ncar public software packages such as sepeli c and sepx4. c c usage since this package contains no user entries, c no usage instructions or argument descriptions c are given here. c c special conditions none c c i/o none c c precision single c c required library none c files c c language fortran c c history developed in the late 1970's by john c. adams c of ncar's scienttific computing division. c c portability fortran 66 c ********************************************************************** Subroutine SEPORT(USOL,IDMN,ZN,ZM,PERTRB) c c this subroutine orthoganalizes the array usol with respect to c the constant array in a weighted least squares norm c Common /SPLP/ KSWX, KSWY, K, L, AIT, BIT, CIT, DIT, MIT, NIT, IS, * MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4 Dimension USOL(IDMN,1), ZN(1), ZM(1) ISTR = IS IFNL = MS JSTR = JS JFNL = NS c c compute weighted inner products c UTE = 0.0 ETE = 0.0 Do 20 I = IS, MS II = I - IS + 1 Do 10 J = JS, NS JJ = J - JS + 1 ETE = ETE + ZM(II) * ZN(JJ) UTE = UTE + USOL(I,J) * ZM(II) * ZN(JJ) 10 Continue 20 Continue c c set perturbation parameter c PERTRB = UTE / ETE c c subtract off constant pertrb c Do 40 I = ISTR, IFNL Do 30 J = JSTR, JFNL USOL(I,J) = USOL(I,J) - PERTRB 30 Continue 40 Continue Return End Subroutine SEPMIN(USOL,IDMN,ZN,ZM,PERTB) c c this subroutine orhtogonalizes the array usol with respect to c the constant array in a weighted least squares norm c Common /SPLP/ KSWX, KSWY, K, L, AIT, BIT, CIT, DIT, MIT, NIT, IS, * MS, JS, NS, DLX, DLY, TDLX3, TDLY3, DLX4, DLY4 Dimension USOL(IDMN,1), ZN(1), ZM(1) c c entry at sepmin occurrs when the final solution is c to be minimized with respect to the weighted c least squares norm c ISTR = 1 IFNL = K JSTR = 1 JFNL = L c c compute weighted inner products c UTE = 0.0 ETE = 0.0 Do 20 I = IS, MS II = I - IS + 1 Do 10 J = JS, NS JJ = J - JS + 1 ETE = ETE + ZM(II) * ZN(JJ) UTE = UTE + USOL(I,J) * ZM(II) * ZN(JJ) 10 Continue 20 Continue c c set perturbation parameter c PERTRB = UTE / ETE c c subtract off constant pertrb c Do 40 I = ISTR, IFNL Do 30 J = JSTR, JFNL USOL(I,J) = USOL(I,J) - PERTRB 30 Continue 40 Continue Return End Subroutine SEPTRI(N,A,B,C,D,U,Z) c c this subroutine solves for a non-zero eigenvector corresponding c to the zero eigenvalue of the transpose of the rank c deficient one matrix with subdiagonal a, diagonal b, and c superdiagonal c , with a(1) in the (1,n) position, with c c(n) in the (n,1) position, and all other elements zero. c Dimension A(N), B(N), C(N), D(N), U(N), Z(N) BN = B(N) D(1) = A(2) / B(1) V = A(1) U(1) = C(N) / B(1) NM2 = N - 2 Do 10 J = 2, NM2 DEN = B(J) - C(J-1) * D(J-1) D(J) = A(J+1) / DEN U(J) = -C(J-1) * U(J-1) / DEN BN = BN - V * U(J-1) V = -V * D(J-1) 10 Continue DEN = B(N-1) - C(N-2) * D(N-2) D(N-1) = (A(N)-C(N-2)*U(N-2)) / DEN AN = C(N-1) - V * D(N-2) BN = BN - V * U(N-2) DEN = BN - AN * D(N-1) c c set last component equal to one c Z(N) = 1.0 Z(N-1) = -D(N-1) NM1 = N - 1 Do 20 J = 2, NM1 K = N - J Z(K) = -D(K) * Z(K+1) - U(K) * Z(N) 20 Continue Return End Subroutine SEPDX(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.LE.2.OR.I.GE.(K-1)) Then If (I.NE.1) Then If (I.EQ.2) Go To 10 If (I.EQ.K-1) Go To 20 If (I.EQ.K) Go To 30 End If c c compute partial derivative approximations at x=a c If (KSWX.NE.1) Then 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 End If c c periodic at x=a c 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 10 If (KSWX.NE.1) Then 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 End If c c periodic at x=a+dlx c 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 End If c c compute partial derivative approximations on the interior c 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 20 If (KSWX.NE.1) Then 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 End If c c periodic at x=b-dlx c 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 30 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 SEPDY(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.LE.2.OR.J.GE.(L-1)) Then If (J.NE.1) Then If (J.EQ.2) Go To 10 If (J.EQ.L-1) Go To 20 If (J.EQ.L) Go To 30 End If c c compute partial derivative approximations at y=c c If (KSWY.NE.1) Then 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 End If c c periodic at x=a c 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 10 If (KSWY.NE.1) Then 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 End If c c periodic at y=c+dly c 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 End If c c compute partial derivative approximations on the interior c 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 20 If (KSWY.NE.1) Then 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 End If c c periodic at y=d-dly c 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 30 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 c c revision history--- c c december 1979 first added to nssl c----------------------------------------------------------------------- End Subroutine COFX(XX,AFUN,BFUN,CFUN) c c subroutine to compute the coefficients of the elliptic equation c solved to fill in the grid. it is passed (along with cofy) to c the elliptic solver sepeli. Include 'cgridgen.inc' I = NINT(XX) DXDI = 0.5 * (SXI(I+1)-SXI(I-1)) AFUN = 1. / DXDI ** 2 BFUN = (-SXI(I+1)+2.*SXI(I)-SXI(I-1)) / DXDI ** 3 CFUN = 0. Return End Subroutine COFY(YY,DFUN,EFUN,FFUN) Include 'cgridgen.inc' J = NINT(YY) DEDJ = 0.5 * (SETA(J+1)-SETA(J-1)) DFUN = 1. / DEDJ ** 2 EFUN = (-SETA(J+1)+2.*SETA(J)-SETA(J-1)) / DEDJ ** 3 FFUN = 0. Return End