program trisolve(a,b,c,r,I,L,u,n)
c Solves the tridiagonal matrix system given below. a,b,c,r
c are input vectors; u is returned. Note a(I) and c(L) are
c not referenced.
c I is starting index; L is ending index
c
c | b(I) c(I) 0 0 0 0 | | u(I) | | r(I) |
c |a(I+1) b(I+1) c(I+1) 0 0 0 | |u(I+1)| |r(I+1)|
c | 0 a(I+2) b(I+2) c(I+2) 0 0 | |u(I+2)| |r(I+2)|
c | .......... | | .. | = | .. |
c | a(L-1) b(L-1) c(L-1)| |u(L-1)| |r(L-1)|
c | a(L) b(L) | | u(L) | | r(L) |
c
Implicit double precision (a-h,o-z)
parameter (NMAX=200)
dimension a(n),b(n),c(n),r(n),u(n),gamma(NMAX)
if (b(I) .eq. 0.0d0) go to 999
bet=b(I)
u(I)=r(I)/bet
c DECOMPOSITION & FORWARD SUBSTITUTION
do 10 j=I+1,L
gamma(j)=c(j-1)/bet
bet=b(j)-a(j)*gamma(j)
if (bet .eq. 0.0d0) go to 999
u(j)=(r(j)-a(j)*u(j-1))/bet
10 continue
c BACKSUBSTITUTION
do 20 j=L-1,1,-1
u(j)=u(j)-gamma(j+1)*u(j+1)
20 continue
go to 1000
999 write(3,666)
666 format('','ALGORITHM FAILS DIVIDE BY 0')
1000 return
end