SUBROUTINE LINEQS(A,B,X,N,NR)
C     =============================
C
C     Solution of the linear system A*X=B
C     by Gaussian elimination with partial pivoting
C
C     Input: A  - matrix of the linear system, with actual size (N x N),
C                and maximum size (NR x NR)
C            B  - the rhs vector
C     Output: X - solution vector
C     Note that matrix A and vector B are destroyed here
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      DIMENSION A(NR,NR),B(NR),X(NR),D(MLEVEL),IP(MLEVEL)
c
      if(n.eq.2) then
         a11=a(1,1)
         a12=a(1,2)
         a21=a(2,1)
         a22=a(2,2)
         x(1)=(a(2,2)*b(1)-a(1,2)*b(2))/
     *        (a(1,1)*a(2,2)-a(1,2)*a(2,1))
         x(2)=(b(2)-a(2,1)*x(1))/a(2,2)
         return
      end if
c
      DO I=1,N
         DO J=1,N
            D(J)=A(J,I)
         END DO
         IM1=I-1
         IF(IM1.GE.1) THEN
            DO J=1,IM1
               IT=IP(J)
               A(J,I)=D(IT)
               D(IT)=D(J)
               JP1=J+1
               DO K=JP1,N
                  D(K)=D(K)-A(K,J)*A(J,I)
               END DO
            END DO
         END IF
         AM=ABS(D(I))
         IP(I)=I
         DO K=I,N
            IF(AM.LT.ABS(D(K))) THEN
               IP(I)=K
               AM=ABS(D(K))
            END IF
         END DO
         IT=IP(I)
         A(I,I)=D(IT)
         D(IT)=D(I)
         IP1=I+1
         IF(IP1.GT.N) GO TO 10
         DO K=IP1,N
            A(K,I)=D(K)/A(I,I)
         END DO
      END DO
C
   10 CONTINUE
      DO I=1,N
         IT=IP(I)
         X(I)=B(IT)
         B(IT)=B(I)
         IP1=I+1
         IF(IP1.GT.N) GO TO 20
         DO J=IP1,N
            B(J)=B(J)-A(J,I)*X(I)
         END DO
      END DO
C
   20 CONTINUE
      DO I=1,N
         K=N-I+1
         SUM=0.
         KP1=K+1
         IF(KP1.LE.N) THEN
            DO J=KP1,N
               SUM=SUM+A(K,J)*X(J)
            END DO
         END IF
         X(K)=(X(K)-SUM)/A(K,K)
      END DO
      RETURN
      END