'**************************************************************
'* Solve banded linear system AX = B By LU decomposition *
'* ---------------------------------------------------------- *
'* SAMPLE RUN: *
'* (Solve banded linear system AX = B, where: *
'* 1 10 0 0 0 1 *
'* 9 2 20 0 0 1 *
'* A = 0 19 3 30 0 B = 1 *
'* 0 0 29 4 40 1 *
'* 0 0 0 39 5 1 ) *
'* *
'* Note: A is given as follows: *
'* 0 0 0 0 0 *
'* 0 10 20 30 40 *
'* 1 2 3 4 5 *
'* 9 19 29 39 0 *
'* 0 0 0 0 0 *
'* *
'* SOLUTION: *
'* *
'* 0.40664649 *
'* 0.05933535 *
'* -0.13892446 *
'* 0.00964672 *
'* 0.12475556 *
'* *
'* Error code: 0 *
'* *
'* ---------------------------------------------------------- *
'* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77 *
'* [BIBLI 18]. *
'* *
'* Basic Release 1.0 By J-P Moreau, Paris. *
'* (www.jpmoreau.fr) *
'**************************************************************
'PROGRAM TEST_NSBSLV
DEFDBL A-H, O-Z
DEFINT I-N
N = 5 'size of linear system
DIM A(N, N), B(N), X(N), IPVT(N), TMP(N), TMP1(N)
FOR I = 1 TO N
FOR J = 1 TO N
A(I, J) = 0#
NEXT J
NEXT I
MU = 1 'one upper subdiagonal
ML = 1 'one lower subdiagonal
' Define banded A matrix
A(2, 2) = 10#: A(2, 3) = 20#: A(2, 4) = 30#: A(2, 5) = 40#
A(3, 1) = 1#: A(3, 2) = 2#: A(3, 3) = 3#: A(3, 4) = 4#: A(3, 5) = 5#
A(4, 1) = 9#: A(4, 2) = 19#: A(4, 3) = 29#: A(4, 4) = 39#
' Example #1 (B(1)..B(5) = ONE)
FOR I = 1 TO N
B(I) = 1#
NEXT I
' Example #2 (results are X1..X5 = ONE)
' B(1)=11#: B(2)=31#: B(3)=52#: B(4)=73#: B(5)=44#
' call LU factorization
GOSUB 1000 'call NSBFAC(A,N,ML,MU,IPVT,IND)
' call appropriate solver
GOSUB 2000 'call NSBSLV(A,N,ML,MU,IPVT,B,X)
' print solution and error code (must be zero)
F$ = "###.########"
CLS
PRINT
PRINT " SOLUTION:"
PRINT
FOR I = 1 TO N
PRINT USING F$; X(I)
NEXT I
PRINT
PRINT " Error code: "; IND
PRINT
END 'of main program
500 'Function IAMAX(TMP,N1)
T = 0#
FOR I1 = 1 TO N1
IF ABS(TMP(I1)) > T THEN
T = ABS(TMP(I1))
IAMAX = I1
END IF
NEXT I1
RETURN
600 'Subroutine SCALE(N1,T,TMP)
FOR I1 = 1 TO N1
TMP(I1) = T * TMP(I1)
NEXT I1
RETURN
700 'Subroutine DAXPY(N1,AA,TMP,TMP1)
FOR I1 = 1 TO N1
TMP1(I1) = TMP1(I1) + AA * TMP(I1)
NEXT I1
RETURN
1000 'Subroutine NSBFAC(B,N,ML,MU,IPVT,IND)
'-------------------------------------------------------------------
' LU factorization of a band matrix (non symmetric) with partial
' pivoting.
' INPUTS:
' B : banded matrix. The correspondance between full matrix
' A(i,j) and band matrix B(k,l) is given by following sequence:
' m=ml+mu+1
' for j=1 to n
' i1=max(1,j-mu)
' i2=min(n,j+ml)
' for i=i1 to i2
' k=i-j+m;
' B(k,j)=A(i,j)
' next i
' next j
' N : size of B
' ML : number of lower diagonals
' MU : number of upper diagonals
' OUTPUTS:
' B : banded matrix storing the LU elements, the ML first lines
' of which are used for pivoting.
' IPVT: integer vector of size N storing the pivoting indices.
' IND : flag = 0, B is non singular
' = k, B may be singular
'-------------------------------------------------------------------
'Note: here matrix B is called A (shared with calling program).
'Labels: 10, 20
M = ML + MU + 1
IND = 0
J0 = MU + 2
'J1=MIN(N,M)-1
J1 = M - 1
IF J1 >= J0 THEN
FOR JZ = J0 TO J1
I0 = M + 1 - JZ
FOR I = I0 TO ML
A(I, JZ) = 0#
NEXT I
NEXT JZ
END IF
JZ = J1
JU = 0
NM1 = N - 1
IF NM1 >= 1 THEN
FOR K = 1 TO NM1
KP1 = K + 1
JZ = JZ + 1
IF JZ <= N AND ML >= 1 THEN
FOR I = 1 TO ML
A(I, JZ) = 0#
NEXT I
END IF
'LM=MIN(ML,N-K)
IF ML <= N - K THEN
LM = ML
ELSE
LM = N - K
END IF
FOR I = 1 TO LM + 1
TMP(I) = A(I + M - 1, K)
NEXT I
'L=IAMAX(TMP,LM+1)+M-1
N1 = LM + 1: GOSUB 500: L = IAMAX + M - 1
IPVT(K) = L + K - M
IF A(L, K) = 0# THEN GOTO 10
IF L <> M THEN
T = A(L, K)
A(L, K) = A(M, K)
A(M, K) = T
END IF
T = -1# / A(M, K)
FOR I = 1 TO LM
TMP(I) = A(I + M, K)
NEXT I
'SCALE(LM,T,TMP)
N1 = LM: GOSUB 600
FOR I = 1 TO LM
A(I + M, K) = TMP(I)
NEXT I
'JU=MIN(MAX(JU,MU+IPVT(K)),N)
IF JU >= MU + IPVT(K) THEN
MAX = JU
ELSE
MAX = MU + IPVT(K)
END IF
IF MAX <= N THEN
JU = MAX
ELSE
JU = N
END IF
MM = M
IF JU >= KP1 THEN
FOR J = KP1 TO JU
L = L - 1
MM = MM - 1
T = A(L, J)
IF L <> MM THEN
A(L, J) = A(MM, J)
A(MM, J) = T
END IF
FOR I = 1 TO LM
TMP(I) = A(I + M, K)
TMP1(I) = A(I + MM, J)
NEXT I
'DAXPY(LM,T,TMP,TMP1)
N1 = LM: AA = T: GOSUB 700
FOR I = 1 TO LM
A(I + MM, J) = TMP1(I)
NEXT I
NEXT J
END IF
GOTO 20
10 IND = K
20 NEXT K
END IF
IPVT(N) = N
IF A(M, N) = 0# THEN IND = N
RETURN
2000 'Subroutine NSBSLV(A,N,ML,MU,IPVT,B,X)
'--------------------------------------------------------------------
' Solve banded linear system AX = B
' INPUTS:
' A : banded matrix as output of LU factorization by NSBFAC
' (see storing mode in NSBFAC subroutine).
' N : order of A ---------------
' ML : number of lower diagonals
' MU : number of upper diagonals
' IPVT: integer vector of size N storing the pivoting indices
' as output of NSBFAC.
' B : second member vector
' OUTPUT:
' X : solution vector
'--------------------------------------------------------------------
FOR I = 1 TO N
X(I) = B(I)
NEXT I
M = ML + MU + 1
NM1 = N - 1
' solve L*y = B
IF ML <> 0 AND NM1 >= 1 THEN
FOR K = 1 TO NM1
'LM=MIN(ML,N-K)
IF ML <= N - K THEN
LM = ML
ELSE
LM = N - K
END IF
L = IPVT(K)
T = X(L)
IF L <> K THEN
X(L) = X(K)
X(K) = T
END IF
FOR I = 1 TO LM
TMP(I) = A(I + M, K)
TMP1(I) = X(I + K)
NEXT I
'DAXPY(LM,T,TMP,TMP1)
N1 = LM: AA = T: GOSUB 700
FOR I = 1 TO LM
X(I + K) = TMP1(I)
NEXT I
NEXT K
END IF
' solve U*y = X
FOR K = N TO 1 STEP -1
X(K) = X(K) / A(M, K)
'LM=MIN(K,M)-1
IF K <= M THEN
LM = K - 1
ELSE
LM = M - 1
END IF
LA = M - LM
LB = K - LM
T = -X(K)
FOR I = 1 TO LM
TMP(I) = A(I + LA - 1, K)
TMP1(I) = X(I + LB - 1)
NEXT I
'DAXPY(LM,T,TMP,TMP1)
N1 = LM: AA = T: GOSUB 700
FOR I = 1 TO LM
X(I + LB - 1) = TMP1(I)
NEXT I
NEXT K
RETURN
'end of file nsbslv.bas