lubksb.f

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

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

Purpose

LUBKSB solves $A\,x = b$ given the LU decomposition of $A$ from LUDCMP (Numerical Recipes' LUBKSB). The solution overwrites B.

How it fits

Called only by THDIFF. No COMMON blocks, no calls.

The method

Forward substitution with $L$ (unpacking the permutation INDX, and skipping leading zeros of $b$), then back-substitution with $U$:

$$y_i = b_{p(i)} - \sum_{j<i} L_{ij}\,y_j,\qquad x_i = \frac{1}{U_{ii}}\Big(y_i - \sum_{j>i} U_{ij}\,x_j\Big).$$

Definitions

subroutineLUBKSB

lines 3–51 · GitHub

Forward and back substitution (above).

Called by THDIFF ×3

Compiler warnings (1)

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

line 33: Inequality 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********************************************************************* 3 SUBROUTINE LUBKSB(A,N,NP,INDX,B) 4C 10/20/26 MHP via Claude 5C code changed to replace implicit typing with IMPLICIT NONE and 6C explicit declarations. 7 IMPLICIT NONE 8 REAL*8 A,B,SUM 9 10 11C .. Scalar Arguments .. 12 INTEGER N,NP 13C .. 14C .. Array Arguments .. 15 DIMENSION A(NP,NP),B(N) 16 INTEGER INDX(N) 17C .. 18C .. Local Scalars .. 19C REAL SUM 20 INTEGER I,II,J,LL 21 SAVE 22C .. 23 II = 0 24 DO 12 I = 1,N 25 LL = INDX(I) 26 SUM = B(LL) 27 B(LL) = B(I) 28 IF (II.NE.0) THEN 29 DO 11 J = II,I - 1 30 SUM = SUM - A(I,J)*B(J) 31 11 11 CONTINUE 32 33 ELSE IF (SUM.NE.0D0) THEN 34 II = I 35 END IF 36 37 B(I) = SUM 38 12 12 CONTINUE 39 DO 14 I = N,1,-1 40 SUM = B(I) 41 IF (I.LT.N) THEN 42 DO 13 J = I + 1,N 43 SUM = SUM - A(I,J)*B(J) 44 13 13 CONTINUE 45 END IF 46 47 B(I) = SUM/A(I,I) 48 14 14 CONTINUE 49 RETURN 50 51 END

1"""Python version of lubksb.f (YREC development, 7be4133).""" 2 3import numpy as np 4 5 6def lubksb(a: np.ndarray, n: int, np_: int, indx: np.ndarray, b: np.ndarray) -> None: 7 """Solve A x = b from the LU decomposition made by LUDCMP. 8 9 *Numerical Recipes*' LUBKSB: forward substitution with L (unscrambling 10 the row permutation ``indx`` and skipping the leading zeros of b), then 11 back-substitution with U. The solution overwrites ``b``. 12 13 Translates: lubksb.f#LUBKSB (7be4133) 14 15 Parameters 16 ---------- 17 a : ndarray, shape (np_ + 1, np_ + 1) 18 The LU decomposition from LUDCMP, ``a[i, j]`` for ``A(I, J)``. 19 n : int 20 The order of the system. 21 np_ : int 22 The declared size of ``a`` (Fortran's ``NP``; renamed so as not 23 to hide NumPy's usual name). 24 indx : ndarray of int, shape (n + 1,) 25 The row permutation from LUDCMP. 26 b : ndarray, shape (n + 1,) 27 In: the right-hand side. Out, in place: the solution. 28 29 """ 30 31 # Forward substitution with L (loop 12, Fortran lines 23–38). 32 ii = 0 33 for i in range(1, n + 1): 34 ll = indx[i] 35 sum_ = b[ll] 36 b[ll] = b[i] 37 if ii != 0: 38 for j in range(ii, i): # loop 11 39 sum_ = sum_ - a[i, j] * b[j] 40 elif sum_ != 0.0: 41 ii = i 42 b[i] = sum_ 43 44 # Back substitution with U (loop 14, Fortran lines 39–48). 45 for i in range(n, 0, -1): 46 sum_ = b[i] 47 if i < n: 48 for j in range(i + 1, n + 1): # loop 13 49 sum_ = sum_ - a[i, j] * b[j] 50 b[i] = sum_ / a[i, i]