test1402.f90

这个程序是对高斯消去法求解线性方程组的fortran源代码

SUBROUTINE Gauss_Expunction(A, B, X, N)
  IMPLICIT NONE
! 变量定义
  INTEGER :: N
  REAL :: A(N,N), B(N)
  REAL :: X(N)
  INTEGER :: I, J
  REAL :: C
  LOGICAL, EXTERNAL :: Check_Root

  ! 消元步
  DO J = 1, N-1
    IF(.NOT. Check_Root(A, J, N)) THEN
      PRINT *, '方程组无解或无唯一解'
      STOP '程序停止于函数Check_Root'
    ENDIF
    CALL Transfer_Max(A, B, J, N)
    DO I = J+1, N
      C = A(I, J)/A(J, J)
      B(I) = B(I) - B(J)*C
      A(I, J:N) = A(I, J:N) - A(J, J:N)*C
    ENDDO
  ENDDO
  ! 回代步
  X(N) = B(N)/A(N, N)
  DO I = N-1, 1, -1
    C = 0.0
    DO J = I+1, N
      C = C + A(I, J)*X(J)
    ENDDO
    X(I) = (B(I) - C)/A(I, I)
  ENDDO

END SUBROUTINE Gauss_Expunction

LOGICAL FUNCTION Check_Root(A, J, N)
  IMPLICIT NONE
! 变量定义
  INTEGER :: N
  REAL :: A(N,N)
  INTEGER :: J
  INTEGER :: I
 数据处理
  DO I = J, N
    IF(A(I, J) /= 0.0) THEN
      Check_Root = .TRUE.
      RETURN
    ENDIF
  ENDDO
  Check_Root = .FALSE.
  RETURN

END FUNCTION Check_Root

SUBROUTINE Transfer_Max(A, B, J, N)
  IMPLICIT NONE
! 变量定义
  INTEGER :: N
  REAL :: A(N,N)
  REAL :: B(N)
  INTEGER :: J
  INTEGER :: I, K 
! 数据处理
  K = J
  DO I = J, N
    IF(ABS(A(I, J)) > ABS(A(J, J))) K = I
  ENDDO
  IF(K == J) RETURN
  CALL Swap_Element(B(J), B(K))
  DO I = J, N
    CALL Swap_Element(A(J, I), A(K, I))
  ENDDO

CONTAINS
! 元素扫描子程序
  SUBROUTINE Swap_Element(A, B)
    IMPLICIT NONE

    REAL :: A, B
    REAL :: C

    C = A
    A = B
    B = C

  END SUBROUTINE Swap_Element

END SUBROUTINE Transfer_Max

! 高斯消去法的范例
PROGRAM TEST1402
  IMPLICIT NONE
! 变量定义
  REAL, ALLOCATABLE :: A(:,:), B(:), X(:)
  INTEGER :: N, I, J, Err
! 屏幕提示并分配存储空间
  WRITE(*, '(1X,A)')'请输入方程组中的方程个数:'
  READ(*, *)N
  ALLOCATE(A(N,N), STAT=Err)
  IF(0 /= Err) STOP '动态分配方程组系数矩阵失败!'
  ALLOCATE(B(N), STAT=Err)
  IF(0 /= Err) STOP '动态分配方程组结果数组失败!'
  ALLOCATE(X(N), STAT=Err)
  IF(0 /= Err) STOP '动态分配方程组变量数组失败!'
! 屏幕提示
  WRITE(*, '(1X,A,I3,A,I3,A)')'请输入',N,'X',N,'个方程组系数:'
  DO I = 1, N
    READ(*, *)(A(I, J), J = 1, N)
  ENDDO
  WRITE(*, '(1X,A,I3,A,I3,A)')'请输入',N,'个方程组的值:'
  READ(*, *)(B(I), I = 1, N)
! 高斯消去法求解
  CALL Gauss_Expunction(A, B, X, N)

  WRITE(*, '(1X,A)')'方程组的解如下:'
  WRITE(*, '(1X,"X(",I2,") = ", F10.6)')(I, X(I), I=1,N)

END PROGRAM TEST1402