CFD Online URL
[Sponsors]
Home > Forums > Main CFD Forum

Subroutine for Block tridiagonal System

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

Reply
 
LinkBack Thread Tools Display Modes
Old   February 26, 2007, 22:39
Default Subroutine for Block tridiagonal System
  #1
Khan
Guest
 
Posts: n/a
I am looking for subroutine for solving block tridiagonal solver(periodic and non periodic). I found in the book of Anderson but some one in this forum suggest that it has an error, which i dont know.

Could any help me by sending the subroutine or help me in finding the error in the subroutine given in the book of anderson.

thanks

  Reply With Quote

Old   February 27, 2007, 02:49
Default Re: Subroutine for Block tridiagonal System
  #2
UT CAA guy
Guest
 
Posts: n/a
Here is the routine out of the Tannehill book "Computational Fluid Mechanics and Heat Transfer"

This will be an exact copy, and it looks like fortran77

c..subroutine to solve non-period block tridiaganol c..system of equations without pivoting strategy c..with the dimension of the block matrices being c..n x n (n is any number greater than 1)

SUBROUTINE NBTRIP(A,B,C,D,IL,IU,ORDER) INTEGER ORDER,ORDSQ DIMENSION A(1),B(1),D(1)

C..A = SUB DIAG. MATRIX C..B = DIAG. MATRIX C..C = SUP DIAG. MATRIX C..D = RIGHT HAND SIDE OF VECTOR C..IL = LOWER VALUE OF INDEX FOR WHICH MATRICES ARE DEFINED C..IU = UPPER VALUE OF ' ' ' ' C.. (SOLUTION IS SOUGHT FOR BTRI(A,B,C)*X=D C.. FOR INDICES OF X BETWEEN IL AND IU (INCLUSIVE). C.. SOLUTION WRITTEN IN D VECTOR (ORIGINAL CONTENTS C.. ARE OVERWRITTEN)). C..ORDER = ORDER OF A,B,C MATRICES AND LENGTH OF D VECTOR C.. AT EACH POINT DENOTED BY INDEX I C.. C..THE MATRICES AND VECTORS ARE STORED IN SINGLE C..SUBSCRIPT FORM

ORDSQ = ORDER**2

C.. C.. FORWARD ELIMINATION C..

I=IL

IOMAT = 1+(I-1)*ORDSQ

IOVEC = 1+(I-2)*ORDER

CALL LUDECO(B(IOMAT),ORDER)

CALL LUSOLV(B(IOMAT),D(IOVEC),D(IOVEC),ORDER)

DO 100 J=1,ORDER

IOMATJ = IOMAT+(J-1)*ORDER

CALL LUSOLV(B(IOMAT),C(IOMAT),C(IOMATJ),ORDER) 100 CONTINUE 200 CONTINUE

I=I+1

IOMAT = 1+(I-1)*ORDSQ

IOVEC = 1+(I-1)*ORDER

I1MAT = IOMAT-ORDSQ

I1VEC = IOVEC-ORDER

CALL MULPUT(A(IOMAT),D(I1VEC),D(IOVEC),ORDER)

DO 300 J=1,ORDER

IOMATJ = IOMAT+(J-1)*ORDER

I1MATJ = I1MAT+(J-1)*ORDER

CALL MULPUT(A(IOMAT),C(I1MATJ),B(IOMATJ),ORDER) 300 CONTINUE

CALL LUDECO(B(IOMAT),ORDER)

CALL LUSOLV(B(IOMAT),D(IOVEC),D(IOVEC),ORDER)

IF(I.EQ.IU) GOTO 500

DO 400 J=1,ORDER

IOMATJ = IOMAT+(J-1)*ORDER

CALL LUSOLV(B(IOMAT),C(IOMATJ),C(IOMATJ),ORDER) 400 CONTINUE

GOTO 200 500 CONTINUE

C.. C.. BACK SUBSTITUTION C..

I=IU 600 CONTINUE

I=I-1

IOMAT=1+(I-1)*ORDSQ

IOVEC=1+(I-1)*ORDER

I1VEC=IOVEC+ORDER

CALL MULPUT(C(IOMAT),D(I1VEC),D(IOVEC),ORDER)

IF(I.GT.IL) GOTO 600 C..

RETURN

END

C.. SUBROUTINE TO CALCULATE L-U DECOMPOSITION C.. OF A GIVEN MATRIX A AND STORE RESULT IN A

SUBROUTINE LUDECO(A,ORDER)

INTEGER ORDER

DIMENSION A(ORDER,1)

DO 8 JC=2,ORDER 8 A(1,JC) = A(1,JC)/A(1,1)

JRJC = 1 10 CONTINUE

JRJC = JRJC+1

JRJCM1 = JRJC-1

JRJCP1 = JRJC+1

DO 14 JR=JRJC,ORDER

SUM = A(JR,JRJC)

DO 12 JM=1,JRJCM1 12 SUM = SUM-A(JR,JM)*A(JM,JRJC) 14 A(JR,JRJC) = SUM

IF (JRJC.EQ.ORDER) RETURN

DO 18 JC=JRJCP1,ORDER

SUM = A(JRJC,JC)

DO 16 JM=1,JRJCM1 16 SUM = SUM-A(JRJC,JM)*A(JM,JC) 18 A(JRJC,JC) = SUM/A(JRJC,JRJC)

GOTO 10

END

C.. SUBROUTINE TO MULTIPLY A VECTOR B BY A MATRIX A C.. SUBTRACT RESULT FROM ANOTHER VECTOR C AND STORE C.. RESULT IN C. THUS VECTOR C IS OVERWRITTEN

SUBROUTINE MULPUT(A,B,C,ORDER)

INTEGER ORDER

DIMENSION A(1),B(1),C(1)

DO 200 JR=1,ORDER

SUM = 0.0

DO 100 JC=1,ORDER

IA = JR+(JC-1)*ORDER 100 SUM = SUM+A(IA)*B(JC) 200 C(JR) = C(JR)-SUM

RETURN

END

C.. SUBROUTINE TO SOLVE LINEAR ALGEBRAIC SYSTEM OF C.. EQUATION A*C=B AND STORE RESULTS IN VECTOR C C.. MATRIX A INPUT IN L-U DECOMPOSITION FORM

SUBROUTINE LUSOLV(A,B,C,ORDER)

INTEGER ORDER

DIMENSION A(ORDER,1),B(1),C(1)

C.. FIRST L(INV)*B

C(1) = C(1)/A(1,1)

DO 14 JR=2,ORDER

JRM1 = JR-1

SUM = B(JR)

DO 12 JM=1,JRM1 12 SUM = SUM-A(JR,JM)*C(JM) 14 C(JR) = SUM/A(JR,JR)

C.. NEXT U(INV) OF L(INV)*B

DO 18 JRJR=1,ORDER

JR = ORDER,JRJR+1

JRP1 = JR+1

SUM = C(JR)

DO 16 JMJM = JRP1,ORDER

JM = ORDER-JMJM+JRP1 16 SUM = SUM-A(JR,JM)*C(JM) 18 C(JR) = SUM

RETURN

END

I guess if I get really bored tomorrow, I'll post the periodic subroutine, which uses the same sub..subroutines.

  Reply With Quote

Old   February 27, 2007, 02:50
Default Re: Subroutine for Block tridiagonal System
  #3
UT CAA guy
Guest
 
Posts: n/a
wow, that posted really funky, i hope you can read it
  Reply With Quote

Old   February 27, 2007, 03:16
Default Re: Subroutine for Block tridiagonal System
  #4
khan
Guest
 
Posts: n/a
Hi! UT CAA

Thanks for positing the code. have you tested it or not??

in a previous message at this forum, there is a discussuion about this code, which says

" Both of subroutine ludeco and lusolv have the same error. It appears that the code has not been tested in the form as in the book before publication of the book."

but no one point out the exact error. Are you sure about it correctness ?

Khan

  Reply With Quote

Old   February 27, 2007, 11:41
Default Re: Subroutine for Block tridiagonal System
  #5
UT CAA guy
Guest
 
Posts: n/a
Oh, I wasn't sure that it was the same code as in the Anderson book...guess it makes sense that he's the second author.

I will give it a look, but you should have a general idea of the process, and the subroutines aren't 'that' long.

So, as and cfd person would, debug, debug, debug
  Reply With Quote

Old   July 6, 2011, 04:44
Default
  #6
New Member
 
vijay kadli
Join Date: Jul 2011
Posts: 3
Rep Power: 5
vijay kadli is on a distinguished road
hi frnd.. i am solving euler implicit equation and i went thru this btdma solver.. but cudn debug!! can u help me?
vijay kadli is offline   Reply With Quote

Old   July 8, 2011, 15:16
Default
  #7
Member
 
Join Date: Jul 2011
Location: US
Posts: 39
Rep Power: 5
Docfreezzzz is on a distinguished road
I have source code up on the following site if you would like it.
Block Tridiagonal Solver
__________________
CFD engineering resource

Last edited by Docfreezzzz; July 22, 2011 at 16:37. Reason: Eliminate direct code paste
Docfreezzzz is offline   Reply With Quote

Old   July 11, 2011, 07:49
Default btdma
  #8
New Member
 
vijay kadli
Join Date: Jul 2011
Posts: 3
Rep Power: 5
vijay kadli is on a distinguished road
thank u.. wil try the same logic with fortran
thanks once again
vijay kadli is offline   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
CFX11 + Fortran compiler ? Mohan CFX 20 March 30, 2011 19:56
block tridiagonal system solver for periodic bc Abdullah Main CFD Forum 8 February 23, 2007 14:15
Tridiagonal System Zhong Lei Main CFD Forum 3 May 24, 2001 15:18
On the block tridiagonal linear system zhanglei Main CFD Forum 2 July 24, 2000 07:15
On the block tridiagonal linear system zhanglei Main CFD Forum 0 July 18, 2000 09:55


All times are GMT -4. The time now is 16:07.