Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w]=hstcrt(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w);
```function [a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w]=hstcrt(a,b,m,mbdcnd,bda,bdb,c,d,n,nbdcnd,bdc,bdd,elmbda,f,idimf,pertrb,ierror,w);
persistent deltax deltay delxsq delysq i id2 id3 id4 ierr1 j mp mperod np nperod s st2 twdelx twdely twdysq ;

if isempty(deltax), deltax=0; end;
if isempty(deltay), deltay=0; end;
if isempty(delxsq), delxsq=0; end;
if isempty(delysq), delysq=0; end;
if isempty(s), s=0; end;
if isempty(st2), st2=0; end;
if isempty(twdelx), twdelx=0; end;
if isempty(twdely), twdely=0; end;
if isempty(twdysq), twdysq=0; end;
%***BEGIN PROLOGUE  HSTCRT
%***PURPOSE  Solve the standard five-point finite difference
%            approximation on a staggered grid to the Helmholtz equation
%            in Cartesian coordinates.
%***LIBRARY   SLATEC (FISHPACK)
%***CATEGORY  I2B1A1A
%***TYPE      SINGLE PRECISION (HSTCRT-S)
%***KEYWORDS  ELLIPTIC, FISHPACK, HELMHOLTZ, PDE
%           Swarztrauber, P. N., (NCAR)
%           Sweet, R., (NCAR)
%***DESCRIPTION
%
%      HSTCRT solves the standard five-point finite difference
%      approximation on a staggered grid to the Helmholtz equation in
%      Cartesian coordinates
%
%      (d/dX,dU/dX) + (d/dY,dU/dY) + LAMBDA*U = F(X,Y)
%
%     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
%     * * * * * * * *    Parameter Description     * * * * * * * * * *
%
%             * * * * * *   On Input    * * * * * *
%
%    A,B
%      The range of X, i.e. A .LE. X .LE. B.  A must be less than B.
%
%    M
%      The number of grid points in the interval (A,B).  The grid points
%      in the X-direction are given by X(I) = A + (I-0.5)dX for
%      I=1,2,...,M where dX =(B-A)/M.  M must be greater than 2.
%
%    MBDCND
%      Indicates the type of boundary conditions at X = A and X = B.
%
%      = 0  If the solution is periodic in X,
%           U(M+I,J) = U(I,J).
%
%      = 1  If the solution is specified at X = A and X = B.
%
%      = 2  If the solution is specified at X = A and the derivative
%           of the solution with respect to X is specified at X = B.
%
%      = 3  If the derivative of the solution with respect to X is
%           specified at X = A  and X = B.
%
%      = 4  If the derivative of the solution with respect to X is
%           specified at X = A  and the solution is specified at X = B.
%
%    BDA
%      A one-dimensional array of length N that specifies the boundary
%      values (if any) of the solution at X = A.  When MBDCND = 1 or 2,
%
%               BDA(J) = U(A,Y(J)) ,          J=1,2,...,N.
%
%      When MBDCND = 3 or 4,
%
%               BDA(J) = (d/dX)U(A,Y(J)) ,    J=1,2,...,N.
%
%    BDB
%      A one-dimensional array of length N that specifies the boundary
%      values of the solution at X = B.  When MBDCND = 1 or 4
%
%               BDB(J) = U(B,Y(J)) ,          J=1,2,...,N.
%
%      When MBDCND = 2 or 3
%
%               BDB(J) = (d/dX)U(B,Y(J)) ,    J=1,2,...,N.
%
%    C,D
%      The range of Y, i.e. C .LE. Y .LE. D.  C must be less
%      than D.
%
%    N
%      The number of unknowns in the interval (C,D).  The unknowns in
%      the Y-direction are given by Y(J) = C + (J-0.5)DY,
%      J=1,2,...,N, where DY = (D-C)/N.  N must be greater than 2.
%
%    NBDCND
%      Indicates the type of boundary conditions at Y = C
%      and Y = D.
%
%      = 0  If the solution is periodic in Y, i.e.
%           U(I,J) = U(I,N+J).
%
%      = 1  If the solution is specified at Y = C and Y = D.
%
%      = 2  If the solution is specified at Y = C and the derivative
%           of the solution with respect to Y is specified at Y = D.
%
%      = 3  If the derivative of the solution with respect to Y is
%           specified at Y = C and Y = D.
%
%      = 4  If the derivative of the solution with respect to Y is
%           specified at Y = C and the solution is specified at Y = D.
%
%    BDC
%      A one dimensional array of length M that specifies the boundary
%      values of the solution at Y = C.   When NBDCND = 1 or 2,
%
%               BDC(I) = U(X(I),C) ,              I=1,2,...,M.
%
%      When NBDCND = 3 or 4,
%
%               BDC(I) = (d/dY)U(X(I),C),     I=1,2,...,M.
%
%      When NBDCND = 0, BDC is a dummy variable.
%
%    BDD
%      A one-dimensional array of length M that specifies the boundary
%      values of the solution at Y = D.  When NBDCND = 1 or 4,
%
%               BDD(I) = U(X(I),D) ,              I=1,2,...,M.
%
%      When NBDCND = 2 or 3,
%
%               BDD(I) = (d/dY)U(X(I),D) ,    I=1,2,...,M.
%
%      When NBDCND = 0, BDD is a dummy variable.
%
%    ELMBDA
%      The constant LAMBDA in the Helmholtz equation.  If LAMBDA is
%      greater than 0, a solution may not exist.  However, HSTCRT will
%      attempt to find a solution.
%
%    F
%      A two-dimensional array that specifies the values of the right
%      side of the Helmholtz equation.  For I=1,2,...,M and J=1,2,...,N
%
%               F(I,J) = F(X(I),Y(J)) .
%
%      F must be dimensioned at least M X N.
%
%    IDIMF
%      The row (or first) dimension of the array F as it appears in the
%      program calling HSTCRT.  This parameter is used to specify the
%      variable dimension of F.  IDIMF must be at least M.
%
%    W
%      A one-dimensional array that must be provided by the user for
%      work space.  W may require up to 13M + 4N + M*INT(log2(N))
%      locations.  The actual number of locations used is computed by
%      HSTCRT and is returned in the location W(1).
%
%
%             * * * * * *   On Output   * * * * * *
%
%    F
%      contains the solution U(I,J) of the finite difference
%      approximation for the grid point (X(I),Y(J)) for
%      I=1,2,...,M, J=1,2,...,N.
%
%    PERTRB
%      If a combination of periodic or derivative boundary conditions is
%      specified for a Poisson equation (LAMBDA = 0), a solution may not
%      exist.  PERTRB is a constant, calculated and subtracted from F,
%      which ensures that a solution exists.  HSTCRT then computes this
%      solution, which is a least squares solution to the original
%      approximation.  This solution plus any constant is also a
%      solution; hence, the solution is not unique.  The value of PERTRB
%      should be small compared to the right side F.  Otherwise, a
%      solution is obtained to an essentially different problem.  This
%      comparison should always be made to insure that a meaningful
%      solution has been obtained.
%
%    IERROR
%      An error flag that indicates invalid input parameters.
%       Except for numbers 0 and  6, a solution is not attempted.
%
%      =  0  No error
%
%      =  1  A .GE. B
%
%      =  2  MBDCND .LT. 0 or MBDCND .GT. 4
%
%      =  3  C .GE. D
%
%      =  4  N .LE. 2
%
%      =  5  NBDCND .LT. 0 or NBDCND .GT. 4
%
%      =  6  LAMBDA .GT. 0
%
%      =  7  IDIMF .LT. M
%
%      =  8  M .LE. 2
%
%      Since this is the only means of indicating a possibly
%      incorrect call to HSTCRT, the user should test IERROR after
%      the call.
%
%    W
%      W(1) contains the required length of W.
%
% *Long Description:
%
%     * * * * * * *   program Specifications    * * * * * * * * * * * *
%
%     Dimension of   BDA(N),BDB(N),BDC(M),BDD(M),F(IDIMF,N),
if isempty(i), i=0; end;
if isempty(id2), id2=0; end;
if isempty(id3), id3=0; end;
if isempty(id4), id4=0; end;
if isempty(ierr1), ierr1=0; end;
if isempty(j), j=0; end;
if isempty(mp), mp=0; end;
if isempty(mperod), mperod=0; end;
if isempty(np), np=0; end;
if isempty(nperod), nperod=0; end;
%     Arguments      W(See argument list)
%
%     Latest         June 1, 1977
%     Revision
%
%     Subprograms    HSTCRT,POISTG,POSTG2,GENBUN,POISD2,POISN2,POISP2,
%     Required       COSGEN,MERGE,TRIX,TRI3,PIMACH
%
%     Special        NONE
%     Conditions
%
%     Common         NONE
%     Blocks
%
%     I/O            NONE
%
%     Precision      Single
%
%     Specialist     Roland Sweet
%
%     Language       FORTRAN
%
%     History        Written by Roland Sweet at NCAR in January , 1977
%
%     Algorithm      This subroutine defines the finite-difference
%                    equations, incorporates boundary data, adjusts the
%                    right side when the system is singular and calls
%                    either POISTG or GENBUN which solves the linear
%                    system of equations.
%
%     Space          8131(decimal) = 17703(octal) locations on the
%     Required       NCAR Control Data 7600
%
%     Timing and        The execution time T on the NCAR Control Data
%     Accuracy       7600 for subroutine HSTCRT is roughly proportional
%                    to M*N*log2(N).  Some typical values are listed in
%                    the table below.
%                       The solution process employed results in a loss
%                    of no more than FOUR significant digits for N and M
%                    as large as 64.  More detailed information about
%                    accuracy can be found in the documentation for
%                    subroutine POISTG which is the routine that
%                    actually solves the finite difference equations.
%
%
%                       M(=N)    MBDCND    NBDCND    T(MSECS)
%                       -----    ------    ------    --------
%
%                        32       1-4       1-4         56
%                        64       1-4       1-4        230
%
%     Portability    American National Standards Institute Fortran.
%                    The machine dependent constant PI is defined in
%                    function PIMACH.
%
%     Required       COS
%     Resident
%     Routines
%
%     Reference      Schumann, U. and R. Sweet,'A Direct Method For
%                    The Solution Of Poisson's Equation With Neumann
%                    Boundary Conditions On A Staggered Grid Of
%                    Arbitrary Size,' J. COMP. PHYS. 20(1976),
%                    PP. 171-182.
%
%     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
%
%***REFERENCES  U. Schumann and R. Sweet, A direct method for the
%                 solution of Poisson's equation with Neumann boundary
%                 conditions on a staggered grid of arbitrary size,
%                 Journal of Computational Physics 20, (1976),
%                 pp. 171-182.
%***ROUTINES CALLED  GENBUN, POISTG
%***REVISION HISTORY  (YYMMDD)
%   801001  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  HSTCRT
%
%
f_shape=size(f);f=reshape([f(:).',zeros(1,ceil(numel(f)./prod([idimf])).*prod([idimf])-numel(f))],idimf,[]);
bda_shape=size(bda);bda=reshape(bda,1,[]);
bdb_shape=size(bdb);bdb=reshape(bdb,1,[]);
bdc_shape=size(bdc);bdc=reshape(bdc,1,[]);
bdd_shape=size(bdd);bdd=reshape(bdd,1,[]);
w_shape=size(w);w=reshape(w,1,[]);
%***FIRST EXECUTABLE STATEMENT  HSTCRT
ierror = 0;
if( a>=b )
ierror = 1;
end;
if( mbdcnd<0 || mbdcnd>4 )
ierror = 2;
end;
if( c>=d )
ierror = 3;
end;
if( n<=2 )
ierror = 4;
end;
if( nbdcnd<0 || nbdcnd>4 )
ierror = 5;
end;
if( idimf<m )
ierror = 7;
end;
if( m<=2 )
ierror = 8;
end;
if( ierror~=0 )
f_shape=zeros(f_shape);f_shape(:)=f(1:numel(f_shape));f=f_shape;
bda_shape=zeros(bda_shape);bda_shape(:)=bda(1:numel(bda_shape));bda=bda_shape;
bdb_shape=zeros(bdb_shape);bdb_shape(:)=bdb(1:numel(bdb_shape));bdb=bdb_shape;
bdc_shape=zeros(bdc_shape);bdc_shape(:)=bdc(1:numel(bdc_shape));bdc=bdc_shape;
bdd_shape=zeros(bdd_shape);bdd_shape(:)=bdd(1:numel(bdd_shape));bdd=bdd_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
return;
end;
nperod = fix(nbdcnd);
mperod = 0;
if( mbdcnd>0 )
mperod = 1;
end;
deltax =(b-a)./m;
twdelx = 1../deltax;
delxsq = 2../deltax.^2;
deltay =(d-c)./n;
twdely = 1../deltay;
delysq = deltay.^2;
twdysq = 2../delysq;
np = fix(nbdcnd + 1);
mp = fix(mbdcnd + 1);
%
%     DEFINE THE A,B,C COEFFICIENTS IN W-ARRAY.
%
id2 = fix(m);
id3 = fix(id2 + m);
id4 = fix(id3 + m);
s =(deltay./deltax).^2;
st2 = 2..*s;
for i = 1 : m;
w(i) = s;
j = fix(id2 + i);
w(j) = -st2 + elmbda.*delysq;
j = fix(id3 + i);
w(j) = s;
end; i = fix(m+1);
%
%     ENTER BOUNDARY DATA FOR X-BOUNDARIES.
%
if( mp~=1 )
if( mp==4 || mp==5 )
for j = 1 : n;
f(1,j) = f(1,j) + bda(j).*twdelx;
end; j = fix(n+1);
w(id2+1) = w(id2+1) + w(1);
else;
for j = 1 : n;
f(1,j) = f(1,j) - bda(j).*delxsq;
end; j = fix(n+1);
w(id2+1) = w(id2+1) - w(1);
end;
if( mp~=1 )
if( mp==3 || mp==4 )
for j = 1 : n;
f(m,j) = f(m,j) - bdb(j).*twdelx;
end; j = fix(n+1);
w(id3) = w(id3) + w(1);
else;
for j = 1 : n;
f(m,j) = f(m,j) - bdb(j).*delxsq;
end; j = fix(n+1);
w(id3) = w(id3) - w(1);
end;
end;
end;
%
%     ENTER BOUNDARY DATA FOR Y-BOUNDARIES.
%
if( np~=1 )
if( np==4 || np==5 )
for i = 1 : m;
f(i,1) = f(i,1) + bdc(i).*twdely;
end; i = fix(m+1);
else;
for i = 1 : m;
f(i,1) = f(i,1) - bdc(i).*twdysq;
end; i = fix(m+1);
end;
if( np~=1 )
if( np==3 || np==4 )
for i = 1 : m;
f(i,n) = f(i,n) - bdd(i).*twdely;
end; i = fix(m+1);
else;
for i = 1 : m;
f(i,n) = f(i,n) - bdd(i).*twdysq;
end; i = fix(m+1);
end;
end;
end;
for i = 1 : m;
for j = 1 : n;
f(i,j) = f(i,j).*delysq;
end; j = fix(n+1);
end; i = fix(m+1);
if( mperod~=0 )
w(1) = 0.;
w(id4) = 0.;
end;
pertrb = 0.;
if( elmbda>=0 )
if( elmbda~=0 )
ierror = 6;
elseif( mp~=2 && mp~=3 && mp~=5 ) ;
if( np~=2 && np~=3 && np~=5 )
%
%     FOR SINGULAR PROBLEMS MUST ADJUST DATA TO INSURE THAT A SOLUTION
%     WILL EXIST.
%
s = 0.;
for j = 1 : n;
for i = 1 : m;
s = s + f(i,j);
end; i = fix(m+1);
end; j = fix(n+1);
pertrb = s./(m.*n);
for j = 1 : n;
for i = 1 : m;
f(i,j) = f(i,j) - pertrb;
end; i = fix(m+1);
end; j = fix(n+1);
pertrb = pertrb./delysq;
end;
end;
end;
%
%     SOLVE THE EQUATION.
%
if( nperod==0 )
[nperod,n,mperod,m,dumvar5,dumvar6,dumvar7,idimf,f,ierr1,dumvar11]=genbun(nperod,n,mperod,m,w(sub2ind(size(w),max(1,1)):end),w(sub2ind(size(w),max(id2+1,1)):end),w(sub2ind(size(w),max(id3+1,1)):end),idimf,f,ierr1,w(sub2ind(size(w),max(id4+1,1)):end));   dumvar5i=find((w(sub2ind(size(w),max(1,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(id2+1,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(id3+1,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(id4+1,1)):end))~=(dumvar11));   w(1-1+dumvar5i)=dumvar5(dumvar5i); w(id2+1-1+dumvar6i)=dumvar6(dumvar6i); w(id3+1-1+dumvar7i)=dumvar7(dumvar7i); w(id4+1-1+dumvar11i)=dumvar11(dumvar11i);
else;
[nperod,n,mperod,m,dumvar5,dumvar6,dumvar7,idimf,f,ierr1,dumvar11]=poistg(nperod,n,mperod,m,w(sub2ind(size(w),max(1,1)):end),w(sub2ind(size(w),max(id2+1,1)):end),w(sub2ind(size(w),max(id3+1,1)):end),idimf,f,ierr1,w(sub2ind(size(w),max(id4+1,1)):end));   dumvar5i=find((w(sub2ind(size(w),max(1,1)):end))~=(dumvar5));dumvar6i=find((w(sub2ind(size(w),max(id2+1,1)):end))~=(dumvar6));dumvar7i=find((w(sub2ind(size(w),max(id3+1,1)):end))~=(dumvar7));dumvar11i=find((w(sub2ind(size(w),max(id4+1,1)):end))~=(dumvar11));   w(1-1+dumvar5i)=dumvar5(dumvar5i); w(id2+1-1+dumvar6i)=dumvar6(dumvar6i); w(id3+1-1+dumvar7i)=dumvar7(dumvar7i); w(id4+1-1+dumvar11i)=dumvar11(dumvar11i);
end;
w(1) = w(id4+1) + 3.*m;
f_shape=zeros(f_shape);f_shape(:)=f(1:numel(f_shape));f=f_shape;
bda_shape=zeros(bda_shape);bda_shape(:)=bda(1:numel(bda_shape));bda=bda_shape;
bdb_shape=zeros(bdb_shape);bdb_shape(:)=bdb(1:numel(bdb_shape));bdb=bdb_shape;
bdc_shape=zeros(bdc_shape);bdc_shape(:)=bdc(1:numel(bdc_shape));bdc=bdc_shape;
bdd_shape=zeros(bdd_shape);bdd_shape(:)=bdd(1:numel(bdd_shape));bdd=bdd_shape;
w_shape=zeros(w_shape);w_shape(:)=w(1:numel(w_shape));w=w_shape;
end %subroutine hstcrt
%DECK HSTCS1
```