ludcmp.f

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

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

Purpose

LUDCMP replaces a square matrix by its LU decomposition with partial pivoting (Crout's method with implicit scaling), as Numerical Recipes' LUDCMP. LUBKSB then solves systems with it.

How it fits

Called only by THDIFF (the diffusion coefficients of Thoul et al., by its name). No COMMON blocks, no calls. Arguments: A(NP,NP) (overwritten by L and U), N, NP, the row permutation INDX, and D = ±1 (the parity of the row swaps).

The method

It first stores each row's scale, $1/\max_k |a_{ik}|$, and stops with STOP 'Singular matrix.' if a row is all zero. Then, for each column $j$ it computes the $U$ entries above the diagonal and the candidates below, picks the pivot with the largest scaled magnitude $|a_{ij}|/\max_k |a_{ik}|$, swaps rows, and divides the sub-diagonal entries by the pivot. A zero pivot is replaced by TINY = $10^{-20}$.

The scale vector VV has the fixed size NMAX = 100 (line 13 ↗; line 28 ↗), and N is not checked against it. THDIFF passes N = $2M + 2$ with its own arrays sized 42 (thdiff.f, lines 48–54 ↗; thdiff.f, line 72 ↗; thdiff.f, line 264 ↗), so N ≤ 42; KNOWN_ISSUES.md (lines 1845–1848) keeps NMAX per routine on purpose (settled from the code, Nov 9, 2026).

Definitions

subroutineLUDCMP

lines 2–97 · GitHub

LU decomposition with partial pivoting (Numerical Recipes).

Called by THDIFF

Compiler warnings (3)

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

line 40: Equality comparison for REAL(8) [-Wcompare-reals]
line 86: Equality comparison for REAL(8) [-Wcompare-reals]
line 94: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE with explicit declarations.

1********************************************************* 2 SUBROUTINE LUDCMP(A,N,NP,INDX,D) 3C 10/20/26 MHP via Claude 4C code changed to replace implicit typing with IMPLICIT NONE and 5C explicit declarations. 6 IMPLICIT NONE 7 REAL*8 A,AAMAX,D,DUM,SUM,TINY,VV 8 9 10C .. Parameters .. 11 INTEGER NMAX 12C REAL TINY 13 PARAMETER (NMAX=100,TINY=1.0D-20) 14C .. 15C .. Scalar Arguments .. 16C REAL D 17 INTEGER N,NP 18C .. 19C .. Array Arguments .. 20 DIMENSION A(NP,NP) 21 INTEGER INDX(N) 22C .. 23C .. Local Scalars .. 24C REAL AAMAX,DUM,SUM 25 INTEGER I,IMAX,J,K 26C .. 27C .. Local Arrays .. 28 DIMENSION VV(NMAX) 29 SAVE 30C .. 31C .. Intrinsic Functions .. 32 INTRINSIC ABS 33C .. 34 D = 1.D0 35 DO 12 I = 1,N 36 AAMAX = 0.D0 37 DO 11 J = 1,N 38 IF (ABS(A(I,J)).GT.AAMAX) AAMAX = ABS(A(I,J)) 39 11 11 CONTINUE 40 IF (AAMAX.EQ.0D0) STOP 'Singular matrix.' 41 VV(I) = 1.D0/AAMAX 42 12 12 CONTINUE 43 DO 19 J = 1,N 44 IF (J.GT.1) THEN 45 DO 14 I = 1,J - 1 46 SUM = A(I,J) 47 IF (I.GT.1) THEN 48 DO 13 K = 1,I - 1 49 SUM = SUM - A(I,K)*A(K,J) 50 13 13 CONTINUE 51 A(I,J) = SUM 52 END IF 53 54 14 14 CONTINUE 55 END IF 56 57 AAMAX = 0.D0 58 DO 16 I = J,N 59 SUM = A(I,J) 60 IF (J.GT.1) THEN 61 DO 15 K = 1,J - 1 62 SUM = SUM - A(I,K)*A(K,J) 63 15 15 CONTINUE 64 A(I,J) = SUM 65 END IF 66 67 DUM = VV(I)*ABS(SUM) 68 IF (DUM.GE.AAMAX) THEN 69 IMAX = I 70 AAMAX = DUM 71 END IF 72 73 16 16 CONTINUE 74 IF (J.NE.IMAX) THEN 75 DO 17 K = 1,N 76 DUM = A(IMAX,K) 77 A(IMAX,K) = A(J,K) 78 A(J,K) = DUM 79 17 17 CONTINUE 80 D = -D 81 VV(IMAX) = VV(J) 82 END IF 83 84 INDX(J) = IMAX 85 IF (J.NE.N) THEN 86 IF (A(J,J).EQ.0D0) A(J,J) = TINY 87 DUM = 1.D0/A(J,J) 88 DO 18 I = J + 1,N 89 A(I,J) = A(I,J)*DUM 90 18 18 CONTINUE 91 END IF 92 93 19 19 CONTINUE 94 IF (A(N,N).EQ.0D0) A(N,N) = TINY 95 RETURN 96 97 END

1"""Python version of ludcmp.f (YREC development, 7be4133).""" 2 3import numpy as np 4 5NMAX = 100 6TINY = 1.0e-20 7 8 9def ludcmp(a: np.ndarray, n: int, np_: int, indx: np.ndarray, d: float) -> float: 10 """LU decomposition with partial pivoting (Crout, implicit scaling). 11 12 *Numerical Recipes*' LUDCMP: ``a`` is replaced by L (below the 13 diagonal, unit diagonal implied) and U; ``indx`` records the row 14 swaps and ``d`` their parity. LUBKSB then solves systems with it. 15 16 Translates: ludcmp.f#LUDCMP (7be4133) 17 18 Parameters 19 ---------- 20 a : ndarray, shape (np_ + 1, np_ + 1) 21 In: the matrix, ``a[i, j]`` for ``A(I, J)``. Out, in place: its 22 LU decomposition. 23 n : int 24 The order of the matrix. 25 np_ : int 26 The declared size of ``a`` (Fortran's ``NP``). 27 indx : ndarray of int, shape (n + 1,) 28 Set in place: the row permutation. 29 d : float 30 Output only (passed in with any value). 31 32 Returns 33 ------- 34 d : float 35 +1 or -1, the parity of the row interchanges. 36 37 """ 38 39 # The scale vector VV has a fixed size NMAX, and n is not checked 40 # against it (THDIFF passes n <= 42). 41 vv = np.zeros(NMAX + 1) 42 imax = 0 43 44 # Each row's scale (Fortran lines 34–42). 45 d = 1.0 46 for i in range(1, n + 1): # loop 12 47 aamax = 0.0 48 for j in range(1, n + 1): # loop 11 49 if abs(a[i, j]) > aamax: 50 aamax = abs(a[i, j]) 51 if aamax == 0.0: 52 raise SystemExit("Singular matrix.") 53 vv[i] = 1.0 / aamax 54 55 # Crout's method, column by column (loop 19, Fortran lines 43–93). 56 for j in range(1, n + 1): 57 # U above the diagonal (Fortran lines 44–55). 58 if j > 1: 59 for i in range(1, j): # loop 14 60 sum_ = a[i, j] 61 if i > 1: 62 for k in range(1, i): # loop 13 63 sum_ = sum_ - a[i, k] * a[k, j] 64 a[i, j] = sum_ 65 66 # The candidates below, and the largest scaled pivot (Fortran 67 # lines 57–73). 68 aamax = 0.0 69 for i in range(j, n + 1): # loop 16 70 sum_ = a[i, j] 71 if j > 1: 72 for k in range(1, j): # loop 15 73 sum_ = sum_ - a[i, k] * a[k, j] 74 a[i, j] = sum_ 75 dum = vv[i] * abs(sum_) 76 if dum >= aamax: 77 imax = i 78 aamax = dum 79 80 # Swap rows (Fortran lines 74–82). 81 if j != imax: 82 for k in range(1, n + 1): # loop 17 83 dum = a[imax, k] 84 a[imax, k] = a[j, k] 85 a[j, k] = dum 86 d = -d 87 vv[imax] = vv[j] 88 89 # Divide by the pivot, TINY if zero (Fortran lines 84–91). 90 indx[j] = imax 91 if j != n: 92 if a[j, j] == 0.0: 93 a[j, j] = TINY 94 dum = 1.0 / a[j, j] 95 for i in range(j + 1, n + 1): # loop 18 96 a[i, j] = a[i, j] * dum 97 98 # Fortran line 94. 99 if a[n, n] == 0.0: 100 a[n, n] = TINY 101 return d