      program elim

c Solve the linear system
c
c   eps*x1 + x2 = 1,  x1 + x2 = 2
c
c via Gaussian elimination.  This is to contrast the
c difference between solving the two equations in the given
c ordering and interchanging the equations before solving.
c For small values of e (say 1.E-5 > eps > 1.E-14) the computation
c of x1 loses a lot of precision if the rows are not interchanged.

      implicit real*4 (a-z)
      double precision eps,exact1,exact2
c read epsilon
      read*,eps
c initialize matrix and right hand side
      b1=1.0
      b2=2.0
      a11=eps
      a12=1.0
      a21=1.0
      a22=1.0
c print header
      print1
c print out exact solution
      exact1=1.d0/(1.d0-eps)
      exact2=(2.d0*eps-1.d0)/(1.d0*eps-1.d0)
      print2,'exact',exact1,exact2
c do Gaussian elimination
      m=a21/a11
      a22=a22-m*a12
      b2=b2-m*b1
      x2=b2/a22
      x1=(b1-a12*x2)/a11
c print out answer and error      
      err1=(exact1-x1)/exact1
      err2=(exact2-x2)/exact2
      print2,'no interchange',x1,x2,err1,err2
c reinitialize matrix and right hand side with rows interchanged
      b1=2.0
      b2=1.0
      a11=1.0
      a12=1.0
      a21=eps
      a22=1.0
c do Gaussian elimination
      m=a21/a11
      a22=a22-m*a12
      b2=b2-m*b1
      x2=b2/a22
      x1=(b1-a12*x2)/a11
c print out answer and error      
      err1=(exact1-x1)/exact1
      err2=(exact2-x2)/exact2
      print2,'interchange',x1,x2,err1,err2
c format statements
1     format(9x,'METHOD         X1        X2    RelErr1   RelErr2')
2     format(a15,3x,2f10.6,1p2e10.2)
      end
