subroutine tridia(ndim,n,a,b,c,d,k,l,errk,errl,erru, & u,ier,err,umax) c General information: c Subroutine tridia solves a tridiagonal system of equations for the c n+1 unknowns u(0), u(1), ..., u(n-1), u(n). This system must c consist of n+1 equations of the form: c c b(0) u(0) + c(0) u(1) = d(0) first equation c a(i) u(i-1) + b(i) u(i) + c(i) u(i+1) = d(i) for i = 1, 2, ..., n-1 c a(n) u(n-1) + b(n) u(n) = d(n) final equation c Subroutine tridia uses the Thomas algorithm. This method does not c use pivoting; it may crash or be inaccurate if the tridiagonal c system is not strictly diagonally dominant or symmetric positive c definite. c Always check the error code ier after using tridia; if ier is 2 or c more, there may be a problem. c Copyright 1996 Leon van Dommelen c Version 1.0 Leon van Dommelen 1/2/97 c Computational procedure: c Subroutine tridia works by reducing the given system to a c bidiagonal system of equations of the form c c u(0) + k(0) u(1) = l(0) first equation c u(i) + k(i) u(i+1) = l(i) for i = 1, 2, ..., n-1 c u(n) = l(n) final equation c c and then solving this system in inverse order. c Arguments: c Avoid typos: implicit none c Input: declared maximum equation number (0 <= ndim): integer ndim c Input: actual last equation number (0 <= n <= ndim): integer n c Input: the 'subdiagonal' coefficients a(i) in the equations: double precision a(0:ndim) c Input: the 'diagonal' coefficients b(i) in the equations: double precision b(0:ndim) c Input: the 'superdiagonal' coefficients c(i) in the equations: double precision c(0:ndim) c Input: the right hand sides d(i) of the equations: double precision d(0:ndim) c Work: the superdiagonal of the reduced system: double precision k(0:ndim) c Work: the right hand side vector of the reduced system: double precision l(0:ndim) c Work: upper bounds for the absolute errors in the k(i): double precision errk(0:ndim) c Work: upper bounds for the absolute errors in the l(i): double precision errl(0:ndim) c Work: upper bounds for the absolute errors in the u(i): double precision erru(0:ndim) c Output: the solution values u(i): double precision u(0:ndim) c Output: an error code: integer ier c Code ier is set equal to the sum of: c 1: if there is only a single equation; c 2: if the equations seem effectively singular; c 4: if the subroutine crashed: no solution was found. c Always check that ier is less than 2 after using tridia. c Output: the estimated maximum error in the coefficients of u: double precision err c This error is based on the round-off error in the coefficients c of the equations. It ignores the possibility of underflow. c A small value of err implies that the solution will be accurate c as long as underflow plays no role. However, the value of err is c typically very pessimistic for reasons that are presently unclear. c Part of it seems to be significant error cancellation in the back c substitution stage. c Output: the maximum absolute coefficient of u: double precision umax c Local variables: c An index to keep track whether tridia has been initialized: integer init c Index over the equations: integer i c The diagonal coefficient of the equation after eliminating a(i): double precision bb c The relative error in bb: double precision relbb c The machine epsilon (the relative error in real numbers): double precision eps save eps c Numerical codes for the various problems that tridia may encounter: integer oneeq,singul,crash parameter (oneeq=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*,'*** tridia: ', & 'Invalid arguments issued to subroutine tridia.' print*,' Correct the program code.' stop endif 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 Initialize the error code: ier=0 if(n.eq.0)ier=oneeq c Reduce the first equation (i=0) by dividing it by b(0) c Check for exact singularity: if(b(0).eq.zero)then ier=ier+crash return endif c The normalized superdiagonal coefficient (only exists for n > 0): if(n.gt.0)k(0)=c(0)/b(0) c The normalized right hand side in array l: l(0)=d(0)/b(0) c Compute error bounds for these values: if(n.gt.0)errk(0)=abs(k(0))*(eps+eps+eps) errl(0)=abs(l(0))*(eps+eps+eps) c Reduce the remaining equations to the desired bidiagonal form: do 100 i=1,n c Substract a(i) times the reduced (i-1)-th equation; c u(i-1) + k(i-1) u(i) = l(i-1), c from the i-th equation; c a(i) u(i-1) + b(i) u(i) + c(i) u(i+1) = d(i), c to eliminate the subdiagonal term: c The new main-diagonal coefficient: bb=b(i)-a(i)*k(i-1) c Check for exact singularity: if(bb.eq.zero)then ier=ier+crash return endif c The new right hand side: l(i)=d(i)-a(i)*l(i-1) c Now divide the equation by the diagonal coefficient: c The reduced superdiagonal coefficient: if(i.lt.n)k(i)=c(i)/bb c The reduced right hand side: l(i)=l(i)/bb c Compute error bounds for these values relbb=(abs(b(i)/bb)+abs(a(i)*k(i-1)/bb))*eps & +abs(a(i)/bb)*errk(i-1)+eps if(i.lt.n)errk(i)=abs(k(i))*(eps+relbb+eps) errl(i)=(abs(d(i)/bb)+abs(a(i)*l(i-1)/bb))*eps & +abs(a(i)/bb)*errl(i-1)+abs(l(i))*(relbb+eps) 100 continue c Solve the final equation: c The equation is trivial: u(n)=l(n) c Compute the error bound: erru(n)=errl(n) err=erru(n) umax=abs(u(n)) c Solve the remaining equations: do 200 i=n-1,0,-1 c Solution of the reduced equation u(i)=l(i)-k(i)*u(i+1) c Compute the error bound: erru(i)=errl(i)+abs(u(i+1))*errk(i) & +abs(k(i))*erru(i+1)+abs(u(i))*eps err=max(err,erru(i)) umax=max(umax,abs(u(i))) 200 continue c Warn for singularity: if(err.gt.umax)ier=ier+singul c Exit: return end