c     Copyright (C) 2008-2014 Vincent Eymet
c
c     KDISTRIBUTION is free software; you can redistribute it and/or modify
c     it under the terms of the GNU General Public License as published by
c     the Free Software Foundation; either version 3, or (at your option)
c     any later version.

c     KDISTRIBUTION is distributed in the hope that it will be useful,
c     but WITHOUT ANY WARRANTY; without even the implied warranty of
c     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
c     GNU General Public License for more details.

c     You should have received a copy of the GNU General Public License
c     along with KDISTRIBUTION; if not, see <http://www.gnu.org/licenses/>
c
      subroutine tridiag(n,alpha,beta,gamma,y,x,valid)
      implicit none
      include 'max.inc'
c
c     Purpose: to solve a tridiagonal linear system
c
c     Inputs:
c       + n: dimension of the system
c       + alpha: array of below-diagonal terms , with alpha(1)=0
c       + beta: array of diagonal terms
c       + gamma: array of above-diagonal terms, with gamma(n)=0
c       + y: array of y values
c
c     Ouputs:
c       + x: array of x values
c       + valid: 0 if solution is not valid; 1 if solution is valid
c

c     input
      integer n
      double precision alpha(1:kmax)
      double precision beta(1:kmax)
      double precision gamma(1:kmax)
      double precision y(1:kmax)
c     output
      double precision x(1:kmax)
      integer valid
c     temp
      integer i
      double precision delta(1:kmax)
      double precision yp(1:kmax)
      double precision yv(1:kmax)
      double precision epsilon
      parameter(epsilon=1.0D-6)
      logical is_nan
c     label
      integer strlen
      character*(Nchar_mx) label
      label='subroutine tridiag'

      do i=1,n
         call test_nan(alpha(i),is_nan)
         if (is_nan) then
            call error(label)
            write(*,*) 'alpha(',i,')=',alpha(i)
            stop
         endif
         call test_nan(beta(i),is_nan)
         if (is_nan) then
            call error(label)
            write(*,*) 'beta(',i,')=',beta(i)
            stop
         endif
         call test_nan(gamma(i),is_nan)
         if (is_nan) then
            call error(label)
            write(*,*) 'gamma(',i,')=',gamma(i)
            stop
         endif
         call test_nan(y(i),is_nan)
         if (is_nan) then
            call error(label)
            write(*,*) 'y(',i,')=',y(i)
            stop
         endif
      enddo ! i

c     Debug
c      write(*,*) 'This is:',label(1:strlen(label))
c      write(*,*) 'y(1)=',y(1)
c      call test_nan(y(1),is_nan)
c      write(*,*) 'is_nan=',is_nan
c      do i=1,n
c         write(*,*) 'alpha(',i,')=',alpha(i)
c      enddo
c      do i=1,n
c         write(*,*) 'beta(',i,')=',beta(i)
c      enddo
c      do i=1,n
c         write(*,*) 'gamma(',i,')=',gamma(i)
c      enddo
c      do i=1,n
c         write(*,*) 'y(',i,')=',y(i)
c      enddo
c     Debug

c     tests
      if (n.lt.2) then
         call error(label)
         write(*,*) 'n=',n
         stop
      endif
      if (alpha(1).ne.0.0D+0) then
         call error(label)
         write(*,*) 'alpha(1)=',alpha(1)
         stop
      endif
      if (gamma(n).ne.0.0D+0) then
         call error(label)
         write(*,*) 'gamma(n)=',gamma(n)
         stop
      endif
c     tests

      if (beta(1).eq.0.0D+0) then
         call error(label)
         write(*,*) 'beta(1)=',beta(1)
         stop
      else
         delta(1)=gamma(1)/beta(1)
         yp(1)=y(1)/beta(1)
      endif
      do i=2,n
         if (beta(i)-delta(i-1)*alpha(i).eq.0.0D+0) then
         call error(label)
            write(*,*) 'beta(',i,')=',beta(i)
            write(*,*) 'delta(',i-1,')*alpha(',i,')=',
     &           delta(i-1)*alpha(i)
            stop
         else
            delta(i)=gamma(i)/(beta(i)-delta(i-1)*alpha(i))
            yp(i)=(y(i)-yp(i-1)*alpha(i))/(beta(i)-delta(i-1)*alpha(i))
         endif
      enddo ! i
c     Debug
c      do i=1,n
c         write(*,*) 'delta(',i,')=',delta(i),' yp(',i,')=',yp(i)
c      enddo
c     Debug

      x(n)=yp(n)
      do i=n-1,1,-1
         x(i)=yp(i)-x(i+1)*delta(i)
      enddo
c     Debug
c      do i=1,n
c         write(*,*) 'x(',i,')=',x(i)
c      enddo
c     Debug

c     Tests
      yv(1)=beta(1)*x(1)+gamma(1)*x(2)
      do i=2,n-1
         yv(i)=alpha(i)*x(i-1)+beta(i)*x(i)+gamma(i)*x(i+1)
      enddo
      yv(n)=alpha(n)*x(n-1)+beta(n)*x(n)
      valid=1
      do i=1,n
         if (dabs(yv(i)-y(i)).gt.epsilon) then
            valid=0
            goto 123
         endif
      enddo
c     Tests
 123  continue
c     Debug
c      write(*,*) 'valid=',valid
c      do i=1,n
c         write(*,*) 'y(',i,')=',y(i),' yv(',i,')=',yv(i)
c      enddo
c     Debug
      return
      end