SUBROUTINE SOLVE_LINEAR_SYSTEM(A, X, B, ERROR)
        IMPLICIT NONE
        ! Array specifications
        DOUBLE PRECISION, DIMENSION (:, :), INTENT (IN)  :: A
        DOUBLE PRECISION, DIMENSION (:),    INTENT (OUT) :: X
        DOUBLE PRECISION, DIMENSION (:),    INTENT (IN)  :: B
        LOGICAL, INTENT (OUT) 		                 :: ERROR

        ! The work area M is A extended with B
        DOUBLE PRECISION, DIMENSION (SIZE (B), SIZE (B) + 1) :: M
        INTEGER, DIMENSION (1)            		     :: MAX_LOC
        DOUBLE PRECISION, DIMENSION (SIZE (B) + 1)           :: TEMP_ROW 
        INTEGER                                  	     :: N, K

        ! Initiate M
        N = SIZE (B)
        M (1:N, 1:N) = A
        M (1:N, N+1) = B

        ! Triangularization phase
        ERROR = .FALSE.
        TRIANG_LOOP: DO K = 1, N - 1
                ! Pivotering
                MAX_LOC = MAXLOC (ABS (M (K:N, K)))
                IF ( MAX_LOC(1) /= 1 ) THEN
                    TEMP_ROW (K:N+1 ) = M (K, K:N+1)
                    M (K, K:N+1) = M (K-1+MAX_LOC(1), K:N+1)
                    M (K-1+MAX_LOC(1), K:N+1) = TEMP_ROW (K:N+1)
                END IF

                IF (M (K, K) == 0.0D0) THEN
                        ERROR  = .TRUE.     ! Singular matrix A
                        EXIT TRIANG_LOOP
                ELSE
                        TEMP_ROW (K+1:N) = M (K+1:N, K) / M (K, K)
                        M (K+1:N, K+1:N+1) = M (K+1:N, K+1:N+1)  &
                           - SPREAD( TEMP_ROW (K+1:N), 2, N-K+1) &
                           * SPREAD( M (K, K+1:N+1), 1, N-K)
                        M (K+1:N, K) = 0 ! These values are not used
                END IF
        END DO TRIANG_LOOP

        IF (M (N, N) == 0.0D0) ERROR  = .TRUE.  ! Singular matrix A

        ! Back substitution
        IF (ERROR) THEN
           X = 0.0D0
        ELSE
           DO K = N, 1, -1
              X (K) = M (K, N+1) - SUM (M (K, K+1:N) * X (K+1:N))
              X (K) = X (K) / M (K, K)
           END DO
        END IF

END SUBROUTINE SOLVE_LINEAR_SYSTEM
