亚当高斯消元法解线性方程组

来源：互联网发布：淘宝摄影棚制作编辑：程序博客网时间：2024/05/21 10:19

PROGRAM hw
IMPLICIT NONE
INTEGER,PARAMETER::DBL=SELECTED_REAL_KIND(p=13)
REAL(KIND=DBL),ALLOCATABLE,DIMENSION(:,:)::a
REAL(KIND=DBL),ALLOCATABLE,DIMENSION(:)::b
REAL(KIND=DBL),ALLOCATABLE,DIMENSION(:)::soln

INTEGER::IO_state
CHARACTER(LEN=128)::IO_msg

INTEGER::error_flag
INTEGER::i,j,n
CHARACTER(LEN=128)::file_name

    ! FILE valid?
    WRITE(*,*) 'Enter the file name'
    DO
        READ(*,*) file_name
        OPEN(UNIT=1,FILE=file_name,IOSTAT=IO_state,IOMSG=IO_msg,STATUS='OLD',ACTION='READ')
        IF(IO_state==0) EXIT
        WRITE(*,*) IO_msg
        WRITE(*,*) 'Re-enter the file name:'
    END DO

    READ(1,*) n
    ALLOCATE(a(n,n),b(n),soln(n))
    DO i=1,n
        READ(1,*,IOSTAT=IO_state)(a(i,j),j=1,n), b(i)
        data_valid:IF(IO_state/=0) THEN
            WRITE(*,*) 'Invalid data, check input file'
            EXIT
        END IF data_valid

    END DO
    WRITE(*,'(/,1X,"Coefficients:")')
    DO i=1,n
        WRITE(*,'(1X,7F11.4)') (a(i,j),j=1,n),b(i)
    END DO
    CALL dsimul(a,b,soln,n,n,error_flag)
    IF(error_flag/=0) THEN
        WRITE(*,*) 'Zero pivot encountered'
    ELSE
        DO i=1,n
            WRITE(*,'(1X, I3,2X,G15.6)') i, soln(i)
        END DO
    END IF
    DEALLOCATE(a,b,soln)
END PROGRAM hw

SUBROUTINE dsimul(a,b,soln,ndim,n,error)
IMPLICIT NONE
INTEGER,PARAMETER::DBL=SELECTED_REAL_KIND(p=13)
REAL(KIND=DBL),PARAMETER::EPSILON=1.0E-12

INTEGER,INTENT(IN)::ndim
REAL(KIND=DBL),INTENT(IN),DIMENSION(ndim,ndim)::a
REAL(KIND=DBL),INTENT(IN),DIMENSION(ndim)::b
REAL(KIND=DBL),INTENT(OUT),DIMENSION(ndim)::soln
INTEGER,INTENT(IN)::n
INTEGER,INTENT(OUT)::error

REAL(KIND=DBL),DIMENSION(n,n)::a1
REAL(KIND=DBL)::factor
INTEGER::irow,ipeak,jrow
REAL(KIND=DBL)::temp
REAL(KIND=DBL),DIMENSION(n)::temp1
    a1=a(1:n,1:n)
    soln=b(1:n)
    mainloop:DO irow=1,n
        ipeak=irow
        max_pivot:DO jrow=irow+1,n
            IF(ABS(a1(jrow,irow))>ABS(a1(ipeak,irow))) THEN
                ipeak=jrow
            END IF
        END DO max_pivot
        ! check for singular
        singular:IF(ABS(a1(ipeak,irow))<EPSILON) THEN
            error=1
            RETURN
        END IF singular
        ! SWAP
        swap_eqn:IF(ipeak/=irow) THEN
            temp1=a1(irow,1:n)
            a1(irow,1:n)=a1(ipeak,1:n)
            a1(ipeak,1:n)=temp1
            temp=soln(irow)
            soln(irow)=soln(ipeak)
            soln(ipeak)=temp
        END IF swap_eqn
        ! eliminate
        eliminate:DO jrow=1,n
            IF(jrow/=irow) THEN
                factor=-a1(jrow,irow)/a1(irow,irow)
                a1(jrow,1:n)=a1(irow,1:n)*factor+a1(jrow,1:n)
                soln(jrow)=soln(irow)*factor+soln(jrow)
            END IF
        END DO eliminate
    END DO mainloop
    ! divide
    divide:DO irow=1,n
    soln(irow)=soln(irow)/a1(irow,irow)
    END DO divide
    error=0
END SUBROUTINE dsimul