SUBROUTINE SOLVE_LINEAR_SYSTEM(A, X, B, ERROR) ! Revised 2001-12-11 in order to avoid an unnecessary signal ! of out-of-index bounds in an empty sum when using index checking. 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 X (N) = M (N, N+1) / M (N, N) DO K = N-1, 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