# Subroutine for Block tridiagonal System

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

 February 26, 2007, 22:39 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

 February 27, 2007, 02:49 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.

 February 27, 2007, 02:50 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

 February 27, 2007, 03:16 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

 February 27, 2007, 11:41 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

 July 6, 2011, 03:44 #6 New Member   vijay kadli Join Date: Jul 2011 Posts: 3 Rep Power: 7 hi frnd.. i am solving euler implicit equation and i went thru this btdma solver.. but cudn debug!! can u help me?

 July 8, 2011, 14:16 #7 Member   Join Date: Jul 2011 Location: US Posts: 39 Rep Power: 7 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 15:37. Reason: Eliminate direct code paste

 July 11, 2011, 06:49 btdma #8 New Member   vijay kadli Join Date: Jul 2011 Posts: 3 Rep Power: 7 thank u.. wil try the same logic with fortran thanks once again

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Mohan CFX 20 March 30, 2011 18:56 Abdullah Main CFD Forum 8 February 23, 2007 14:15 Zhong Lei Main CFD Forum 3 May 24, 2001 14:18 zhanglei Main CFD Forum 2 July 24, 2000 06:15 zhanglei Main CFD Forum 0 July 18, 2000 08:55

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