SUBROUTINE LOES_LIN_EKV_SYSTEM(A, X, B, FEL) ! Reviderad 2001-12-03 för att undvika irrelevant signal ! för indexöverskridande IMPLICIT NONE ! Fältdeklarationer DOUBLE PRECISION, DIMENSION (:, :), INTENT (IN) :: A DOUBLE PRECISION, DIMENSION (:), INTENT (OUT) :: X DOUBLE PRECISION, DIMENSION (:), INTENT (IN) :: B LOGICAL, INTENT (OUT) :: FEL ! Arbetsarean M är A utvidgad med 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 ! Initiera M N = SIZE (B) M (1:N, 1:N) = A M (1:N, N+1) = B ! Triangulariseringsfas FEL = .FALSE. TRIANG_SLINGA: 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 FEL = .TRUE. ! Singulär matris A EXIT TRIANG_SLINGA 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 ! Dessa värden användes ej END IF END DO TRIANG_SLINGA IF (M (N, N) == 0.0D0) FEL = .TRUE. ! Singulär matris A ! Återsubstitution IF (FEL) 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)) ) / M (K, K) END DO END IF END SUBROUTINE LOES_LIN_EKV_SYSTEM