数值分析—全主元消去法解线性方程组—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_qzy(n,A,B,x)
print*,‘x=’,x
write(333,*)‘x=’,x
end program main

subroutine Gauss_qzy(n,A,B,x)
implicit none
integer::n,i,j,k,m,imax,jmax
real8::amaxh,amaxl,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
!选最大元素所行
amaxl=abs(A(k,k))
imax=k
do i=k+1,n
if(abs(A(i,k))>amaxl) then
amaxl=abs(A(i,k))
imax=i
endif
enddo
!选最大元素所在列
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(abs(amaxl)<1e-17)then
print*,‘ERROR!’
stop
elseif(imax/=k)then
!交换行
do j=k,n
call swap(A(k,j),A(imax,j))
end do
call swap(B(k),B(imax))
end if
if(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_qzy

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.0000000000000053 -2.0000000000000009 2.0000000000000022