'************************************************************** '* SOLVE A LINEAR SYSTEM BY DIRECT FACTORIZATION * '* ---------------------------------------------------------- * '* SAMPLE RUN: * '* (Solve linear system A X = B, where: * '* * '* 1 0 0 0 0 1 1 * '* 1 1 0 0 0 -1 0 * '* A = -1 1 1 0 0 1 B = 1 * '* 1 -1 1 1 0 -1 0 * '* -1 1 -1 1 1 1 1 * '* 1 -1 1 -1 1 -1 0 ) * '* * '* LINEAR SYSTEM AX = B: * '* * '* 1.0000 0.0000 0.0000 0.0000 0.0000 1.0000 1.0000 * '* 1.0000 1.0000 0.0000 0.0000 0.0000 -1.0000 0.0000 * '* -1.0000 1.0000 1.0000 0.0000 0.0000 1.0000 1.0000 * '* 1.0000 -1.0000 1.0000 1.0000 0.0000 -1.0000 0.0000 * '* -1.0000 1.0000 -1.0000 1.0000 1.0000 1.0000 1.0000 * '* 1.0000 -1.0000 1.0000 -1.0000 1.0000 -1.0000 0.0000 * '* * '* SOLUTION: * '* .34375 * '* .3125 * '* .375 * '* .25 * '* .5 * '* .65625 * '* * '* ---------------------------------------------------------- * '* Ref.: From Numath Library By Tuan Dang Trong in Fortran 77 * '* [BIBLI 18]. * '* * '* Basic Release By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '************************************************************** 'PROGRAM TEST_DLITTL DEFDBL A-H, O-Z DEFINT I-N N = 6 DIM A(N, N), B(N), X(N), W(N, N), Z(N) FOR J = 1 TO N FOR I = 1 TO N READ A(I, J) NEXT I NEXT J DATA 1#, 1#,-1#, 1#,-1#, 1# DATA 0#, 1#, 1#,-1#, 1#,-1# DATA 0#, 0#, 1#, 1#,-1#, 1# DATA 0#, 0#, 0#, 1#, 1#,-1# DATA 0#, 0#, 0#, 0#, 1#, 1# DATA 1#,-1#, 1#,-1#, 1#,-1# FOR I = 1 TO N READ B(I) NEXT I DATA 1#,0#,1#,0#,1#,0# F$ = "###.####" F1$ = " ###.####" CLS PRINT PRINT " LINEAR SYSTEM AX=B:" PRINT FOR I = 1 TO N FOR J = 1 TO N PRINT USING F$; A(I, J); NEXT J PRINT USING F1$; B(I) NEXT I GOSUB 1000 'call DLITTL(A,B,N,X,W,Z) PRINT PRINT " SOLUTION:" FOR I = 1 TO N PRINT " "; X(I) NEXT I PRINT END 'of main program 1000 'SUBROUTINE DLITTL(A,B,N,X,W,Z) '---------------------------------------------------------------------- ' SOLVE A LINEAR SYSTEM BY DIRECT FACTORIZATION ' A*X = L*U*X = B ' L(LOWER TRIANGULAR MATRIX WITH L(I,I)=1) ' U(UPPER TRIANGULAR MATRIX),THE SUB-DIAGONAL TERMS OF L AND THE ' TERMS U(I,J) ARE STORED IN THE SAME MATRIX, W. THE RESOLUTION ' IS OBTAINED BY SOLVING TWO TRIANGULAR SYSTEMS: ' L*Z = B AND U*X = Z. ' THE PARTIAL PIVOTING IS MADE BY CHOOSING THE BIGGEST ELEMENT OF ' EACH COLUMN OF THE TRANSFORMATION MATRIX. ' INPUTS: ' A = TABLE OF SIZE (N,N) R*8 ' B = SECOND MEMBER VECTOR (N) R*8 ' N = ORDER OF LINEAR SYSTEM I*4 ' OUTPUTS: ' X = SOLUTION VECTOR (N) R*8 ' WORKING ZONE: ' W = TABLE OF SIZE (N,N) R*8 ' Z = AUXILIARY VECTOR (N) R*8 ' NOTE: ' MESSAGE '** DLITTL ** NO UNIQUE SOLUTION' IS GIVEN WHEN A ' NULL PIVOT IS FOUND. '----------------------------------------------------------------------- IP = 1 AMAX = ABS(A(1, 1)) ' SEEK MAX. ELEMENT OF FIRST COLUMN FOR J = 2 TO N IF ABS(A(J, 1)) < AMAX THEN GOTO 2 AMAX = A(J, 1) IP = J 2 NEXT J IF AMAX = 0# THEN GOTO 32 IF IP <> 1 THEN ' EXCHANGE LINES OF MATRIX A AND ORDER OF UNKNOWNS ' LINKED TO B FOR J = 1 TO N T = A(1, J) A(1, J) = A(IP, J) A(IP, J) = T NEXT J T = B(1) B(1) = B(IP) B(IP) = T END IF ' FIRST LINE OF U AND FIRST COLUMN OF L W(1, 1) = A(1, 1) FOR J = 2 TO N W(1, J) = A(1, J) W(J, 1) = A(J, 1) / W(1, 1) NEXT J ' SECOND LINE OF U AND SECOND COLUMN OF L ' (N-1)TH LINE OF U AND (N-1)TH COLUMN OF L FOR I = 2 TO N - 1 FOR J = I TO N ' SEEK PIVOT AMAX = 0# SUM = 0# FOR K = 1 TO I - 1 SUM = SUM + W(J, K) * W(K, I) NEXT K T = ABS(A(J, I) - SUM) IF T < AMAX THEN GOTO 10 AMAX = T IP = J 10 NEXT J IF AMAX = 0# THEN GOTO 32 IF IP <> I THEN ' EXCHANGE LINES OF MATRICES A , W AND ORDER OF UNKNOWNS ' LINKED TO B FOR J = 1 TO N T = A(I, J) A(I, J) = A(IP, J) A(IP, J) = T T = W(I, J) W(I, J) = W(IP, J) W(IP, J) = T NEXT J T = B(I) B(I) = B(IP) B(IP) = T END IF ' CALCULATE THE U(I,I) SUM = 0# FOR K = 1 TO I - 1 SUM = SUM + W(I, K) * W(K, I) NEXT K W(I, I) = A(I, I) - SUM IF W(I, I) = 0# THEN GOTO 32 FOR J = I + 1 TO N SUM = 0# FOR K = 1 TO I - 1 SUM = SUM + W(I, K) * W(K, J) NEXT K ' CALCULATE THE U(I,J) W(I, J) = A(I, J) - SUM SUM = 0# FOR K = 1 TO I - 1 SUM = SUM + W(J, K) * W(K, I) NEXT K ' CALCULATE THE L(I,J) W(J, I) = (A(J, I) - SUM) / W(I, I) NEXT J NEXT I SUM = 0# FOR K = 1 TO N - 1 SUM = SUM + W(N, K) * W(K, N) NEXT K ' CALCULATE U(N,N) W(N, N) = A(N, N) - SUM IF W(N, N) = 0# THEN GOTO 32 ' SOLVE SYSTEM Z = L(-1)*B Z(1) = B(1) FOR I = 2 TO N SUM = 0# FOR K = 1 TO I - 1 SUM = SUM + W(I, K) * Z(K) NEXT K Z(I) = B(I) - SUM NEXT I ' SOLVE SYSTEM X = U(-1)*Z X(N) = Z(N) / W(N, N) FOR I = N - 1 TO 1 STEP -1 SUM = 0# FOR K = I + 1 TO N SUM = SUM + W(I, K) * X(K) NEXT K X(I) = (Z(I) - SUM) / W(I, I) NEXT I ' END OF RESOLUTION OF A*X = L*(U*X) = B RETURN 32 PRINT ' ** DLITTL ** NO UNIQUE SOLUTION.' RETURN ' End of file tdlittl.bas