SUBROUTINE TRIDIA(N,DIAGS,B,X)

c     Subroutine to solve a tridiagonal system of N equations of the form
c        DIAGS(1,I)*X(I-1) + DIAGS(2,I)*X(I) + DIAGS(3,I)*X(I+1) = B(I)
c     for i = 1, 2, 3, ... n, to give the unknowns x(1), x(2), ..., x(n).
c     The coefficients DIAGS(1,1) and DIAGS(3,N) of the nonexisting
c     unknowns can be set to zero.

c     Note that the equations and arrays start from 1.
c     However, because of the peculiarities of Fortran array storage and
c     argument passing, the subroutine will work fine if DIAGS, B, and X
c     start from I=0 instead of I=1, as long as the correct number of
c     equations is specified (i.e., N+1 instead of N).

c     In the subroutine, the matrix DIAGS is altered: therefore, to
c     solve the same system but with different right hand side B again,
c     either restore the matrix DIAGS, or call TRIDI2 below instead.

c     No pivoting is performed, which is fine for diagonally dominant or
c     symmetric positive definite matrices.

      IMPLICIT NONE

c Input:

c     number of equations
      INTEGER N
c     the coefficients of the unknowns
      DOUBLE PRECISION DIAGS(3,N)
c     the right hand sides of the equations
      DOUBLE PRECISION B(N)

c Output:

c     the values of the unknowns
      DOUBLE PRECISION X(N)

c Local:

c     equation index
      INTEGER I
c     amount of the previous equation to substract to eliminate DIAGS(1,I)
      DOUBLE PRECISION FAC

c Executable:

c     Reduce the matrix to bidiagonal form:
      DO 10 I=2,N
c        Multiple of equation i-1 to substract to eliminate diags(1,i)
         FAC=DIAGS(1,I)/DIAGS(2,I-1)
c        Substract FAC times equation i-1 from equation i
         DIAGS(2,I)=DIAGS(2,I)-FAC*DIAGS(3,I-1)
         B(I)=B(I)-FAC*B(I-1)
c        Use the now free storage location DIAGS(1,I) to store FAC:
         DIAGS(1,I)=FAC
 10   CONTINUE

c     Solve the bidiagonal system:
      X(N)=B(N)/DIAGS(2,N)
      DO 20 I=N-1,1,-1
         X(I)=(B(I)-DIAGS(3,I)*X(I+1))/DIAGS(2,I)
 20   CONTINUE

      RETURN
      END



      SUBROUTINE TRIDI2(N,DIAGS,B,X)

c     Version of TRIDIA which can be used when TRIDIA has changed matrix
c     DIAGS.  To solve the SAME matrix with different right hand side
c     vectors, the most efficient procedure is to call TRIDIA to solve
c     the first system, and TRIDI2 for the remaining right hand side
c     vectors.

      IMPLICIT NONE

c Input:

c     number of equations
      INTEGER N
c     the coefficients of the unknowns
      DOUBLE PRECISION DIAGS(3,N)
c     the right hand sides of the equations
      DOUBLE PRECISION B(N)

c Output:

c     the values of the unknowns
      DOUBLE PRECISION X(N)

c Local:

c     equation index
      INTEGER I
c     amount of the previous equation to substract to eliminate DIAGS(1,I)
      DOUBLE PRECISION FAC

c Executable:

c     Transform B:
      DO 10 I=2,N
c        as computed by TRIDIA
         FAC=DIAGS(1,I)
c        transform
         B(I)=B(I)-FAC*B(I-1)
 10   CONTINUE

c     Solve the bidiagonal system:
      X(N)=B(N)/DIAGS(2,N)

      DO 20 I=N-1,1,-1
         X(I)=(B(I)-DIAGS(3,I)*X(I+1))/DIAGS(2,I)
 20   CONTINUE

      RETURN
      END