ctridi.f

32 lines · 1 definitions · on GitHub (7be4133)

Python version: ctridi.py, 38 lines, on the right (switch to Python above the code).

Variables on this page (6)

a, b, c, d, gama, u

Purpose

Solves a tridiagonal linear system by the Thomas algorithm (Numerical Recipes' TRIDAG, per its comment), for the new abundance of a species after rotational mixing.

How it fits

Inputs: N, and the coefficients A, B, C and right-hand side D in TRIDI_MOD, set by CCOEFT.
Who calls it: MIXCOM.
Output: the solution U in TRIDI_MOD (GAMA is a work array there). The same module serves the angular-momentum problem: DCOEFT fills it for TRIDIA, and GETQUA also uses it.

The method

For $a_j u_{j-1} + b_j u_j + c_j u_{j+1} = d_j$: forward elimination with $\beta_1 = b_1$, $\gamma_j = c_{j-1}/\beta_{j-1}$, $\beta_j = b_j - a_j\gamma_j$, $u_j = (d_j - a_j u_{j-1})/\beta_j$; then back substitution $u_j \leftarrow u_j - \gamma_{j+1}u_{j+1}$. No pivoting; a zero $\beta$ stops the run (with the message "TRIDIA: singular matrix").

Design

TRIDIA is the same solver with extra bookkeeping of the change in angular momentum; this copy is for composition.

Definitions

subroutineCTRIDI

lines 2–32 · GitHub

Thomas algorithm for the tridiagonal system in TRIDI_MOD; the solution goes to U.

Uses modules TRIDI_MOD

Called by MIXCOM

Compiler warnings (1)

gfortran -Wall -Wextra, compiling only; nothing was run.

line 24: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: /TRIDI/ to TRIDI_MOD, IMPLICIT NONE.

1C $$$$$$ 2 SUBROUTINE CTRIDI(N) 3C 10/14/26 MHP via Claude 4C code changed to convert COMMON/TRIDI/ into TRIDI_MOD. 5 USE TRIDI_MOD 6C 10/20/26 MHP via Claude 7C code changed to replace implicit typing with IMPLICIT NONE and 8C explicit declarations. 9 IMPLICIT NONE 10 REAL*8 BET 11 INTEGER J,N 12 13 SAVE 14 15C SOLVES A TRIDIAGONAL MATRIX SYSTEM FOR THE FRACTIONAL ABUNDANCE 16C OF A SPECIES. NEW ABUNDANCE IS STORED IN ARRAY U. 17C THIS ROUTINE IS FROM NUMERICAL RECIPES, P.40. 18 19 BET = B(1) 20 U(1) = D(1)/BET 21 DO 1 J=2,N 22 GAMA(J) = C(J-1)/BET 23 BET = B(J) - A(J)*GAMA(J) 24 IF( BET.EQ.0 ) STOP '#TRIDIA:SINGULAR MATRIX' 25 U(J) = (D(J) - A(J)*U(J-1))/BET 26 1 1 CONTINUE 27 DO 2 J=N-1,1,-1 28 U(J) = U(J) - GAMA(J+1)*U(J+1) 29 2 2 CONTINUE 30 31 RETURN 32 END

1"""Python version of ctridi.f (YREC 7be4133): the tridiagonal solver 2(Thomas algorithm) for the composition mixing step.""" 3 4from pyrec import tridi_mod 5 6 7def ctridi(n: int) -> None: 8 """Solve the tridiagonal system a, b, c, d of tridi_mod (set by 9 ccoeft) by forward elimination and back substitution, without 10 pivoting; the solution goes to tridi_mod.u. 11 12 Translates: ctridi.f#CTRIDI (7be4133) 13 14 Arguments: n, the number of points. None is returned. 15 16 """ 17 18 a = tridi_mod.a 19 b = tridi_mod.b 20 c = tridi_mod.c 21 d = tridi_mod.d 22 u = tridi_mod.u 23 gama = tridi_mod.gama 24 25 # Forward elimination (Fortran lines 19–26). 26 bet = b[1] 27 u[1] = d[1] / bet 28 for j in range(2, n + 1): 29 gama[j] = c[j - 1] / bet 30 bet = b[j] - a[j] * gama[j] 31 # Question (see the page): only an exact zero is caught. 32 if bet == 0: 33 raise SystemExit('#TRIDIA:SINGULAR MATRIX') 34 u[j] = (d[j] - a[j] * u[j - 1]) / bet 35 36 # Back substitution (Fortran lines 27–29). 37 for j in range(n - 1, 0, -1): 38 u[j] = u[j] - gama[j + 1] * u[j + 1]