% FIGURE 1.5.4
      function X = DBAND(A,N,NLD,NUD,B)
%
%  FUNCTION DBAND SOLVES THE LINEAR SYSTEM A*X=B, WHERE A IS A
%    BAND MATRIX.
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    A      - THE N BY NUD+2*NLD+1 BAND MATRIX.
%             FIRST NLD COLUMNS = SUBDIAGONALS
%             NEXT  1           = MAIN DIAGONAL
%             NEXT  NUD         = SUPERDIAGONALS
%             NEXT  NLD         = WORKSPACE (FILL-IN)
%
%    N      - THE SIZE OF MATRIX A.
%
%    NLD    - NUMBER OF NONZERO LOWER DIAGONALS
%             IN A, I.E., NUMBER OF DIAGONALS
%             BELOW THE MAIN DIAGONAL.
%
%    NUD    - NUMBER OF NONZERO UPPER DIAGONALS
%             IN A, I.E., NUMBER OF DIAGONALS
%             ABOVE THE MAIN DIAGONAL.
%
%    X      -                                   AN N-VECTOR CONTAINING
%                                               THE SOLUTION.
%
%    B      - THE RIGHT HAND SIDE N-VECTOR.
%
%
%-----------------------------------------------------------------------
      MD = NLD+1;
%                              ZERO TOP NLD DIAGONALS (WORKSPACE)
      for I=1:N
      for J=NUD+1:NUD+NLD
         A(I,MD+J) = 0.0;
      end
      end
%                              BEGIN FORWARD ELIMINATION
      for I=1:N-1
%                              SEARCH FROM AII ON DOWN FOR
%                              LARGEST POTENTIAL PIVOT, ALI
         BIG = abs(A(I,MD));
         L = I;
         for J=I+1:min(I+NLD,N);
            if (abs(A(J,MD+I-J)) > BIG)
               BIG = abs(A(J,MD+I-J));
               L = J;
            end
         end
%                              IF LARGEST POTENTIAL PIVOT IS ZERO, 
%                              MATRIX IS SINGULAR
         if (BIG == 0.0)
            error('***** THE MATRIX IS SINGULAR *****')
         end
%                              SWITCH ROW I WITH ROW L, TO BRING
%                              UP LARGEST PIVOT
         for K=I:min(I+NUD+NLD,N)
            TEMP = A(L,MD+K-L);
            A(L,MD+K-L) = A(I,MD+K-I);
            A(I,MD+K-I) = TEMP;
         end
%                              SWITCH B(I) AND B(L)
         TEMP = B(L);
         B(L) = B(I);
         B(I) = TEMP;
         for J=I+1:min(I+NLD,N)
%                              CHOOSE MULTIPLIER TO ZERO AJI
            LJI = A(J,MD+I-J)/A(I,MD);
            if (LJI ~= 0.0)
%                              SUBTRACT LJI TIMES ROW I FROM ROW J
               for K=I:min(I+NUD+NLD,N)
                  A(J,MD+K-J) = A(J,MD+K-J) - LJI*A(I,MD+K-I);
               end
%                              SUBTRACT LJI TIMES B(I) FROM B(J)
               B(J) = B(J) - LJI*B(I);
            end
         end
      end
      if (A(N,MD) == 0.0)
         error('***** THE MATRIX IS SINGULAR *****')
      end
%                              SOLVE U*X = B USING BACK SUBSTITUTION.
      X(N) = B(N)/A(N,MD);
      for I=N-1:-1:1
         SUM = 0.0;
         for J=I+1:min(I+NUD+NLD,N)
            SUM = SUM + A(I,MD+J-I)*X(J);
         end
         X(I) = (B(I)-SUM)/A(I,MD);
      end