C======================================================================
C NPAC $Id: tridiag_psi.f90,v 1.4 1995/11/21 20:42:41 hwyau Exp $
C======================================================================
C
      SUBROUTINE tridiag_psi(alpha,beta,rhs,u,temp,N)
C _____________________________________________________________________
C|                                                                     |
C| Subroutine, to solve a tridiagonal matrix, for the psi Poisson      |
C| equation.                                                           |
C|                                                                     |
C| The rhs(N,N) is the RHS, u(N,N) the solution, and beta(N,N) and     |
C| temp(N,N) workspaces arrays.                                        |
C| The tridiagonal matrix in this problem has the simple form of:      |
C|   diagonal elements = alpha +2.0                                    |
C|   off-diagonal      =  -1.0                                         |
C|                                                                     |
C| H W Yau.  17th of April, 1995.                                      |
C| Northeast Parallel Architectures Center.                            |
C| Syracuse University.                                                |
C|_____________________________________________________________________|
C|                                                                     |
C| Edit record:                                                        |
C| 17/Apr/1995: First edit, code modified from Numerical Recipes.      |
C| 19/Apr/1995: Customised for solving psi Poisson equation.           |
C| 15/Nov/1995: Rewrite, to avoid INDEPENDENT.                         |
C|              Need to add another 2-D array, beta().                 |
C|              Added HPF directives.                                  |
C| 17/Nov/1995: Option to not update boundary elements.                |
C| 21/Nov/1995: Proper handling of boundary elements.                  |
C|_____________________________________________________________________|
C
      IMPLICIT NONE
C
C .. Arguments ..
      INTEGER,INTENT(IN)              :: N
      REAL,INTENT(IN)                 :: alpha
      REAL,DIMENSION(:,:),INTENT(IN)  :: rhs
      REAL,DIMENSION(:,:),INTENT(OUT) :: u,temp,beta
C
C .. HPF directives ..
!HPF$ DISTRIBUTE *(*,BLOCK) :: beta,rhs,u,temp
C
C .. Local variables ..
      INTEGER i,k,j1,jN
      REAL a,b,c
C______________________________________________________________________
C
C Begin Executable
C______________________________________________________________________
C
      j1 = 2                        ! Start element to update
      jN = N-1                      ! End index to update
C
C Define matrix elements.
      a = -1.0                      ! Lower diagonal elements.
      b = alpha + 2.0               ! Diagonal elements.
      c = -1.0                      ! Upper diagonal elements.
C
C Do LU decomposition and forward substitution.
      u(j1,j1:jN)    = ( rhs(j1,j1:jN) - a*u(1,j1:jN) )/b
      beta(j1,j1:jN) = b
      DO k=j1+1,jN
        temp(k,j1:jN) = c/beta(k-1,j1:jN)
        beta(k,j1:jN) = b - a*temp(k,j1:jN)
        u(k,j1:jN)    = (rhs(k,j1:jN) - a*u(k-1,j1:jN))/beta(k,j1:jN)
      END DO
      temp(N,j1:jN) = c/beta(jN,j1:jN)
C
C Do backwards substitution.
      DO k=jN,j1,-1
        u(k,j1:jN) = u(k,j1:jN) - temp(k+1,j1:jN)*u(k+1,j1:jN)
      END DO
C
      RETURN
C
      END SUBROUTINE tridiag_psi
C
C======================================================================

