subroutine satis2(ndim,n,u0,r10,r20,u1,r11,r21,u2,r12,r22,
     & erru,errr,u,ier,err)

c General information:

c     Subroutine satis2 finds the combination "u" of three vectors "u0",
c     "u1", and "u2" so that two linear equations become satisfied.  Any
c     other linear equation satisfied by "u0", "u1", and "u2" remains
c     satisfied by "u".

c     Copyright 1996 Leon van Dommelen
c     Version 1.0 Leon van Dommelen 1/2/97

c Computational procedure:

c     The weighted average "u" can be written as:
c        u = u0 + f (u1-u0) + g (u2-u0)
c     with f and g to be determined.  The combined residuals in the
c     two equations are:
c        r1 = r10 + f (r11-r10) + g (r12-r10) = 0
c        r2 = r20 + f (r21-r20) + g (r22-r20) = 0
c     from which f and g can be determined.

c Arguments:

c     Avoid typos:
      implicit none

c     Input: declared maximum index of the vectors (1 <= ndim):
      integer ndim

c     Input: actual maximum index of the vectors (1 <= n <= ndim):
      integer n

c     Input: the first vector to use:
      double precision u0(0:ndim)

c     Input: the residuals in the two equations for vector u0:
      double precision r10,r20

c     Input: the second vector to use:
      double precision u1(0:ndim)

c     Input: the residuals in the two equations for vector u1:
      double precision r11,r21

c     Input: the third vector to use:
      double precision u2(0:ndim)

c     Input: the residuals in the two equations for vector u2:
      double precision r12,r22

c     Input: the typical absolute error in u:
      double precision erru
c        If zero, the machine round-off error is used to find erru.

c     Input: the typical absolute error in the residuals:
      double precision errr
c        If zero, the machine round-off error is used to find errr.

c     Output: the vector u producing zero residuals r1 and r2:
      double precision u(0:ndim)

c     Output: error code:
      integer ier
c        Always check the value of ier after satis2.  If ier is nonzero,
c     there may be a problem.

c     Output: estimated error in u based on the values of erru and errr:
      double precision err
c        A small value can give some confidence in the accuracy of the
c     obtained solution.

c Local variables:

c     An index to keep track whether satis2 has been initialized:
      integer init

c     The matrix of the system of equations for f and g:
      double precision a11,a12,a21,a22

c     Determinant of the matrix
      double precision det

c     The solution of the system:
      double precision f,g

c     Index over the coefficients of the vectors:
      integer i

c     The machine epsilon (the relative error in real numbers):
      double precision eps
      save eps

c     The error in the residuals
      double precision dr

c     The error in the u values
      double precision du

c     Estimated error in the determinant:
      double precision dd

c     Estimated errors in f and g:
      double precision df,dg

c     Maximum absolute coefficient in u0, u1, and u2:
      double precision umax

c     Maximum absolute coefficient in u0:
      double precision u0max

c     Maximum absolute coefficient in u1 - u0:
      double precision u10max

c     Maximum absolute coefficient in u2 - u0:
      double precision u20max

c     Numerical codes for the various problems that may occur:
      integer onedim,singul,crash
      parameter (onedim=1,singul=2,crash=4)

c     Constants defined to make changing precision easier:
      double precision zero,one,two
      parameter (zero=0.d0,one=1.d0,two=2.d0)

c Data statements:

      data init/0/

c Executable statements:

c     Check the arguments:
      if(ndim.lt.0 .or. n.lt.0 .or. n.gt.ndim)then
         print*,'*** satis2: invalid arguments'
         stop
      endif

c     Initialize the error code:
      ier=0
c     Cannot satisfy two equations using a scalar:
      if(n.lt.1)ier=onedim

c     Initialize the first time around:
      if(init.eq.0)then

c        Find the machine epsilon as the absolute error in the value 1:
         eps=two
 10      eps=eps/two
         if(one+eps.ne.one)goto 10

c        Mark as initialized:
         init=1
      endif

c     Solve the equations
c        r1 = r10 + f (r11-r10) + g (r12-r10) = 0
c        r2 = r20 + f (r21-r20) + g (r22-r20) = 0
c     for f and g:

c     Find the matrix of the system of equations:
      a11=r11-r10
      a12=r12-r10
      a21=r21-r20
      a22=r22-r20

c     Find the determinant of the matrix:
      det=a11*a22-a12*a21
      if(det.eq.zero)then
         ier=ier+crash
         return
      endif

c     Find f and g:
      f=(r20*a12-r10*a22)/det
      g=(r10*a21-r20*a11)/det

c     Combine the solutions:
      do 100 i=0,n
         u(i)=u0(i)+f*(u1(i)-u0(i))+g*(u2(i)-u0(i))
 100  continue

c     Estimate the error in the solution roughly:

c     Find the limiting u-values:
      umax=zero
      u0max=zero
      u10max=zero
      u20max=zero
      do 200 i=0,n
         umax=max(umax,abs(u0(i)),abs(u1(i)),abs(u2(i)))
         u0max=max(u0max,abs(u0(i)))
         u10max=max(u10max,abs(u1(i)-u0(i)))
         u20max=max(u20max,abs(u2(i)-u0(i)))
 200  continue

c     Error in the u-values:
      du=erru
      if(du.le.zero)du=eps*umax

c     Error in the residuals:
      dr=errr
      if(dr.le.zero)dr=eps*
     & max(abs(r10),abs(r11),abs(r12),abs(r20),abs(r21),abs(r22))

c     Error in the determinant:
      dd=two*dr*(abs(a11)+abs(a12)+abs(a21)+abs(a22))
      if(dd.ge.abs(det))
     & ier=ier-(ier/singul-ier/(2*singul)*2)*singul+singul

c     Error in f and g:
      df=abs(f)*dd/abs(det)
     & +dr*(abs(a12)+abs(a22)+two*(abs(r10)+abs(r20)))/abs(det)
      dg=abs(g)*dd/abs(det)
     & +dr*(abs(a21)+abs(a11)+two*(abs(r10)+abs(r20)))/abs(det)

c     Final error estimate:
      err=du+df*u10max+dg*u20max+two*du*(df+dg)
      if(err.gt.u0max+abs(f*u10max)+abs(g*u20max))
     & ier=ier-(ier/singul-ier/(2*singul)*2)*singul+singul

c     All done:
      return
      end