tridiag_gs.f

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

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

Variables on this page (1)

json

Purpose

TRIDIAG_GS solves a tridiagonal system $a_j x_{j-1} + b_j x_j + c_j x_{j+1} = e'_j$, $j = 1, \dots, n$ (the Thomas algorithm), with the matrix and right-hand side as arguments. A zero pivot stops the run, with the message "#TRIDIA:SINGULAR MATRIX."

How it fits

Who calls it: GRSETT (the new run of $X$ in gravitational settling, by its comment) and MICRODIFF_RUN (the new run of Li).
Output: the solution EX.
Modules: DIMENSIONS_MOD (the array size JSON) only.

Definitions

subroutineTRIDIAG_GS

lines 2–31 · GitHub

The Thomas algorithm for a tridiagonal system, with arrays of size JSON (5000, from DIMENSIONS_MOD) as arguments.

Uses modules DIMENSIONS_MOD

Called by GRSETT, MICRODIFF_RUN

Compiler warnings (1)

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

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

Change log

The changes since rp_1

Changed since rp_1. Declarations only: JSON from DIMENSIONS_MOD, IMPLICIT NONE.

1C $$$$$$ 2 SUBROUTINE TRIDIAG_GS(A,B,C,EX_PRIME,NPT,EX) 3C 10/15/26 MHP via Claude 4C code changed to take JSON from DIMENSIONS_MOD instead of declaring 5C PARAMETER(JSON=5000) here. 6 USE DIMENSIONS_MOD 7C 10/20/26 MHP via Claude 8C code changed to replace implicit typing with IMPLICIT NONE and 9C explicit declarations. 10 IMPLICIT NONE 11 REAL*8 A,B,BET,C,EX,EX_PRIME,GAMA 12 INTEGER J,NPT 13 DIMENSION A(JSON),B(JSON),C(JSON),EX_PRIME(JSON),EX(JSON), 14 * * GAMA(JSON) 15 SAVE 16 17 BET = B(1) 18 EX(1) = EX_PRIME(1)/BET 19 DO 1 J=2,NPT 20 GAMA(J) = C(J-1)/BET 21 BET = B(J) - A(J)*GAMA(J) 22 IF( BET.EQ.0 ) STOP '#TRIDIA:SINGULAR MATRIX' 23 EX(J) = (EX_PRIME(J) - A(J)*EX(J-1))/BET 24 1 1 CONTINUE 25 DO 2 J=NPT-1,1,-1 26 EX(J) = EX(J) - GAMA(J+1)*EX(J+1) 27 2 2 CONTINUE 28C 911 FORMAT(1P6E10.2) 29 30 RETURN 31 END

1"""Python version of tridiag_gs.f (YREC 7be4133): the general tridiagonal 2solver (Thomas algorithm), used by gravitational settling and the Li run.""" 3 4import numpy as np 5 6from pyrec import dimensions_mod 7 8# The Fortran's local work array (SAVE; every element read is set first). 9gama = np.zeros(dimensions_mod.json + 1) 10 11 12def tridiag_gs(a: np.ndarray, b: np.ndarray, c: np.ndarray, ex_prime: np.ndarray, npt: int, ex: np.ndarray) -> None: 13 """Solve a_j x_{j-1} + b_j x_j + c_j x_{j+1} = ex_prime_j, j = 1..npt, 14 by forward elimination and back substitution, without pivoting. 15 16 Translates: tridiag_gs.f#TRIDIAG_GS (7be4133) 17 18 Arguments: the diagonals a, b, c and the right-hand side ex_prime 19 (arrays, index 0 unused); npt, the number of points; ex, the solution, 20 set in place. None is returned. 21 22 """ 23 24 # Forward elimination (Fortran lines 17–24). 25 bet = b[1] 26 ex[1] = ex_prime[1] / bet 27 for j in range(2, npt + 1): # loop 1 28 gama[j] = c[j - 1] / bet 29 bet = b[j] - a[j] * gama[j] 30 if bet == 0: 31 raise SystemExit('#TRIDIA:SINGULAR MATRIX') 32 ex[j] = (ex_prime[j] - a[j] * ex[j - 1]) / bet 33 34 # Back substitution (Fortran lines 25–27). 35 for j in range(npt - 1, 0, -1): # loop 2 36 ex[j] = ex[j] - gama[j + 1] * ex[j + 1]