Hi, folks:

Can somebody post a code for Periodic Tridiagonal Solver?

Simply a tridiagonal equation Ax=b with tridiagonal elements and corner elements A(1,n) and A(n,1) non zero due to periodic boundary condition.

Thanks,

Wen

Hi,

You should have a look at www.netlib.org. They have probably a solver.

Tom

Thanks, I think I figured out. Here is my subroutine:

Subroutine cyclic (alpha, beta,gama, b, x, N) !solve the following cyclic tridiagonal equation

!************************************************* ************!

!* *!

!* (B1 C1 A1 ) (x1) (b1) *!

!* (A2 B2 C2 ) (x2) (b2) *!

!* ( A3 B3 C3 ) (x3) (b3) *!

!* ( A4 B4 C4 ) (x4)==== *!

!* ( A5 B5 C5 ) ... ==== ... *!

!* ( ... ... ... ) ... ... *!

!* ( AN-1 BN-1 CN-1) (xN-1) (bN-1)*!

!* (CN AN BN ) (xN) (bN) *!

!* *!

!* *!

!************************************************* ************! ! where A are alpha, B are beta, C are gama c----dummy variables--

Integer N

double precision alpha(N), beta(N),gama(N),b(N),x(N) c----local variables--

double precision betaPrime(N),u(N),z(N),c,fact

c=-beta(1) !the minus sign makes sure betaPrime(1) non zero

betaPrime(1) = beta(1) - c !the first tirdiagonal element

do i=2,N-1,1

betaPrime(i)=beta(i)

enddo

betaPrime(N) = beta(N) - alpha(1) * gama(N) / c !the last tridiagonal element

call tri_ge(alpha, betaPrime, gama, b, x, N) !solve for x

u(1) = c; !first u,

do i=2,N-1

u(i)=0.d0

enddo

u(N) = gama(N) !last U is same as alpha

call tri_ge(alpha, betaPrime, gama, u, z, N) !solve for z

fact = (x(1) + alpha(1) * x(N) / c)

& / (1.0 + z(1) + alpha(1) * z(N) / c)

do i=1,N

x(i) =x(i)- fact * z(i) !construct final results

enddo

Return

End

Subroutine tri_ge(alpha,beta,gama,b, x, N) c----basically same as subroutine trig()----but allows diagonal variables not equal to unit

!************************************************* ************!

!* *!

!* (B1 C1 ) (x1) (b1) *!

!* (A2 B2 C2 ) (x2) (b2) *!

!* ( A3 B3 C3 ) (x3) (b3) *!

!* ( A4 B4 C4 ) (x4)==== *!

!* ( A5 B5 C5 ) ... ==== ... *!

!* ( ... ... ... ) ... ... *!

!* ( An-1 Bn-1 Cn-1) (xn-1) (bn-1)*!

!* ( An Bn ) (xn) (bn) *!

!* *!

!* *!

!************************************************* ************! ! where A are alpha, B are beta, C are gama c----dummy variables---

integer N

double precision alpha(N),beta(N), gama(N), b(N), x(N) c-----local variables

double precision betaPrime(N), bPrime(N)

double precision coeff

integer i

!Perform forward elimination

betaPrime(1) = beta(1)

bPrime(1) = b(1)

do i=2,N

coeff = alpha(i) / betaPrime(i-1)

betaPrime(i) = beta(i) - coeff * gama(i-1)

bPrime(i) = b(i) - coeff * bPrime(i-1)

enddo

! Perform back substitution

x(N) = bPrime(N) / betaPrime(N)

do i=N-1,1,-1

x(i) = (bPrime(i) - gama(i) * x(i+1)) / betaPrime(i)

enddo

return

End