'**************************************************************************** '* Solve a symmetric linear system using subroutine SYMSOL * '* ------------------------------------------------------------------------ * '* SAMPLE RUN: * '* ( Solve symmetric linear system: * '* X1 -2X2 +3X3 - 4X4 + 5X5 - 6X6 = -3 * '* -2X1 +3X2 -5X3 + 7X4 - 4X5 + 8X6 = 7 * '* 3X1 -5X2 +5X3 + 6X4 - 2X5 + X6 = 8 * '* -4X1 +7X2 +6X3 - 9X4 + 20X5 + 0.5X6 = 20.5 * '* 5X1 -4X2 -2X3 + 20X4 + 11X5 +0.3333X6 = 30.3333 * '* -6X1 +8X2 + X3 +0.5X4 +0.3333X5 + 13X6 = 16.8333 ) * '* * '* SYMMETRIC SYSTEM TO SOLVE: * '* 1.000000 -2.000000 3.000000 -4.000000 5.000000 -6.000000 -3.000000 * '* -2.000000 3.000000 -5.000000 7.000000 -4.000000 8.000000 7.000000 * '* 3.000000 -5.000000 5.000000 6.000000 -2.000000 1.000000 8.000000 * '* -4.000000 7.000000 6.000000 -9.000000 20.000000 0.500000 20.500000 * '* 5.000000 -4.000000 -2.000000 20.000000 11.000000 0.333300 30.333300 * '* -6.000000 8.000000 1.000000 0.500000 0.333300 13.000000 16.833300 * '* * '* SOLUTIONS: * '* 1.00000008 * '* 0.99999954 * '* 0.99999976 * '* 0.99999977 * '* 1.00000013 * '* 1.00000031 * '* * '* ------------------------------------------------------------------------ * '* Reference: From Numath Library By Tuan Dang Trong in Fortran 77 * '* [BIBLI 18]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '**************************************************************************** 'PROGRAM TEST_SYMSOL DEFDBL A-H, O-Z DEFINT I-N N = 6 DIM A(N, N), B(N), X(N), D(N) 'define upper half part of matrix A A(1, 1) = 1#: A(1, 2) = -2#: A(1, 3) = 3#: A(1, 4) = -4#: A(1, 5) = 5#: A(1, 6) = -6# A(2, 2) = 3#: A(2, 3) = -5#: A(2, 4) = 7#: A(2, 5) = -4#: A(2, 6) = 8# A(3, 3) = 5#: A(3, 4) = 6#: A(3, 5) = -2#: A(3, 6) = 1# A(4, 4) = -9#: A(4, 5) = 20#: A(4, 6) = .5 A(5, 5) = 11#: A(5, 6) = .3333 A(6, 6) = 13# 'define 2nd member B(1) = -3#: B(2) = 7#: B(3) = 8#: B(4) = 20.5: B(5) = 30.3333: B(6) = 16.8333 F\$ = "####.######" 'optional to print full data FOR I = 2 TO N FOR J = 1 TO I - 1 A(I, J) = A(J, I) NEXT J NEXT I CLS PRINT PRINT " SYMMETRIC SYSTEM TO SOLVE:" FOR I = 1 TO N FOR J = 1 TO N PRINT USING F\$; A(I, J); NEXT J PRINT " "; PRINT USING F\$; B(I) NEXT I GOSUB 1000 'call SYMSOL(N,A,B,X,D) G\$ = "####.########" 'print solutions PRINT " SOLUTIONS:" FOR I = 1 TO N PRINT USING G\$; X(I) NEXT I END 1000 'SUBROUTINE SYMSOL (N,A,B,X,D) '----------------------------------------------------------------------- ' SOLVE A SYMMETRIC LINEAR SYSTEM A.X = B ' A IS SUPPOSED WELL BALANCED AND CAN BE NOT POSITIVE DEFINITE, ' THE COEFFICIENTS OF A ARE STORED IN THE UPPER TRIANGULAR HALF ' OF THE MATRIX AND ARE LOST DURING THE PROCESS. ' THE USED METHOD IS BASED UPON THE DECOMPOSITION OF THE INITIAL ' MATRIX INTO A PRODUCT OF 3 MATRICES: A = (U TRANSPOSE).D.U ' INPUTS: ' N NUMBER OF LINES OF THE LINEAR SYSTEM ' A TABLE OF DIMENSION N X N STORING THE COEFFICIENTS ' OF THE SYSTEM ' B VECTOR OF DIMENSION N (SECOND MEMBER) ' OUTPUT: ' X SYSTEM SOLUTION ' WORKING ZONE: ' D DIAGONAL D OF DECOMPOSITION '----------------------------------------------------------------------- FOR I = 1 TO N L = I - 1 M = I + 1 S = A(I, I) IF I = 1 THEN GOTO 2 FOR K = 1 TO L S = S - A(K, I) * D(K) * A(K, I) NEXT K 2 D(I) = S IF M > N THEN GOTO 6 FOR J = M TO N S = A(I, J) IF I = 1 THEN GOTO 4 FOR K = 1 TO L S = S - A(K, I) * D(K) * A(K, J) NEXT K 4 A(I, J) = S / D(I) NEXT J NEXT I 6 FOR I = 1 TO N S = B(I) L = I - 1 IF I = 1 THEN GOTO 8 FOR K = 1 TO L S = S - A(K, I) * X(K) NEXT K 8 X(I) = S NEXT I FOR L = 1 TO N I = N - L + 1 S = X(I) / D(I) M = I + 1 IF M > N THEN GOTO 10 FOR K = M TO N S = S - A(I, K) * X(K) NEXT K 10 X(I) = S NEXT L RETURN 'end of file tsymsol.bas