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