'********************************************************* '* Solving a linear system AX = B by LU decomposition * '* * '* Basic Version By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ----------------------------------------------------- * '* Reference: * '* * '* "Numerical Recipes By W.H. Press, B. P. Flannery, * '* S.A. Teukolsky and W.T. Vetterling, Cambridge * '* University Press, 1986" [BIBLI 08]. * '* ----------------------------------------------------- * '* SAMPLE RUN: * '* * '* Input file (test_lu.dat): * '* * '* 4 * '* 8 2 3 12 25.0 * '* 2 4 7 0.25 13.25 * '* 3 7 3 5 18.0 * '* 12 0.25 5 2 19.25 * '* * '* Output file (test_lu.lst): * '* * '* --------------------------------------------------- * '* LINEAR SYSTEM TO BE SOLVED: * '* --------------------------------------------------- * '* N=4 * '* * '* 8.000000 2.000000 3.000000 12.00000 25.00000 * '* 2.000000 4.000000 7.000000 0.250000 13.25000 * '* 3.000000 7.000000 3.000000 5.000000 18.00000 * '* 12.00000 0.250000 5.000000 2.000000 19.25000 * '* * '* System solution: * '* * '* X1= 1.000000 * '* X2= 1.000000 * '* X3= 1.000000 * '* X4= 1.000000 * '* --------------------------------------------------- * '* * '* Subroutines used: * '* 1000 LUDCMP: LU decomposition of a real square matrix * '* 2000 LUBKSB: Solve linear system AX = B, A being the * '* output of LUDCMP. * '* * '********************************************************* DEFINT I-N DEFDBL A-H, O-Z CLS : PRINT INPUT " Data file name (without .dat): ", nom\$ in\$ = nom\$ + ".dat" ou\$ = nom\$ + ".lst" OPEN in\$ FOR INPUT AS #1 OPEN ou\$ FOR OUTPUT AS #2 INPUT #1, N 'size of linear system n2 = (N + 1) * (N + 1) DIM A(n2), B(N + 1), temp(N + 2), INDX(N + 1) '------------------------------------------------------------------ 'NOTA: Matrix A of size N x N is stored in a vector, line by line. ' The element (i,j) of the matrix is at location i(n+1) + j in ' the vector. The access is faster. '------------------------------------------------------------------ PRINT #2, "--------------------------------------------------" PRINT #2, " LINEAR SYSTEM TO BE SOLVED:" PRINT #2, "--------------------------------------------------" PRINT #2, PRINT #2, " N="; N PRINT #2, FOR I = 1 TO N FOR J = 1 TO N + 1 INPUT #1, temp(J) 'read a line NEXT J FOR J = 1 TO N A(I * (N + 1) + J) = temp(J) NEXT J B(I) = temp(N + 1) FOR J = 1 TO N + 1 PRINT #2, USING "###.######"; temp(J); NEXT J PRINT #2, NEXT I CLOSE #1 GOSUB 1000 'LUDCMP(A,n,INDX,D,rc); IF rc = 0 THEN GOSUB 2000 'LUBKSB(A,n,INDX,B); IF rc = 1 THEN PRINT #2, " The system matrix is singular !" ELSE PRINT #2, PRINT #2, " System solution:" PRINT #2, FOR I = 1 TO N PRINT #2, USING " X## = ##.######"; I; B(I) PRINT #2, NEXT I END IF PRINT #2, "--------------------------------------------------" PRINT #2, " End of file "; ou\$ CLOSE #2 PRINT PRINT " Results in file "; ou\$ PRINT END 'of main program '*************************************************************** '* Given an N x N matrix A, this routine replaces it by the LU * '* decomposition of a rowwise permutation of itself. A and N * '* are input. INDX is an output vector which records the row * '* permutation effected by the partial pivoting; D is output * '* as -1 or 1, depending on whether the number of row inter- * '* changes was even or odd, respectively. This routine is used * '* in combination with LUBKSB to solve linear equations or to * '* invert a matrix. Return code is 1, if matrix is singular. * '*************************************************************** 1000 'SUBROUTINE LUDCMP NMAX = 100: TINY = 1.5E-16 DIM VV(N) D = 1: CODE = 0: N1 = N + 1 FOR I = 1 TO N AMAX = 0# FOR J = 1 TO N IF ABS(A(I * N1 + J)) > AMAX THEN AMAX = ABS(A(I * N1 + J)) NEXT J IF (AMAX < TINY) THEN CODE = 1 RETURN END IF VV(I) = 1# / AMAX NEXT I 'i loop FOR J = 1 TO N FOR I = 1 TO J - 1 SUM = A(I * N1 + J) FOR K = 1 TO I - 1 SUM = SUM - A(I * N1 + K) * A(K * N1 + J) NEXT K A(I * N1 + J) = SUM NEXT I 'i loop AMAX = 0# FOR I = J TO N SUM = A(I * N1 + J) FOR K = 1 TO J - 1 SUM = SUM - A(I * N1 + K) * A(K * N1 + J) NEXT K A(I * N1 + J) = SUM DUM = VV(I) * ABS(SUM) IF DUM >= AMAX THEN IMAX = I AMAX = DUM END IF NEXT I 'i loop IF J <> IMAX THEN FOR K = 1 TO N DUM = A(IMAX * N1 + K) A(IMAX * N1 + K) = A(J * N1 + K) A(J * N1 + K) = DUM NEXT K 'k loop D = -D VV(IMAX) = VV(J) END IF INDX(J) = IMAX IF ABS(A(J * N1 + J)) < TINY THEN A(J * N1 + J) = TINY END IF IF J <> N THEN DUM = 1# / A(J * N1 + J) FOR I = J + 1 TO N A(I * N1 + J) = A(I * N1 + J) * DUM NEXT I END IF NEXT J 'j loop RETURN 'subroutine LUDCMP '****************************************************************** '* Solves the set of N linear equations A . X := B. Here A is * '* input, not as the matrix A but rather as its LU decomposition, * '* determined by the routine LUDCMP. INDX is input as the permuta-* '* tion vector returned by LUDCMP. B is input as the right-hand * '* side vector B, and returns with the solution vector X. A, N and* '* INDX are not modified by this routine and can be used for suc- * '* cessive calls with different right-hand sides. This routine is * '* also efficient for plain matrix inversion. * '****************************************************************** 2000 'SUBROUTINE LUBKSB II = 0: N1 = N + 1 FOR I = 1 TO N LL = INDX(I) SUM = B(LL) B(LL) = B(I) IF II <> 0 THEN FOR J = II TO I - 1 SUM = SUM - A(I * N1 + J) * B(J) NEXT J ELSEIF SUM <> 0 THEN II = I END IF B(I) = SUM NEXT I 'i loop FOR I = N TO 1 STEP -1 SUM = B(I) IF I < N THEN FOR J = I + 1 TO N SUM = SUM - A(I * N1 + J) * B(J) NEXT J END IF B(I) = SUM / A(I * N1 + I) NEXT I 'i loop RETURN 'End of file lu.bas