bandw.f

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

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

Purpose

This file solves a banded linear system $A\,\mathbf{x} = \mathbf{b}$ by Gaussian elimination without pivoting. It is the solver for the implicit angular-momentum transport in DADCOEFT (its comments speak of $\Omega$ and its radial derivatives).

How it fits

Inputs: the arguments A(NMAX,10) (the band), NM (the number of unknowns) and B (the right-hand side).
Who calls it: only DADCOEFT.
Output: the solution, returned in B; A is overwritten.
No COMMON blocks.

Storage. Row $i$ of the matrix is stored in A(i,1:10), with the diagonal in column 5: A(i,5+k) is the coefficient of unknown $x_{i+k}$, for $k = -4, \dots, 5$. The band thus has 4 subdiagonals and 5 superdiagonals. A is declared A(NMAX,10) with NMAX = 8000 (line 9 ↗; line 12 ↗), and DADCOEFT declares the same (dadcoeft.f, line 39 ↗; dadcoeft.f, line 92 ↗); the two must agree, since the first dimension sets where each column starts. KNOWN_ISSUES.md (lines 1845–1848) keeps NMAX per routine on purpose.

The method

Forward elimination, row by row: divide row $i$ and $b_i$ by the pivot $a_{ii}$, then remove $x_i$ from the next 4 rows, $$a_{jk} \leftarrow a_{jk} - a_{ji}\,a_{ik}, \qquad b_j \leftarrow b_j - a_{ji}\,b_i, \qquad j = i+1, \dots, i+4.$$ The last 4 rows are handled separately, since their bands are cut off by the end of the matrix. A pivot smaller than $10^{-20}$ in absolute value stops the run (STOP 999), after printing the row (for the last row it stops without printing).
Back substitution: with unit diagonal, for $i = N, \dots, 2$, $$b_{i-l} \leftarrow b_{i-l} - b_i\,a_{i-l,\,i}, \qquad l = 1, \dots, \min(5, i-1).$$

Design

No pivoting: the method relies on the matrix being diagonally dominant enough.

About 100 lines of an earlier version (blocks of 4 rows per shell, "DIAGONALIZE OMEGA(II) TERM") remain as comments.

Definitions

subroutineBANDW

lines 3–186 · GitHub

Gaussian elimination without pivoting for a band matrix with 4 subdiagonals and 5 superdiagonals (diagonal in column 5 of A); the solution replaces B.

Called by DADCOEFT

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE, and the unused JSON parameter removed, which settled half of a former question on the dimensions (the other half is now in Storage, above).

1C $$$$$$ 2C SUBROUTINE BANDW(A,NM,M1,M2,B) ! KC 2025-05-31 3 SUBROUTINE BANDW(A,NM,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 INTEGER II,IMJ,IMJ5,JJ,NMAX 9 PARAMETER (NMAX=8000) 10C INTEGER*4 M1,M2,NM ! KC 2025-05-31 11 INTEGER*4 NM 12 REAL*8 A(NMAX,10),B(NMAX),TINY 13 PARAMETER (TINY=1.D-20) 14 INTEGER*4 I,J,K 15 SAVE 16C INITIAL STEP: DIAGONALIZE TERMS INVOLVING OMEGA(1) 17C AND THEN STEP THROUGH ITS FIRST THROUGH THIRD 18C DERIVATIVES. THEN PROCEED TO THE SECOND SHELL AND 19C REPEAT, GOING TO THE BOTTOM. THE REMAINDER OF THE 20C TERMS CAN THEN BE REMOVED THROUGH BACKSUBSTITUTION. 21C DO I = 1,NM-5,4 22C INITIAL STEP: DIAGONALIZE OMEGA(II) TERM. 23C I1 = I+1 24C I2 = I+2 25C I3 = I+3 26C I4 = I+4 27C PIVOT = A(I,5) 28 DO I = 1,NM-5 29 IF(ABS(A(I,5)).LT.TINY)THEN 30 WRITE(*,911)I,(A(I,J),J=1,10) 31 STOP 999 32 ENDIF 33 911 911 FORMAT(I5,1P10E12.3) 34 B(I) = B(I)/A(I,5) 35 DO J = 6,10 36C DO J = 5,10 37 A(I,J) = A(I,J)/A(I,5) 38 END DO 39 DO J = I+1,I+4 40 IMJ = I - J 41 IMJ5 = 5 + IMJ 42C WRITE(*,*)I,J,IMJ5 43 B(J) = B(J) - A(J,IMJ5)*B(I) 44 DO K = IMJ5+1,10+IMJ 45C WRITE(*,*)I,J,K,IMJ5,K+J-I 46 A(J,K) = A(J,K)-A(J,IMJ5)*A(I,K+J-I) 47 END DO 48 END DO 49 END DO 50 DO I = NM-4,NM-1 51 II = NM-I 52 IF(ABS(A(I,5)).LT.TINY)THEN 53 WRITE(*,911)I,(A(I,J),J=1,10) 54 STOP 999 55 ENDIF 56 B(I) = B(I)/A(I,5) 57C WRITE(*,*)I,B(I),A(I,5) 58C DO J = 5,5+II 59 DO J = 6,5+II 60 A(I,J) = A(I,J)/A(I,5) 61 END DO 62C WRITE(*,*)I,(J,A(I,J),J=6,5+II) 63 DO J = I+1,NM 64 IMJ = I - J 65 IMJ5 = 5 + IMJ 66 B(J) = B(J) - A(J,IMJ5)*B(I) 67C WRITE(*,*)I,J,B(I),B(J),IMJ5,A(J,IMJ5) 68 DO K = IMJ5+1,IMJ5+II 69 A(J,K) = A(J,K)-A(J,IMJ5)*A(I,K+J-I) 70C WRITE(*,*)I,J,K,A(J,K),IMJ5,A(J,IMJ5), 71C * K+J-I,A(I,K+J-I) 72 END DO 73 END DO 74 END DO 75 IF(ABS(A(NM,5)).LT.TINY) STOP 999 76 B(NM) = B(NM)/A(NM,5) 77C WRITE(*,*)NM,B(NM),A(NM,5) 78C WRITE(*,912)((A(I,JJ),JJ=1,10),B(I),I=1,NM) 79C 912 FORMAT(1P11E12.3) 80 81C D2W/DR2 TERM 82C PIVOT = A(I+2,3) 83C B(I1) = B(I1) - A(I1,4)*B(I) 84C DO J = 5,9 85C DO J = 4,9 86C A(I1,J) = A(I1,J)-A(I1,4)*A(I,J+1) 87C END DO 88C DW/DR (II) TERM 89C PIVOT = A(I+2,3) 90C B(I2) = B(I2)-A(I2,3)*B(I) 91C DO J = 4,8 92C DO J = 3,8 93C A(I2,J) = A(I2,J)-A(I2,3)*A(I,J+2) 94C END DO 95C OMEGA(II+1) TERM 96C PIVOT = A(I+4,1) 97C B(I4) = B(I4)-A(I4,1)*B(I) 98C DO J = 1,6 99C DO J = 2,6 100C A(I4,J)= A(I4,J)-A(I4,1)*A(I,J+4) 101C END DO 102C 3 NONZERO ENTRIES ABOVE I FOR THE NEXT TWO 103C COLUMNS 104C DO J = I1,I2 105C DIAGONALIZE ROW 106C IF(ABS(A(J,5)).LT.TINY)THEN 107C WRITE(*,911)J,(A(J,K),K=1,10) 108C STOP999 109C ENDIF 110C B(J)= B(J)/A(J,5) 111C K1 = 10 - J + I 112C DO K = 6,K1 113C DO K = 5,K1 114C A(J,K) = A(J,K)/A(J,5) 115C END DO 116C CORRECT RHS 117C DO K = J+1,J+3 118C K2 = 5 - K + J 119C B(K) = B(K) - A(K,K2)*B(J) 120C CORRECT LHS 121C K3 = K2 + 1 122C K4 = K3 + K1 - 6 123C DO JJ = K3,K4 124C DO JJ = K3-1,K4 125C A(K,JJ) = A(K,JJ)-A(K,K2)*A(J,JJ+K-J) 126C A(K,JJ) = A(K,JJ)-A(K,K2)*A(J,6+JJ-K3) 127C END DO 128C END DO 129C END DO 130C D3W/DR3 TERM 131C K1 = 7 132C CORRECT RHS 133C DO J = I3+1,I3+2 134C K2 = 5 - J + I3 135C B(J) = B(J) - A(J,K2)*B(I3) 136C CORRECT LHS 137C K3 = K2 + 1 138C DO JJ = K3,K3+1 139C DO JJ = K3-1,K3+1 140C A(J,JJ) = A(J,JJ) - A(J,K2)*A(I3,JJ+J-I3) 141C END DO 142C END DO 143C END DO 144C FINAL 2 ENTRIES ARE A BIT SIMPLER. 145C NM1 = NM - 1 146C IF(ABS(A(NM1,5)).LT.TINY)STOP999 147C B(NM1) = B(NM1)/A(NM1,5) 148C A(NM1,6) = A(NM1,6)/A(NM1,5) 149C IF(ABS(A(NM,5)).LT.TINY)STOP999 150C B(NM) = B(NM)/A(NM,5) 151C WRITE(*,912)((A(I,JJ),JJ=1,10),B(I),I=1,NM) 152C 912 FORMAT(1P11E12.3) 153C BACKSUBSTITUTION 154C DO I = NM,6,-4 155 DO I = NM,6,-1 156C REMOVE ENTRIES IN LAST COLUMN 157C DO J = I-1,I-5,-1 158C B(J) = B(J) - B(I)*A(J,5+I-J) 159 DO J = -1,-5,-1 160 B(I+J) = B(I+J) - B(I)*A(I+J,5-J) 161 END DO 162C DO J = I-2,I-5,-1 163C B(J) = B(J) - B(I-1)*A(J,4+I-J) 164C END DO 165C DO J = I-3,I-5,-1 166C B(J) = B(J) - B(I-2)*A(J,3+I-J) 167C END DO 168C DO J = I-4,I-5,-1 169C B(J) = B(J) - B(I-3)*A(J,2+I-J) 170C END DO 171C DO J = I,I-3,-1 172C DO K = J-1,I-5,-1 173C KK = 5 + I - K 174C B(K) = B(K)-B(J)*A(K,KK) 175C END DO 176C END DO 177 END DO 178 DO I = 5,2,-1 179 JJ = -(I-1) 180 DO J = -1,JJ,-1 181 B(I+J) = B(I+J) - B(I)*A(I+J,5-J) 182 END DO 183 END DO 184C B(1) = B(1) - B(2)*A(1,6) 185 RETURN 186 END

1"""Python version of bandw.f (YREC 7be4133): a banded linear solver for 2the angular-momentum transport in dadcoeft.""" 3 4import sys 5import numpy as np 6 7NMAX = 8000 8TINY = 1.0e-20 9 10 11def _print_row(a: np.ndarray, i: int) -> None: 12 """Print row i of the band, FORMAT 911: I5, 1P10E12.3. 13 14 """ 15 16 print("%5d" % i + "".join("%12.3E" % a[i, j] for j in range(1, 11))) 17 18 19def bandw(a: np.ndarray, nm: int, b: np.ndarray) -> None: 20 """Solve A x = b by Gaussian elimination without pivoting, for a band 21 matrix with 4 subdiagonals and 5 superdiagonals: a[i, 5 + k] is the 22 coefficient of x[i + k], k = -4..5 (diagonal in column 5). 23 24 Translates: bandw.f#BANDW (7be4133) 25 26 Arguments: a (NMAX+1 by 11, index 0 unused), overwritten; nm, the 27 number of unknowns; b (NMAX+1), the right-hand side, replaced in place 28 by the solution. A pivot below TINY stops the run (STOP 999). 29 30 Returns nothing. 31 32 """ 33 34 # Forward elimination, full band: rows 1..nm-5 (Fortran lines 28–49). 35 # Question (see the page): no pivoting. 36 for i in range(1, nm - 5 + 1): 37 if abs(a[i, 5]) < TINY: 38 _print_row(a, i) 39 sys.exit(999) 40 b[i] = b[i] / a[i, 5] 41 for j in range(6, 11): 42 a[i, j] = a[i, j] / a[i, 5] 43 for j in range(i + 1, i + 4 + 1): 44 imj = i - j 45 imj5 = 5 + imj # column of x_i in row j 46 b[j] = b[j] - a[j, imj5] * b[i] 47 for k in range(imj5 + 1, 10 + imj + 1): 48 a[j, k] = a[j, k] - a[j, imj5] * a[i, k + j - i] 49 50 # The last rows, whose band is cut off by the end of the matrix 51 # (Fortran lines 50–74). 52 for i in range(nm - 4, nm): 53 ii = nm - i 54 if abs(a[i, 5]) < TINY: 55 _print_row(a, i) 56 sys.exit(999) 57 b[i] = b[i] / a[i, 5] 58 for j in range(6, 5 + ii + 1): 59 a[i, j] = a[i, j] / a[i, 5] 60 for j in range(i + 1, nm + 1): 61 imj = i - j 62 imj5 = 5 + imj 63 b[j] = b[j] - a[j, imj5] * b[i] 64 for k in range(imj5 + 1, imj5 + ii + 1): 65 a[j, k] = a[j, k] - a[j, imj5] * a[i, k + j - i] 66 67 # The last row stops without printing (Fortran lines 75–76). 68 if abs(a[nm, 5]) < TINY: 69 sys.exit(999) 70 b[nm] = b[nm] / a[nm, 5] 71 72 # Back substitution (unit diagonal; Fortran lines 155–183). 73 for i in range(nm, 5, -1): 74 for j in range(-1, -6, -1): 75 b[i + j] = b[i + j] - b[i] * a[i + j, 5 - j] 76 for i in range(5, 1, -1): 77 jj = -(i - 1) 78 for j in range(-1, jj - 1, -1): 79 b[i + j] = b[i + j] - b[i] * a[i + j, 5 - j]