% FIGURE 1.2.1
      function [X,IPERM,A] = DLINEQ(A,N,B)
%
%  FUNCTION DLINEQ SOLVES THE LINEAR SYSTEM A*X=B
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    A      - THE N BY N COEFFICIENT MATRIX.    THE DIAGONAL AND UPPER
%                                               TRIANGLE OF A CONTAINS U
%                                               AND THE LOWER TRIANGLE
%                                               OF A CONTAINS THE LOWER
%                                               TRIANGLE OF L, WHERE
%                                               PA = LU, P BEING THE
%                                               PERMUTATION MATRIX
%                                               DEFINED BY IPERM.
%
%    N      - THE SIZE OF MATRIX A.
%
%    X      -                                   AN N-VECTOR CONTAINING
%                                               THE SOLUTION.
%
%    B      - THE RIGHT HAND SIDE N-VECTOR.
%
%    IPERM  -                                   AN N-VECTOR CONTAINING
%                                               A RECORD OF THE ROW
%                                               INTERCHANGES MADE.  IF
%                                               J = IPERM(K), THEN ROW
%                                               J ENDED UP AS THE K-TH
%                                               ROW.
%
%-----------------------------------------------------------------------
%                              INITIALIZE IPERM = (1,2,3,...,N)
      for K=1:N
         IPERM(K) = K;
      end
%                              BEGIN FORWARD ELIMINATION
      for I=1:N-1
%                              SEARCH FROM A(I,I) ON DOWN FOR
%                              LARGEST POTENTIAL PIVOT, A(L,I)
         BIG = abs(A(I,I));
         L = I;
         for J=I+1:N
            if (abs(A(J,I)) > BIG)
               BIG = abs(A(J,I));
               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=1:N
            TEMP = A(L,K);
            A(L,K) = A(I,K);
            A(I,K) = TEMP;
         end
%                              SWITCH B(I) AND B(L)
         TEMP = B(L);
         B(L) = B(I);
         B(I) = TEMP;
%                              SWITCH IPERM(I) AND IPERM(L)
         ITEMP = IPERM(L);
         IPERM(L) = IPERM(I);
         IPERM(I) = ITEMP;
         for J=I+1:N
%                              CHOOSE MULTIPLIER TO ZERO A(J,I)
            LJI = A(J,I)/A(I,I);
            if (LJI ~= 0.0)
%                              SUBTRACT LJI TIMES ROW I FROM ROW J
               for K=I+1:N
                  A(J,K) = A(J,K) - LJI*A(I,K);
               end
%                              SUBTRACT LJI TIMES B(I) FROM B(J)
               B(J) = B(J) - LJI*B(I);
            end
%                              SAVE LJI IN A(J,I).  IT IS UNDERSTOOD,
%                              HOWEVER, THAT A(J,I) IS REALLY ZERO.
            A(J,I) = LJI;
         end
      end
      if (A(N,N) == 0.0)
         error('***** THE MATRIX IS SINGULAR *****')
      end
%                              SOLVE U*X = B USING BACK SUBSTITUTION.
      X(N) = B(N)/A(N,N);
      for I=N-1:-1:1
         SUM = 0.0;
         for J=I+1:N
            SUM = SUM + A(I,J)*X(J);
         end
         X(I) = (B(I)-SUM)/A(I,I);
      end