程序编写如下：
real,allocatable::a(:,:),b(:),c(:) print*,'输入未知数个数n' read*,n
allocate(a(n,n)) allocate(b(n)) allocate(c(n))
print*,'输入系数矩阵a' call input(a,n)
print*,'输入等值矩阵b' read*,b
print*,'联立方程组' call output(a,b,n)
call Gauss_jordan(a,b,c,n) print*,\求解\do i=1,n
print 10,i,c(i) enddo
10 format('x',i1,'=',f8.4) deallocate(a) deallocate(b) deallocate(c)
end
subroutine input(a,n) real a(n,n) do i=1,n
read*,(a(i,j),j=1,n) enddo end
subroutine Gauss_jordan(a,b,c,n) dimension a(n,n),b(n),c(n) call up(a,b,n) call low(a,b,n) forall(i=1:n) c(i)=b(i)/a(i,i) endforall end
subroutine output(a,b,n) real a(n,n),b(n) do i=1,n
print 10,a(1,1),i do j=2,n
if(a(i,j)>0)then print 20,a(i,j),j
else
print 30,abs(a(i,j)),j endif enddo
print 40,b(i) enddo
10 format(f5.2,'x',i1\\) 20 format('+',f5.2,'x',i1\\) 30 format('-',f5.2,'x',i1\\) 40 format('=',f8.4) end
subroutine up(a,b,n) real a(n,n),b(n) do i=1,n-1 do j=i+1,n
p=a(j,i)/a(i,i)
a(j,i:n)=a(j,i:n)-a(i,i:n)*p b(j)=b(j)-b(i)*p enddo enddo end
subroutine low(a,b,n)
real a(n,n),b(n) do i=n,2,-1
do j=i-1,1,-1
p=a(j,i)/a(i,i)
a(j,1:i)=a(j,1:i)-a(i,1:i)*p b(j)=b(j)-b(i)*p enddo enddo end
方程组运行结果