C C REAL TRIDIAGONAL SOLUTION BY POWER-OF-THREE CYCLIC REDUCTION C C CALL TRICYC(A,B,C,D,N) C C SOLVES TX=D WHERE TRIDIAGONAL MATRIX HAS THE FORM C C | B0 C0 | C | A1 B1 C1 | C | A2 B2 C2 | C | A3 B3 C3 | C | | C | | C | | C | | C | AN-1 BN-1 | C C C ANSWER OVERWRITES "D". ALL INPUTS WILL BE OVERWRITTEN C DURING CALCULATION. A0 AND CN-1 WILL BE MODIFIED. C A AND C MAY NOT BE EQUIVALENCED. C C AUTHOR: STEWART A. LEVIN MOBIL RESEARCH AND DEVELOPMENT CORP. C DATE: JULY 28, 1989 C C SUBROUTINE TRICYC(A,B,C,D,N) IMPLICIT NONE INTEGER N REAL A(0:N-1) ! SUBDIAGONAL COEFFICIENTS REAL B(0:N-1) ! DIAGONAL COEFFICIENTS REAL C(0:N-1) ! SUPERDIAGONAL COEFFICIENTS REAL D(0:N-1) ! RIGHT-HAND SIDE VECTOR C C LOCAL VARIABLES (STACK) C INTEGER MAXSTK ! MAXIMUM STACK DEPTH PARAMETER (MAXSTK=20) INTEGER ISHORT ! LENGTH AT WHICH TO SWITCH TO SCALAR SOLVER PARAMETER (ISHORT=25) ! EXTERNAL SCTRIL INTEGER ISTK ! STACK POINTER INTEGER NM1 ! CURRENT VALUE OF N-1 INTEGER ND3 ! CURRENT VALUE OF (N-1)/3 INTEGER NREM(MAXSTK) ! CURRENT AND PRIOR VALUES OF (N-1) MOD 3 INTEGER I,J ! STRIDE AND LOOP INCREMENT REAL TA110,TA210 ! LOOP TEMPORARIES REAL TA130,TA230 ! "" REAL TA150,TA250 ! "" REAL TA020,TC020 ! "" REAL TA040,TC040 ! "" REAL TA060,TC060 ! "" REAL TB110,TB210 ! "" REAL TB130,TB230 ! "" REAL TB150,TB250 ! "" REAL TA061 ! "" REAL TC110,TC210 ! "" REAL TC130,TC230 ! "" REAL TC150,TC250 ! "" REAL TD110,TD210 ! "" REAL TD130,TD230 ! "" REAL TD150,TD250 ! "" REAL TEMP10,TEMP30,TEMP50 ! "" REAL TEMP31 ! "" C C INITIALIZE REDUCTION C ISTK=1 NM1=N-1 I=1 ! INITIALLY COMPACT VECTORS A(0)=0.0 ! CONVENIENCE C(NM1)=0.0 ! "" C C DO WHILE(NM1.GT.ISHORT) 101 CONTINUE IF(NM1.LE.ISHORT)GOTO 102 ND3=NM1/3 NREM(ISTK)=NM1-3*ND3 IF(NREM(ISTK).EQ.2) THEN CDIR$ IVDEP C$DIR NO_RECURRENCE DO 10 J=0,ND3 TB110=B((3*J+1)*I) TB210=B((3*J+2)*I) TA210=A((3*J+2)*I) TC110=C((3*J+1)*I) TEMP10=1.0/(TB110*TB210-TC110*TA210) TB110=TEMP10*TB110 TB210=TEMP10*TB210 TA210=TEMP10*TA210 TC110=TEMP10*TC110 TD110=D((3*J+1)*I) TD210=D((3*J+2)*I) D((3*J+1)*I)=TD110*TB210-TC110*TD210 D((3*J+2)*I)=TD210*TB110-TA210*TD110 TA110=A((3*J+1)*I) A((3*J+2)*I)=TA110*TA210 A((3*J+1)*I)=TA110*TB210 TC210=C((3*J+2)*I) C((3*J+1)*I)=TC110*TC210 C((3*J+2)*I)=TC210*TB110 10 CONTINUE C B(0)=B(0)-C(0)*A(I) D(0)=D(0)-C(0)*D(I) C(0)=C(0)*C(I) CDIR$ IVDEP C$DIR NO_RECURRENCE DO 20 J=1,ND3 TC020=C((3*J)*I) TA020=A((3*J)*I) B((3*J)*I)=B((3*J)*I)-TC020*A((3*J+1)*I) 1 -TA020*C((3*J-1)*I) D((3*J)*I)=D((3*J)*I)-TC020*D((3*J+1)*I) 1 -TA020*D((3*J-1)*I) A((3*J)*I)=TA020*A((3*J-1)*I) C((3*J)*I)=TC020*C((3*J+1)*I) 20 CONTINUE ENDIF IF(NREM(ISTK).EQ.1) THEN CDIR$ IVDEP C$DIR NO_RECURRENCE DO 30 J=0,ND3-1 TB130=B((3*J+1)*I) TB230=B((3*J+2)*I) TA230=A((3*J+2)*I) TC130=C((3*J+1)*I) TEMP30=1.0/(TB130*TB230-TC130*TA230) TB130=TEMP30*TB130 TB230=TEMP30*TB230 TA230=TEMP30*TA230 TC130=TEMP30*TC130 TD130=D((3*J+1)*I) TD230=D((3*J+2)*I) D((3*J+1)*I)=TD130*TB230-TC130*TD230 D((3*J+2)*I)=TD230*TB130-TA230*TD130 TA130=A((3*J+1)*I) A((3*J+2)*I)=TA130*TA230 A((3*J+1)*I)=TA130*TB230 TC230=C((3*J+2)*I) C((3*J+1)*I)=TC130*TC230 C((3*J+2)*I)=TC230*TB130 30 CONTINUE TEMP31=1.0/B(NM1*I) A(NM1*I)=TEMP31*A(NM1*I) D(NM1*I)=TEMP31*D(NM1*I) C B(0)=B(0)-C(0)*A(I) D(0)=D(0)-C(0)*D(I) C(0)=C(0)*C(I) CDIR$ IVDEP C$DIR NO_RECURRENCE DO 40 J=1,ND3 TC040=C((3*J)*I) TA040=A((3*J)*I) B((3*J)*I)=B((3*J)*I)-TC040*A((3*J+1)*I) 1 -TA040*C((3*J-1)*I) D((3*J)*I)=D((3*J)*I)-TC040*D((3*J+1)*I) 1 -TA040*D((3*J-1)*I) A((3*J)*I)=TA040*A((3*J-1)*I) C((3*J)*I)=TC040*C((3*J+1)*I) 40 CONTINUE ENDIF IF(NREM(ISTK).EQ.0) THEN CDIR$ IVDEP C$DIR NO_RECURRENCE DO 50 J=0,ND3-1 TB150=B((3*J+1)*I) TB250=B((3*J+2)*I) TA250=A((3*J+2)*I) TC150=C((3*J+1)*I) TEMP50=1.0/(TB150*TB250-TC150*TA250) TB150=TEMP50*TB150 TB250=TEMP50*TB250 TA250=TEMP50*TA250 TC150=TEMP50*TC150 TD150=D((3*J+1)*I) TD250=D((3*J+2)*I) D((3*J+1)*I)=TD150*TB250-TC150*TD250 D((3*J+2)*I)=TD250*TB150-TA250*TD150 TA150=A((3*J+1)*I) A((3*J+2)*I)=TA150*TA250 A((3*J+1)*I)=TA150*TB250 TC250=C((3*J+2)*I) C((3*J+1)*I)=TC150*TC250 C((3*J+2)*I)=TC250*TB150 50 CONTINUE C B(0)=B(0)-C(0)*A(I) D(0)=D(0)-C(0)*D(I) C(0)=C(0)*C(I) CDIR$ IVDEP C$DIR NO_RECURRENCE DO 60 J=1,ND3-1 TC060=C((3*J)*I) TA060=A((3*J)*I) B((3*J)*I)=B((3*J)*I)-TC060*A((3*J+1)*I) 1 -TA060*C((3*J-1)*I) D((3*J)*I)=D((3*J)*I)-TC060*D((3*J+1)*I) 1 -TA060*D((3*J-1)*I) A((3*J)*I)=TA060*A((3*J-1)*I) C((3*J)*I)=TC060*C((3*J+1)*I) 60 CONTINUE TA061=A(NM1*I) B(NM1*I)=B(NM1*I)-TA061*C((NM1-1)*I) D(NM1*I)=D(NM1*I)-TA061*D((NM1-1)*I) A(NM1*I)=TA061*A((NM1-1)*I) ENDIF NM1=ND3 ISTK=ISTK+1 I=I*3 GOTO 101 102 CONTINUE C C SOLVE SMALL SYSTEM C CALL SCTRIC(A,B,C,D,I,NM1) C C END REDUCE PHASE C C BACK SUBSTITUTE C C DO WHILE(ISTK.GT.1) 201 CONTINUE IF(ISTK.LE.1) GOTO 202 ISTK=ISTK-1 I=I/3 ND3=NM1 NM1=3*ND3+NREM(ISTK) IF(NREM(ISTK).EQ.2) THEN CDIR$ IVDEP C$DIR NO_RECURRENCE DO 70 J=0,ND3-1 D((3*J+1)*I)=D((3*J+1)*I) 1 -A((3*J+1)*I)*D((3*J )*I) 2 +C((3*J+1)*I)*D((3*J+3)*I) D((3*J+2)*I)=D((3*J+2)*I) 1 +A((3*J+2)*I)*D((3*J )*I) 2 -C((3*J+2)*I)*D((3*J+3)*I) 70 CONTINUE D((NM1-1)*I)=D((NM1-1)*I)-A((NM1-1)*I)*D((NM1-2)*I) D(NM1*I)=D(NM1*I)+A(NM1*I)*D((NM1-2)*I) ENDIF IF(NREM(ISTK).EQ.1) THEN CDIR$ IVDEP C$DIR NO_RECURRENCE DO 80 J=0,ND3-1 D((3*J+1)*I)=D((3*J+1)*I) 1 -A((3*J+1)*I)*D((3*J )*I) 2 +C((3*J+1)*I)*D((3*J+3)*I) D((3*J+2)*I)=D((3*J+2)*I) 1 +A((3*J+2)*I)*D((3*J )*I) 2 -C((3*J+2)*I)*D((3*J+3)*I) 80 CONTINUE D(NM1*I)=D(NM1*I)-A(NM1*I)*D((NM1-1)*I) ENDIF IF(NREM(ISTK).EQ.0) THEN CDIR$ IVDEP C$DIR NO_RECURRENCE DO 90 J=0,ND3-1 D((3*J+1)*I)=D((3*J+1)*I) 1 -A((3*J+1)*I)*D((3*J )*I) 2 +C((3*J+1)*I)*D((3*J+3)*I) D((3*J+2)*I)=D((3*J+2)*I) 1 +A((3*J+2)*I)*D((3*J )*I) 2 -C((3*J+2)*I)*D((3*J+3)*I) 90 CONTINUE ENDIF GOTO 201 202 CONTINUE C C END BACK SUBSTITUTION PHASE C RETURN END C C SCALAR TRIDIAGONAL SOLVER WITH ARBITRARY INTER-ELEMENT SPACING C SUBROUTINE SCTRIC(A,B,C,D,I,NM1) IMPLICIT NONE INTEGER I,NM1 REAL A(I,0:NM1) REAL B(I,0:NM1) REAL C(I,0:NM1) REAL D(I,0:NM1) INTEGER J REAL TEMP,TAMP TAMP=1.0/B(1,0) D(1,0)=TAMP*D(1,0) IF(NM1.EQ.0) RETURN C(1,0)=TAMP*C(1,0) C CDIR$ NOVECTOR C$DIR SCALAR DO 10 J=1,NM1-1,1 TEMP=1.0/(B(1,J)-A(1,J)*C(1,J-1)) D(1,J)=TEMP*(D(1,J)-A(1,J)*D(1,J-1)) C(1,J)=TEMP*C(1,J) 10 CONTINUE TAMP=1.0/(B(1,NM1)-A(1,NM1)*C(1,NM1-1)) D(1,NM1)=TAMP*(D(1,NM1)-A(1,NM1)*D(1,NM1-1)) C CDIR$ NOVECTOR C$DIR SCALAR DO 20 J=NM1-1,0,-1 D(1,J)=D(1,J)-C(1,J)*D(1,J+1) 20 CONTINUE RETURN END