高斯消元法Fortran90程序
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
本文末给出Gauss-Jordan消去法的Fortran90源程序。
!/************************************************************* !程序：Gauss_Jordan消去法
!过程：Gauss_Jordan(aa,b,n,sgn)
!作用：aa为方阵，b为aa的逆，n为aa的阶
! sgn为标识符，1表示求逆成功，0表示求逆失败
!调用格式为：call Gauss_Jordan(aa,b,n,sgn)
!*************************************************************/ subroutine Gauss_Jordan(aa,b,n,sgn)
implicit none
integer(4):: n,sgn
real(8):: aa(n,n),b(n,n)
integer(4):: i,j,k
real(8),allocatable:: a(:,:)
real(8):: t
allocate(a(n,n))
a=aa ! a代替aa进行运算
sgn=1
! 初始化b为单位阵
do i=1,n
do j=1,n
if(i==j) then
b(i,j)=1
else
b(i,j)=0
end if
end do
end do
! Gauss_Jordan消去法过程
do k=1,n
if(a(k,k)==0) then
sgn=0;EXIT
end if
! 化第k行使得a(k,k)为1
t=1.0d0/a(k,k)
do j=k,n
a(k,j)=a(k,j)*t
end do
do j=1,n
b(k,j)=b(k,j)*t
end do
! 完成第k列的计算
do i=1,n
if(i/=k)then
t=a(i,k)
do j=k,n
a(i,j)=a(i,j)-a(k,j)*t
end do
do j=1,n
b(i,j)=b(i,j)-b(k,j)*t
end do
end if
end do
end do
end subroutine Gauss_Jordan
program equation_solve
implicit none
integer n,sgn
parameter(n=2)
real t
integer i,j,k
real aa(n,n),b(n,n),c(n,1),x(n,1)
real,allocatable:: a(:,:)
read*,aa,c
allocate(a(n,n))
a=aa ! a代替aa进行运算
sgn=1
! 初始化b为单位阵
do i=1,n
do j=1,n
if(i==j) then
b(i,j)=1
else
b(i,j)=0
end if
end do
end do
! Gauss_Jordan消去法过程
do k=1,n
if(a(k,k)==0) then
sgn=0
EXIT
end if
! 化第k行使得a(k,k)为1
t=1.0/a(k,k)
do j=k,n
a(k,j)=a(k,j)*t
end do
do j=1,n
b(k,j)=b(k,j)*t
end do ! 完成第k列的计算
do i=1,n
if(i/=k)then
t=a(i,k)
do j=k,n
a(i,j)=a(i,j)-a(k,j)*t
end do
do j=1,n
b(i,j)=b(i,j)-b(k,j)*t
end do
end if
end do
end do
x=matmul(b,c)
write(*,20) ((b(i,j),j=1,n),i=1,n)
20 format(5X,2F9.3)
write(*,10) (x(i,1),i=1,n)
10 format(5X,F5.3)
end
!正确的编程，但是需要能求得出逆矩阵的矩阵方程
! write(*,20) ((b(i,j),j=1,n),i=1,n)
! 20 format(5X,2F9.3) !print*," aa的逆矩阵 ",b !x=matmul(b,c)
!write(*,10) (x(i,1),i=1,n)
!10 format(5X,F5.3) !print*," aa的逆矩阵 ",b !end