高斯消元法Fortran90程序

本文末给出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的逆矩阵 \

!x=matmul(b,c)

!write(*,10) (x(i,1),i=1,n)

!10 format(5X,F5.3) !print*,\" aa的逆矩阵 \

!end