'************************************************************** '* 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