CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Looking for Periodic Tridiagonal Solver

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Display Modes
Old   August 26, 2004, 17:20
Default Looking for Periodic Tridiagonal Solver
  #1
Wen Long
Guest
 
Posts: n/a
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
  Reply With Quote

Old   August 26, 2004, 19:46
Default Re: Looking for Periodic Tridiagonal Solver
  #2
Tom
Guest
 
Posts: n/a
Hi,

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

Tom
  Reply With Quote

Old   August 26, 2004, 22:29
Default Re: Looking for Periodic Tridiagonal Solver
  #3
Wen Long
Guest
 
Posts: n/a
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
  Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
OpenCL linear solver for OpenFoam 1.7 (alpha) will come out very soon qinmaple OpenFOAM Announcements from Other Sources 4 August 10, 2012 11:00
Question about the MIGAL solver bearcat Phoenics 0 February 4, 2010 19:02
block tridiagonal system solver for periodic bc Abdullah Main CFD Forum 8 February 23, 2007 14:15
Symmetry plane error in solver Santiago Orrego. CFX 6 January 31, 2007 08:09
Block Tridiagonal Solver Abdulhafid M. Elfaghi Main CFD Forum 2 December 23, 2006 13:20


All times are GMT -4. The time now is 14:44.