!==================================
subroutine GAUSS1(n,A,b,ic)

real :: A(4,4),b(4),tmp
integer :: i,j,n,k,ic,itmp,imax

eps=1.0E-10
ic=0

do k=1,n

 itmp=k
 do i=K+1,n
   if (abs(A(i,k))>abs(A(itmp,k))) then
      itmp=i
   end if
 end do
 
 imax=itmp
 
 do j=k,n
   tmp=A(imax,j)
   A(imax,j)=A(k,j)
   A(k,j)=tmp
 end do
 tmp=b(imax)
 b(imax)=b(k)
 b(k)=tmp
 
 if (abs(A(k,k))<eps) then
   ic=1
   return
 end if
      
 do j=k+1,n
   A(k,j)=A(k,j)/A(k,k)
 end do
 b(k)=b(k)/A(k,k)
 A(k,k)=1.0
 
 do i=K+1,n
   do j=k+1,n
     A(i,j)=A(i,j)-A(i,k)*A(k,j)
   end do
   b(i)=b(i)-A(i,k)*b(k)
   A(i,k)=0.0
 end do
 
end do

return
end subroutine GAUSS1
!==================================
subroutine GAUSS2(n,A,b,x)

real :: A(4,4),b(4),x(4),SUM,tmp
integer :: m,n,k

do k=n,1,-1
  
  tmp=0.0
  do m=k+1,n
    tmp=tmp+A(k,m)*x(m)
  end do
  SUM=tmp
  
  x(k)=b(k)-SUM
  
end do
return
end subroutine GAUSS2
!==================================