数值分析—行主元消去法解线性方程组—FORTRAN程序

数值分析—行主元消去法解线性方程组—FORTRAN程序

program main

implicit none

real8,dimension( :,: ),allocatable::A

real8,dimension( : ),allocatable::B, x

integer::i,n

print*,‘请输入线性方程组的元数n:’

read*,n

!为动态数组分配内存

allocate(A(n,n),B(n),x(n))

open(111,file=‘A.txt’)

open(222,file=‘B.txt’)

open(333,file=‘ans.txt’)

do i=1,n

read(111,)A(i,1:n)!按行输入系数矩阵A

end do

print,“A=”

do i=1,n

print’(1X,3F6.1)’,A(i,1:N)

end do

do i=1,n

read(222,)B(i)

end do

print,“B=”

do i=1,n

print"(1X,3F6.1)",B(i)

end do

call Gauss_hzy(n,A,B,x)

print*,‘x=’,x

write(333,*)‘x=’,x

end program main

subroutine Gauss_hzy(n,A,B,x)

implicit none

integer::n,i,j,k,m,jmax

real8::amaxh,s

real8,dimension(n,n)::A,C

real8,dimension(n)::B,x

integer,dimension(n)::l

do k=1,n

l(k)=k

enddo

do k=1,n-1

!选最大元素所在列

amaxh=abs(A(k,k))

jmax=k

do j=k+1,n

if(abs(A(k,j))>amaxh)then

amaxh=abs(A(k,j))

jmax=j

end if

end do

!判断矩阵A是否奇异

if(abs(amaxh)<1e-17)then

print,‘ERROR!’

stop

elseif(jmax/=k)then

!交换列

do i=1,n

l(k)=jmax

call swap(A(i,k),A(i,jmax))

end do

end if

!消元

do i=k+1,n

C(i,k)=A(i,k)/A(k,k)

B(i)=B(i)-C(i,k)*B(k)

do j=k+1,n

A(i,j)=A(i,j)-C(i,k)*A(k,j)

end do

end do

end do

!回带求解

x(n)=B(n)/A(n,n)

do i=n-1,1,-1

s=0

do j=n,i+1,-1

s=s+A(i,j)*x(j)

end do

x(i)=(B(i)-s)/A(i,i)

end do

!恢复解的正常次序

do k=n,1,-1

if(l(k)/=k)then

m=l(k)

call swap(x(k),x(m))

end if

end do

end subroutine Gauss_hzy

subroutine swap(a,b)

implicit none

real*8::a,b,c

c=a

a=b

b=c

end subroutine swap

!输入：n=3

!输入A= 0 2 1

! 3 5 -5

! 2 4 -2

!输入B=-2 1 2

!输出：x= 7.00000191 -2.00000048 2.00000072

!输入：n=3

!输入A= 1 1 1

! 1 -1 1

! 1 1 -1

!输入B= 1 3 3

!输出x= 3.0000000000000000 -1.0000000000000000 -1.0000000000000000