'***********************************************************
 '*           Test Gauss method for band matrix             *
 '*                                                         *
 '*                    Basic version by J-P Moreau, Paris   *
 '*                          (in simple precision)          *
 '*                            (www.jpmoreau.fr)            *
 '* ------------------------------------------------------- *
 '* Reference for original C version:                       *
 '*                                                         *
 '* "Numerical Algorithms with C, By Gisela Engeln-Muellges *
 '*  and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].    *
 '* ------------------------------------------------------- *
 '* SAMPLE RUN:                                             *
 '*                                                         *
 '* Input file: bandmat.dat                                 *
 '*                                                         *
 '* 5                                                       *
 '* 1                                                       *
 '* 1                                                       *
 '*  0.0   1.0   -2.0                                       *
 '* 10.0   5.0   -7.5                                       *
 '*  2.22 -3.25   0.00456                                   *
 '*  2.0   4.0   -3.3                                       *
 '* -1.0   3.0    0.0                                       *
 '*                                                         *
 '* Output file: bandmat.lst                                *
 '*                                                         *
 '* --------------------------------------------------      *
 '*  Banded matrix                                          *
 '* --------------------------------------------------      *
 '* Dimension     :  5                                      *
 '* Subdiagonals  :  1                                      *
 '* Superdiagonals:  1                                      *
 '*                                                         *
 '* Input condensed band matrix:                            *
 '*  0.000000  1.000000 -2.000000                           *
 '* 10.000000  5.000000 -7.500000                           *
 '*  2.220000 -3.250000  0.004560                           *
 '*  2.000000  4.000000 -3.300000                           *
 '* -1.000000  3.000000  0.000000                           *
 '* Input uncondensed band matrix:                          *
 '*  1.000000 -2.000000  0.000000  0.000000  0.000000       *
 '* 10.000000  5.000000 -7.500000  0.000000  0.000000       *
 '*  0.000000  2.220000 -3.250000  0.004560  0.000000       *
 '*  0.000000  0.000000  2.000000  4.000000 -3.300000       *
 '*  0.000000  0.000000  0.000000 -1.000000  3.000000       *
 '* --------------------------------------------------      *
 '* Band without pivot                                      *
 '* --------------------------------------------------      *
 '* Inverse (transposed):                                   *
 '* -0.005941 -0.502971 -0.343235  0.236714  0.078905       *
 '*  0.100594  0.050297  0.034324 -0.023671 -0.007890       *
 '* -0.231916 -0.115958 -0.386526  0.266570  0.088857       *
 '*  0.000365  0.000182  0.000608  0.344408  0.114803       *
 '*  0.000401  0.000201  0.000669  0.378849  0.459616       *
 '* Determinant = -562.7040405273438                        *
 '* --------------------------------------------------      *
 '* Band with pivot                                         *
 '* --------------------------------------------------      *
 '* Inverse (transposed):                                   *
 '* -0.005924 -0.503080 -0.343285  0.236748  0.078916       *
 '*  0.100592  0.050308  0.034328 -0.023675 -0.007892       *
 '* -0.231375 -0.115820 -0.385713  0.266009  0.088670       *
 '*  0.000370  0.000013  0.000501  0.344482  0.114827       *
 '*  0.000407  0.000014  0.000552  0.378930  0.459643       *
 '* Determinant =  562.7040405273438                        *
 '* System Solution:                                        *
 '* -0.135929 -0.568566 -0.693616  1.202494  0.734165       *
 '*                                                         *
 '* --------------------------------------------------      *
 '* End of file bandmat.lst                                 *
 '*                                                         *
 '* Note: in basic, the results for the method with pivot   *
 '*       (bandsol) are slightly erroneous.                 *
 '***********************************************************

'-----------------------------------------------------------------
' Nota: Matrices (n,n) are put in line in a vector of size n*n
'                        1 2 3
'       Example:  Matrix 4 5 6 is stored as (1,2,3,4,5,6,7,8,9)
'                        7 8 9
'-----------------------------------------------------------------
' Global variables:
'  iud, ld, n, i, k, icode, isign, idim, mode : integer
'  a            : vector(n*n)  storing a real matrix n*n
'  packmat,save : vector(n*n)  storing a real matrix n*idim
'  b            : real vector(n)
'  iperm        : integer vector(n)
'  determ       : REAL   {matrix determinant}
'  solution, xmatinv : vectors(n*n) storing real matrices n*n

defint i-n

'open input and output files

  open "bandmat.dat" for input as #1
  open "bandmat.lst" for output as #2

  cls
  MACHEPS=1e-12

  print #2, "--------------------------------------------------"
  print #2, "  Banded matrix"
  print #2, "--------------------------------------------------"

  input #1, n

  if n < 1 then
    nom$=" n must be > 0"
    gosub 5000
  end if

  dim a(n*n),b(n),iperm(n)

  input #1, ld
  if ld < 0 then
    nom$=" ld must be > 0"
    gosub 5000
  end if

  input #1, iud
  if iud < 0 then
    nom$=" ud must be > 0"
    gosub 5000
  end if

  print #2, " Dimension     : "; n
  print #2, " Subdiagonals  : "; ld
  print #2, " Superdiagonals: "; iud
  print #2,

  if ld + iud >= n then
    nom$=" ld + ud must be < n"
    gosub 5000
  end if

  if ld<iud then
    min=ld
  else
    min=iud
  end if

  idim = ld + iud + 1 + min

  dim packmat(n*n), save(n*n), solution(n*n), tmp(n*n), xmatinv(n*n)

  'Set vector packmat to zero
  for i=0 to n*(ld+iud+1)-1
    packmat(i)=0#
  next i

  'read input matrix in condensed form and save in "save" vector
  for i=0 to n*(ld+iud+1)-1
    input #1, packmat(i)
    save(i)=packmat(i)
  next i
  'close input file
  close #1

  print #2, " Input condensed band matrix:"

  'print packmat, ld+iud+1 items per line
  j=0
  for i=0 to n*(ld+iud+1)-1
    j=j+1
    print #2, using " ##.######"; packmat(i);
    if j=ld+iud+1 then
      j=0
      print #2,
    end if
  next i

  'Set a to zero
  for i=0 to n*n-1
    a(i)=0#
  next i

  for i=0 to n-1
    for k=0 to n-1
      if ld+k-i>-1 then
        if ld+k-i<ld+iud+1 then
          a(i*n+k)=packmat(i*(ld+iud+1)+ld+k-i)
        else
          a(i*n+k)=0#
        end if
      end if
    next k
  next i

  print #2, " Input uncondensed band matrix:"

  'print a matrix, n items per line
  j=0
  for i=0 to n*n-1
    j=j+1
    print #2, using " ##.######"; a(i);
    if j=n then
      j=0
      print #2,
    end if
  next i

  print #2, "--------------------------------------------------"
  print #2, " Band without pivot"
  print #2, "--------------------------------------------------"
  print #2, " Inverse (transposed):"

  mode = 0
  for ii = 0 to n-1
    'Set b vector to zero
    for j=0 to n-1
      b(j)=0#
    next j
    b(ii) = 1#
    'call band method without pivot
    'bando (mode, n, ld, iud, packmat, b, icode)
    gosub 1000
    if icode<>0 then
      nom$= " bando return code<>0"
      gosub 5000
    end if
    'print b vector, n items per line
    for j=0 to n-1
      print #2, using " ##.######";b(j);
    next j
    print #2,
    mode = 2
  next ii

  determ = 1#
  for i = 0 to n-1
    determ = determ * packmat(i*(ld+iud+1)+ld)
  next i

  print #2, " Determinant = "; determ

  print #2, "--------------------------------------------------"
  print #2, " Band with pivot"
  print #2, "--------------------------------------------------"
  print #2, " Inverse (transposed):"

  mode = 0
  for ii = 0 to n-1
    'Set b vector to zero
    for j=0 to n-1
      b(j)=0#
    next j
    b(ii) = 1#
    'call band method with pivot
    'band (mode, n, ld, iud, save, b, iperm, isign, icode)
    gosub 2000

    if icode<>0 then
      nom$=" band return code <> 0"
      gosub 5000
    end if
    'print b vector, n items per line
    for j=0 to n-1
      print #2, using " ##.######";b(j);
    next j
    print #2,
    solution(ii)=0#
    for k=0 to n-1
      xmatinv(ii*n+k)=b(k)
    next k
    mode = 2
  next ii

  determ = 1#*isign
  for i = 0 to n-1
    determ = determ * packmat(i*(ld+iud+1)+ld)
  next i

  print #2, " Determinant = "; determ

  print #2, " System Solution:"
  for k=0 to n-1
    for i=0 to n-1
      solution(k) = solution(k) + xmatinv(i*n+k)
    next i
  next k

  'print solution, n items per line
  j=0
  for i=0 to n-1
    j=j+1
    print #2, using " ##.######"; solution(i);
    if j=n then
      j=0
      print #2,
    end if
  next i

  print #2,
  print #2, "--------------------------------------------------"
  print #2, "End of file bandmat.lst
  close #2

  print
  print " Results in bandmat.lst."

END 'of main program

1000 'bando subroutine
'======================================================================
'*                                                                    *
'* The procedure bando solves a linear banded system:  packmat * x = b*
'* Here packmat is a nonsingular n x n matrix in condensed form, i.e. *
'* represented in an ld+ud+1 x n matrix for its ld lower and ud upper *
'* co-diagonals. b denotes the right hand side of the system, and x   *
'* is the solution.                                                   *
'*                                                                    *
'* bando uses the Gauss algorithm without column pivot search.        *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Applications:                                                    *
'*   =============                                                    *
'*      Solve linear systems with nonsingular banded system matrices. *
'*      Particularly useful for large sparse and banded and diagonally*
'*      dominant matrices.                                            *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Control parameter:                                               *
'*   ==================                                               *
'*      mode     integer mode                                         *
'*               calling modus for band:                              *
'*       = 0     factor matrix and solve linear system                *
'*       = 1     factor only, store factors in packmat                *
'*       = 2     solve linear system only; for this call, the factors *
'*               must already be available in packmat, such as when   *
'*               many systems are solved for differing right hand     *
'*               sides and the same system matrix.                    *
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer n  ( n > 2 )                                 *
'*               Dimension of packmat(n*n), size of b(n)              *
'*      ld       integer ld; ( ld >= 0 )                              *
'*               number of lower co-diagonals                         *
'*      iud      integer iud ( iud >= 0 )                             *
'*               number of upper co-diagonals                         *
'*      packmat  REAL packmat(n,n) stored in a vector of size n*n     *
'*               mode = 0, 1:                                         *
'*               Matrix of the system in comndensed form.             *
'*               Each row has length at least ld + 1 + ud + min(ld,ud)*
'*               where the columns 0, .., ld-1 denote the lower       *
'*               co-diagonals, column ld stores the diagonal and the  *
'*               columns ld+1, .., ld+ud contain the upper            *
'*               co-diagonals.                                        *
'*               If A is the original uncondensed band matrix, then : *
'*               A[i][k] = packmat[i][ld+k-i],                        *
'*                           for k,i inside the band                  *
'*               mode = 2:                                            *
'*               LU factors in condensed form.                        *
'*      b        REAL vector b(n)         ( for mode = 0, 2 )         *
'*               Right hand side                                      *
'*                                                                    *
'*   Output parameters:                                               *
'*   ==================                                               *
'*      packmat  REAL packmat(n,n);       ( for mode = 0, 1 )         *
'*               LU factorization in condensed form                   *
'*      b        REAL vector b(n)         ( for mode = 0, 2 )         *
'*               solution vector for the system                       *
'*                                                                    *
'*   Return value icode:                                              *
'*   ==================                                               *
'*      = 0      all ok                                               *
'*      = 1      n < 3 or other incorrect input parameter             *
'*      = 2      lack of space (not used here)                        *
'*      = 3      Matrix is numerically singular                       *
'*      = 4      wrong value for mode                                 *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Subroutines called:                                              *
'*   ==================                                               *
'*      1500 banodec: Factor matrix                                   *
'*      1750 banosol: solve linear system                             *
'*                                                                    *
'*=====================================================================
  'wrong parameters
  if n < 3 or ld < 0 or iud < 0 or n < ld + 1 + iud then
    icode = 1
    return
  end if

  if mode=0 then  'factor and solve system
    'banodec (n, ld, iud, packmat, icode)
    gosub 1500
    if icode = 0 then
      'banosol (n, ld, iud, packmat, b, icode)
      gosub 1750
    end if
  elseif mode=1 then  'factor only
    'banodec (n, ld, iud, packmat, icode)
    gosub 1500
  elseif mode=2 then  'solve only
    'banosol (n, ld, iud, packmat, b, icode)
    gosub 1750
  else
    icode = 3
  end if
RETURN  'bando subroutine


1500 'banodec subroutine
'======================================================================
'*                                                                    *
'*  The procedure banodec factors a condensed banded matrix packmat.  *
'*  banddec uses the Gauss algorithm without column pivot search.     *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer n;  ( n > 2 )                                *
'*               Dimension of packmat(n*n), size of b(n)              *
'*      ld       integer ld; ( ld >= 0 )                              *
'*               number of lower co-diagonals                         *
'*      iud      integer iud; ( iud >= 0 )                            *
'*               number of upper co-diagonals                         *
'*      packmat  REAL packmat(n,n) stored in a vector of size n*n     *
'*               Matrix of the system in comndensed form.             *
'*               Each row has length at least ld + 1 + ud + min(ld,ud)*
'*               where the columns 0, .., ld-1 denote the lower       *
'*               co-diagonals, column ld stores the diagonal and the  *
'*               columns ld+1, .., ld+ud contain the upper            *
'*               co-diagonals.                                        *
'*                                                                    *
'*   Output parameters:                                               *
'*   ==================                                               *
'*      packmat  REAL packmat(n,n) stored in a vector of size n*n     *
'*               LU factorization in condensed form                   *
'*                                                                    *
'*   Return value code:                                               *
'*   =================                                                *
'*      = 0      all ok                                               *
'*      = 1      n < 3 or other incorrect input parameter             *
'*      = 2      Matrix is numerically singular                       *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Constants used :     MACHEPS                                     *
'*   ================                                                 *
'*                                                                    *
'*=====================================================================
' kend, kjend, jm, jk, k, j, i, m : INTEGER;
  m = iud + ld + 1                   'dimension #2 of packmat

  if ld = 0 then                     'Matrix already in upper
    icode=0                          'triangular form
    return
  end if

  for i = 0 to n - 2                 'loop over all rows
    if ld+1 < n-i then
      imin=ld+1
    else
      imin=n-i
    end if
    kend  = imin
    if iud+1 < n-i then
      imin=iud+1
    else
      imin=n-i
    end if
    kjend = imin

    if ABS(packmat(i*m+ld)) < MACHEPS then    'LU decompsition does
      icode=2                                 'not exist
      return
    end if

    for k=1 to kend-1                       'loop over all rows below row i
      packmat((k+i)*m+ld-k) = packmat((k+i)*m+ld-k) / packmat(i*m+ld)

      for j = 1 to kjend-1
        jk = j + ld - k
        jm = j + ld
        packmat((k+i)*m+jk) = packmat((k+i)*m+jk) - packmat((k+i)*m+ld-k) * packmat(i*m+jm)
      next j
    next k   'k loop
  next i  'i loop

  icode = 0
RETURN


1750 'banosol subroutine
'======================================================================
'*                                                                    *
'*  The procedure banosol solves a factored linear banded system in   *
'*  condensed form using banodec.                                     *
'*                                                                    *
'*====================================================================*
'*                                                                    *
'*   Input parameters:                                                *
'*   ================                                                 *
'*      n        integer n;  ( n > 2 )                                *
'*               Dimension of packmat(n*n), size of b(n)              *
'*      ld       integer ld; ( ld >= 0 )                              *
'*               number of lower co-diagonals                         *
'*      iud      integer iud ( iud >= 0 )                             *
'*               number of upper co-diagonals                         *
'*      packmat  REAL packmat(n,n) stored in a vector of size n*n     *
'*               Matrices for the factored system in comndensed form. *
'*                                                                    *
'*   Output parameters:                                               *
'*   ==================                                               *
'*      b        REAL b(n)                                            *
'*               solution vector of linear system                     *
'*                                                                    *
'*   Return value icode:                                              *
'*   ==================                                               *
'*      = 0      all ok                                               *
'*      = 1      n < 3 or other incorrect input parameter             *
'*                                                                    *
'*====================================================================*
'  kend, i, k, m : INTEGER;
  m = ld + iud + 1
                     'Invalid input parameter
  if n < 3 or ld < 0 or iud < 0 or m > n then
    icode = 1
    return
  end if

  for i = 0 to n - 2
    if ld+1 < n-i then
      imin=ld+1
    else
      imin=n-i
    end if
    kend = imin
    for k = 1 to kend-1
      b(k+i) = b(k+i) - packmat((k+i)*m+ld-k) * b(i)
    next k
  next i

  for i = n - 1 to 0 step -1     'back substitution
    if iud+1 < n-i then
      imin=iud+1
    else
      imin=n-i
    end if
    kend = imin
    for k = 1 to kend-1
      b(i) = b(i) - packmat(i*m+k+ld) * b(i+k)
    next k
    b(i) = b(i) / packmat(i*m+ld)
  next i

  icode = 0
RETURN  'banosol subroutine

2000 'band subroutine
'======================================================================
'*                                                                    *
'*  The procedure band solves a linear banded system:  save * x = b.  *
'*  Here save is a nonsingular n x n matrix in condensed form, i.e.   *
'*  represented in an ld+ud+1 x n matrix for its ld lower and ud upper*
'*  co-diagonals. b denotes the right hand side of the system, and x  *
'*  is the solution.                                                  *
'*                                                                    *
'*  band uses the Gauss algorithm with column pivot search.           *
'*  The result of pivoting are min( ud, ld) additional columns,       *
'*  so that save needs all in all a n x (ld+1+ud+min(ld,ud)) matrix.  *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Applications:                                                    *
'*   -------------                                                    *
'*      Solve linear systems with nonsingular banded system matrices. *
'*      Particularly useful for large sparse and banded matrices with *
'*      n >> ld+1+ud.                                                 *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Control parameter:                                               *
'*   ------------------                                               *
'*      mode     mode:INTEGER;                                        *
'*               calling modus for band:                              *
'*       = 0     factor matrix and solve linear system                *
'*       = 1     factor only, store factors in save                   *
'*       = 2     solve linear system only; for this call, the factors *
'*               must already be available in save, such as when      *
'*               many systems are solved for differing right hand     *
'*               sides and the same system matrix.                    *
'*                                                                    *
'*   Input parameters:                                                *
'*   -----------------                                                *
'*      n        integer n  ( n > 2 )                                 *
'*               Dimension of save, size of b                         *
'*      ld       integer ld ( ld >= 0 )                               *
'*               number of lower co-diagonals                         *
'*      iud      integer iud ( iud >= 0 )                             *
'*               number of upper co-diagonals                         *
'*      save     REAL  save(n,n) stored in a vector of size n*n       *
'*               mode = 0, 1:                                         *
'*               Matrix of the system in condensed form.              *
'*               Each row has length at least ld + 1 + ud + min(ld,ud)*
'*               where the columns 0, .., ld-1 denote the lower       *
'*               co-diagonals, column ld stores the diagonal and the  *
'*               columns ld+1, .., ld+ud contain the upper            *
'*               co-diagonals.                                        *
'*               If A is the original uncondensed band matrix, then:  *
'*               A[i][k] = save[i][ld+k-i],                           *
'*                         for k,i inside the band                    *
'*               mode = 2:                                            *
'*               LU factors in condensed form.                        *
'*      b        REAL vector b(n)         ( for mode = 0, 2 )         *
'*               Right hand side                                      *
'*      iperm    integer vector iperm(n)  ( for mode = 2 )            *
'*               row permutation vector                               *
'*      isign    integer isign;           ( for mode = 2 )            *
'*               sign of iperm. The determinant of A can be computed  *
'*               as the product of the diagonal entries times isign.  *
'*                                                                    *
'*   Output parameters:                                               *
'*   ------------------                                               *
'*      save     REAL save(n,n)        ( for mode = 0, 1 )            *
'*               LU factorization in condensed form                   *
'*      iperm    integer iperm(n)      ( for mode = 0, 1 )            *
'*               row permutation vector                               *
'*      b        REAL   b(n)           ( for mode = 0, 2 )            *
'*               solution vector for the system                       *
'*      isign    integer isign;        ( for mode = 0, 1 )            *
'*               sign of iperm                                        *
'*                                                                    *
'*   Return value code:                                               *
'*   -----------------                                                *
'*      = 0      all ok                                               *
'*      = 1      n < 3 or other incorrect input parameter             *
'*      = 2      lack of space ( not used here )                      *
'*      = 3      Matrix is numerically singular                       *
'*      = 4      wrong calling modus                                  *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Subroutines called:                                              *
'*   ------------------                                               *
'*      3000 banddec: Factor matrix                                   *
'*      4000 bandsol: solve linear system                             *
'*                                                                    *
'======================================================================
  if n < 1 or ld < 0 or iud < 0 then  'wrong parameters
    icode=1
    return
  end if

  if mode=0 then  'factor and solve system
    'banddec (n, ld, ud, save, iperm, isign,icode)
    gosub 3000
    if icode = 0 then
      'bandsol (n, ld, ud, save, b, iperm, icode)
      gosub 4000
    end if
  elseif mode=1 then  'factor only
    'banddec (n, ld, ud, save, iperm, isign, icode)
    gosub 3000
  elseif mode=2 then  'solve only
    'bandsol (n, ld, ud, save, b, iperm, icode)
    gosub 4000
  else
    icode = 4
  end if

RETURN 'band subroutine

3000 'banddec subroutine
'======================================================================
'*                                                                    *
'*  The procedure banddec factors a condensed banded matrix save.     *
'*  banddec uses the Gauss algorithm with column pivot search.        *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Input parameters:                                                *
'*   ----------------                                                 *
'*      n        integer n;  ( n > 2 )                                *
'*               Dimension of save, size of b                         *
'*      ld       integer ld; ( ld >= 0 )                              *
'*               number of lower co-diagonals                         *
'*      ud       integer ud; ( ud >= 0 )                              *
'*               number of upper co-diagonals                         *
'*      save     REAL save(n,n) stored in a vector of size n*n        *
'*               Matrix of the system in comndensed form.             *
'*               Each row has length at least ld + 1 + ud + min(ld,ud)*
'*               where the columns 0, .., ld-1 denote the lower       *
'*               co-diagonals, column ld stores the diagonal and the  *
'*               columns ld+1, .., ld+ud contain the upper            *
'*               co-diagonals.                                        *
'*      iperm    integer vector iperm(n)                              *
'*               row permutation vector                               *
'*      isign    integer isign;                                       *
'*               sign of iperm. The determinant of A can be computed  *
'*               as the product of the diagonal entries times isign.  *
'*                                                                    *
'*   Output parameters:                                               *
'*   ------------------                                               *
'*      save     REAL save(n,n) stored in a vector of size n*n        *
'*               LU factorization in condensed form                   *
'*      iperm    integer vector iperm(n)                              *
'*               row permutation vector                               *
'*      isign    integer isign;                                       *
'*               sign of iperm                                        *
'*                                                                    *
'*   Return value code:                                               *
'*   -----------------                                                *
'*      = 0      all ok                                               *
'*      = 1      n < 3 or other incorrect input parameter             *
'*      = 2      not used                                             *
'*      = 3      Matrix is numerically singular                       *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Constants used :     MACHEPS                                     *
'*   ----------------                                                 *
'*                                                                    *
'======================================================================
' j0, mm, iup, istart, iend, istep, kstart,kend, kjend, km, jm, jk:INTEGER
' k, j, i, i1, m : INTEGER
' piv : REAL
' tmp : vector(n)

  m = ld + iud + 1                      'second dimension of original save

  if ld < 0 or iud < 0 or n < 1 then
    icode=1                             'Invalid input parameter
    return
  end if

  mm = ld + 1 + iud + min               'second dimension of extended save

  if ld <= iud then
    iup=1                               'iup = 0 ==> transform into
  else                                  'lower triangular matrix
    iup=0
  end if

  for i = 0 to n-1
  'initialize needed extra columns using a temporary matrix
    for k=0 to mm-1
      if k<m then
        tmp(i*mm+k)=save(i*m+k)
      else
        tmp(i*mm+k)=0#
      end if
    next k
  next i

  'copy tmp to save
  for i=0 to n*mm-1
    save(i)=tmp(i)
  next i

  isign = 1                                  'initialize isign
  if iup=1 then                              'Find start, end and
    istart = 0 : iend = n-1 : istep = 1      'directions of loops
    kstart=1                                 'depending of iup
  else
    istart = n-1 : iend = 0 : istep = -1
    kstart = -1
  end if

  for i=istart to iend step istep
    'C code: kend = (up ? min (ld+1, n-i) : max (-i-1, -ud-1));
    if ld+1<n-i then
      imin=ld+1
    else
      imin=n-i
    end if
    if -i-1 > -iud-1 then
      imax=-i-1
    else
      imax=-iud-1
    end if
    if iup = 1 then
      kend = imin
    else
      kend = imax
    end if
    j0 = 0
    piv = ABS (save(i*mm+ld))                'choose pivot

    for k=kstart to kend step istep
      if ABS(save((k+i)*mm+ld-k)) > piv  then
        piv = ABS(save((k+i)*mm+ld-k))
        j0 = k
      end if
    next k

    'piv is too small, matrix is singular
    if piv < MACHEPS then
      icode = 3
      print #2, " BANDDEC - code="; icode
      close #2
      end
    end if

    iperm(i) = j0
    'C code: kjend :=  (up ? min (j0+ud+1, n-i) : max ( -i-1, j0-ld-1));
    if j0+iud+1<n-i then
      imin=j0+iud+1
    else
      imin=n-i
    end if
    if -i-1 > j0-ld-1 then
      imax=-i-1
    else
      imax=j0-ld-1
    end if
    if iup=1 then
      kjend = imin
    else
      kjend = imax
    end if

    if j0 <> 0 then
      isign = -isign                            'swap rows
      k = 0
3200  km = k + ld
        if km < 0 then km = km + mm
        tt=save(i*mm+km)
        save(i*mm+km)=save((i+j0)*mm+k+ld-j0)
        save((i+j0)*mm+k+ld-j0)=tt
        k = k + istep
      if k<>kjend then goto 3200
    end if

    for k=kstart to kend step istep
      save((k+i)*mm+ld-k) = save((k+i)*mm+ld-k)/ save(i*mm+ld)
      j=kstart
      for j=kstart to kjend step istep
        jk = j + ld - k                'loop over all columns right of i
        jm = j + ld
                                       'additional columns from pivoting
        if jk < 0 then jk = jk + mm    'are stored to right of column
        if jm < 0 then jm = jm + mm    'ud+ld+1

        save((k+i)*mm+jk) = save((k+i)*mm+jk) - save((k+i)*mm+ld-k)*save(i*mm+jm)

      next j
    next k
  next i

  if iup=1 then
    i1 = n-1
  else
    i1 = 0
  end if

  piv = ABS(save(i1*mm+ld))                    ' choose pivot
                                               'If piv = 0, matrix is
  if piv < MACHEPS then                        'singular
    icode=3
    print #2, " BANDDEC - code="; icode
    close #2
    end
  end if
  iperm(iend) = 0
  icode = 0
RETURN  'banddec subroutine


4000 'bandsol subroutine
'======================================================================
'*                                                                    *
'*  The procedure bandsol solves a factored linear banded system in   *
'*  condensed form using banddec.                                     *
'*                                                                    *
'* ------------------------------------------------------------------ *
'*                                                                    *
'*   Input parameters:                                                *
'*   -----------------                                                *
'*      n        integer n;  ( n > 2 )                                *
'*               Dimension of save, size of b                         *
'*      ld       integer ld; ( ld >= 0 )                              *
'*               number of lower co-diagonals                         *
'*      iud      integer iud ( iud >= 0 )                             *
'*               number of upper co-diagonals                         *
'*      save     REAL save(n,n) stored in a vector of size n*n        *
'*               Matrices for the factored system in comndensed form. *
'*      iperm    integer vector iperm(n)                              *
'*               row permutation vector                               *
'*                                                                    *
'*   Output parameters:                                               *
'*   ------------------                                               *
'*      b        REAL vector b(n)                                     *
'*               solution vector of linear system                     *
'*                                                                    *
'*   Return value code:                                               *
'*   -----------------                                                *
'*      = 0      all ok                                               *
'*      = 1      n < 3 or other incorrect input parameter             *
'*                                                                    *
'======================================================================
' i, k, is, mm, iup, istart, iend, m, istep, kstart, kend, km : INTEGER;

  m = ld + iud + 1                          'second dimension of save

  if ld < 0 or iud < 0 or n < 1 then        'invalid input
    icode = 1
    return
  end if

  mm = ld + iud + 1 + min               'mm = max. column number

  if ld <= iud then
    iup=1                               'iup = 0 ==> transform into
  else                                  'lower triangular matrix
    iup=0
  end if

  if iup=1 then
    istart = 0 : iend = n-1 : istep = 1    'determine bounds and
    kstart = 1 : is = -1                    'direction of loop
  else                                     'depending on iup
    istart = n-1 : iend = 0 : istep = -1
    kstart = -1 : is = 1
  end if

  for i=istart to iend-istep step istep
    if iperm(i) <> 0 then
      tt=b(i) : b(i)=b(i+iperm(i)) : b(i+iperm(i))=tt
    end if

    'C code: kend = (up ? min (ld+1, n-i) : max (-i-1, -ud-1));
    if ld+1<n-i then
      imin=ld+1
    else
      imin=n-i
    end if
    if -i-1 > -iud-1 then
      imax=-i-1
    else
      imax=-iud-1
    end if
    if iup=1 then
      kend = imin
    else
      kend = imax
    end if

    for k=kstart to kend-istep step istep
      b(k+i) = b(k+i) - save((k+i)*mm+ld-k) * b(i)
    next k
  next i

  for i=iend to istart+is+istep step -istep
    'C code: kend :=  (up ? min (ld+ud+1, n-i) : max (-i-1, -ud-ld));
    if ld+iud+1<n-i then
      imin=ld+iud+1
    else
      imin=n-i
    end if
    if -i-1 > -iud-ld then
      imax=-i-1
    else
      imax=-iud-ld
    end if
    if iup=1 then
      kend = imin
    else
      kend = imax
    end if

    for k=kstart to kend-istep step istep
      km = k + ld                       'update and
      if km < 0 then km = km + mm       'back substitute
      b(i) = b(i) - save(i*mm+km) * b(i+k)
    next k
    if ABS(save(i*mm+ld)) > MACHEPS then
      b(i) = b(i) / save(i*mm+ld)
    else
      b(i)=0#
    end if
  next i
  icode = 0
RETURN  'bandsol subroutine


5000 'screen error message and stop
  print
  print nom$
  print
  input R$
  end
return

'end of file tband.bas