program agaus0
dimension a(4,4),b(4),x(4),js(4)
double precision a,b,x
data a/0.2368,0.1968,0.1582,1.1161,0.2471,0.2071,1.1675,0.1254,0.2568,1.2168,0.1768,0.1397,1.2671,0.2271,0.1871,0.1490/
data b/1.8471,1.7471,1.6471,1.5471/
n=4
call agaus(a,b,n,x,l,js)
if(l.ne.0)then
  write(*,10)(i,x(i),i=1,4)
end if
10format(1x,'x(',i2,')=',d15.6)
end


 subroutine agaus(a,b,n,x,l,js)
!a--双精度实型二维数组，体积为n*n，输入参数。存放方程组的系数矩阵，返回时将被破坏。
!b--双精度实型一维数组，长度为n，输入参数。存放方程组右端向量，返回时将被破坏。
!n--整型变量，输入参数。存放方程组的阶数。
!x--双精度实型一维数组，长度为n，输出参数，返回方程组的解向量。
!l--整型变量，输出参数。若返回l=0，说明方程组的系数矩阵奇异，求解失败；若l/=0，说明正常返回。
!Js--整型一维数组，长度为n，为本子程序的工作数组。
dimension a(n,n),x(n),b(n),js(n)
double precision a,b,x,t
l=1
do 50 k=1,n-1
do 210  i=k,n
do 210  j=k,n
  if(abs(a(I,j)).gt.d)then
  d=abs(a(I,j))
js(k)=j
is=i
end if
210 continue
If(d+1.0.eq.1.0)then
l=0
else
  if(js(k).ne.k)then
   do 220 i=1,n
     t=a(i,k)
     a(i,k)=a(i,js(k))
     a(i,js(k))=t
220  continue
   end if
if(is.ne.k)then
 do 230 j=k,n
   t=a(k,j)
a(k,j)=a(is,j)
a(is,j)=t
230  continue
t=b(k)
b(k)=b(is)
b(is)=t
end if
end if
if(l.eq.0)then
write(*,100)
return
end if
do 10 j=k+1,n
a(k,j)=a(k,j)/a(k,k)
10  continue
b(k)=b(k)/a(k,k)
do 30 i=k+1,n
  do 20 j=k+1,n
    a(i,j)=a(i,j)-a(i,k)*a(k,j)
 20  continue
    b(i)=b(i)-a(i,k)*b(k)
30  continue
50  continue
   If(abs(a(n,n))+1.0.eq.1.0)then
     l=0
      write(*,100)
      return
    end if
  x(n)=b(n)/a(n,n)
do 70  i=n-1,1,-1
t=0.0
do 60  j=i+1,n
t=t+a(i,j)*x(j)
60  continue
x(i)=b(i)-t
70 continue
100  format(1x,'fail')
js(n)=n
do  150 k=n,1,-1
  if(js(k).ne.k)then
    t=x(k)
    x(k)=x(js(k))
    x(js(k))=t
  end if
150 continue
return

end