**==mtx_solve_eqn.f  
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      SUBROUTINE  MTX_SOLVE_EQN (A,B,N,M)
C
C Pupose: MTX_SOLVE_EQN solves matrix equation for X unknowns. A is overwritten.
C	  Solution of X written into B
C
C      IMPLICIT NONE
C
      INTEGER  N, M
      REAL     A(N,N) , B(N,M)
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Local Variables:
C
      INTEGER  I, J, K, L, NEW_K
      REAL     TEMP, COEFF_MAX
C
C------------------------------------------------------------------------------
C
C========  Reduction of matrix A
C
      DO 80 K = 1 , N - 1
C
C========  Search for largest coefficient of A in first column of reduced system
C
         NEW_K = K
         COEFF_MAX = 0.0
         DO 20 I = K , N
            IF (ABS(COEFF_MAX).LT.ABS(A(I,K))) THEN
               COEFF_MAX = A(I,K)
               NEW_K = I
            END IF
 20      CONTINUE
C
C========  Exchange row number k with row number new_k , if necessary
C
         IF (NEW_K.NE.K) THEN
            DO 30 J = K , N
               TEMP = A(K,J)
               A(K,J) = A(NEW_K,J)
               A(NEW_K,J) = TEMP
 30         CONTINUE
            DO 40 J = 1 , M
               TEMP = B(K,J)
               B(K,J) = B(NEW_K,J)
               B(NEW_K,J) = TEMP
 40         CONTINUE
         END IF
C
C========  Actual reduction
C
         DO 70 I = K + 1 , N
            DO 50 J = K + 1 , N
               A(I,J) = A(I,J) - A(K,J)*A(I,K)/A(K,K)
 50         CONTINUE
            DO 60 J = 1 , M
               B(I,J) = B(I,J) - B(K,J)*A(I,K)/A(K,K)
 60         CONTINUE
 70      CONTINUE
C
 80   CONTINUE
C
C========  Backsubstitution
C
      DO 110 J = 1 , M
         B(N,J) = B(N,J)/A(N,N)
         IF (N.GT.1) THEN
            DO 100 I = N - 1 , 1 , -1
               DO 90 L = I + 1 , N
                  B(I,J) = B(I,J) - A(I,L)*B(L,J)
 90            CONTINUE
               B(I,J) = B(I,J)/A(I,I)
 100        CONTINUE
         END IF
 110  CONTINUE
C
      RETURN
      END