Документ взят из кэша поисковой машины. Адрес оригинального документа : http://crydee.sai.msu.ru/f90/code/solve1.f90 Дата изменения: Tue Jan 9 15:32:08 1996 Дата индексирования: Mon Oct 1 21:51:53 2012 Кодировка: Поисковые слова: vallis

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

! The work area M is A extended with B
REAL, DIMENSION (SIZE (B), SIZE (B) + 1) :: M
INTEGER, DIMENSION (1) :: MAX_LOC
REAL, 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) 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) ERROR = .TRUE. ! Singular matrix A

! Back substitution
IF (ERROR) THEN
X = 0.0
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