'********************************************************
'*  Solving a symmetric linear system by Gauss method   *
'*                                                      *
'*                 Basic version by J-P Moreau, Paris   *
'*                          (www.jpmoreau.fr)           *
'* ---------------------------------------------------- *
'* Reference:                                           *
'*                                                      *
'*   "ALGEBRE Algorithmes et programmes en Pascal       *
'*    de Jean-Louis Jardrin - Dunod BO-PRE 1988"        *
'*    [BIBLI 10].                                       *
'* ---------------------------------------------------- *
'* SAMPLE RUN:                                          *
'*                                                      *
'* Input file (matsym.dat):                             *
'*                                                      *
'* 4                                                    *
'* 8  2  3  12     25                                   *
'*    4  7   0.25  13.25                                *
'*       3   5     18                                   *
'*           2     19.25                                *
'*                                                      *
'* Output file (matsym.lst):                            *
'*                                                      *
'* Size of symmetric system: 4                          *
'*                                                      *
'* SYSTEM TO BE SOLVED:                                 *
'*                                                      *
'*  8.0000  2.0000  3.0000 12.0000 25.0000              *
'*  2.0000  4.0000  7.0000  0.2500 13.2500              *
'*  3.0000  7.0000  3.0000  5.0000 18.0000              *
'* 12.0000  0.2500  5.0000  2.0000 19.2500              *
'*                                                      *
'* System solution:                                     *
'*                                                      *
'*   x1 =   1.000000                                    *
'*   x2 =   1.000000                                    *
'*   x3 =   1.000000                                    *
'*   x4 =   1.000000                                    *
'*                                                      *
'********************************************************
defint i-n

NMAX  = 30
NMAXP = 31

dim A(NMAX,NMAXP), X(NMAX)

  CLS

  open "matsym.dat" for input as #1
  open "matsym.lst" for output as #2

  input #1, n

  print #2,
  print #2, "  Size of symmetric system: "; n

  print #2,
  print #2, "  SYSTEM TO BE SOLVED:"
  print #2,

  'read symmetric system (upper triangle and right hand terms)
  for i=1 to n
    for j=i to n+1
      input #1, A(i,j)
      A(j,i)=A(i,j)
    next j
  next i

  Close #1

  'write to output file symmetric system (complete)
  for i=1 to n
    for j=1 to n+1
      print #2, using " ##.######"; A(i,j);
    next j
    print #2,
  next i

  eps=1e-10  'desired precision

  ' call routine RSLSG(eps,n,A,it,X)
  gosub 1000

  if it=0 then
    print #2, "  *** ERROR ***"
  else
    print #2,
    print #2, "  System solution:"
    print #2,
    for i=1 to n
      print #2,	using "    x(##) = ###.######"; i; X(i)
    next i
  end if
  close #2
  print
  print " Results in file matsym.lst."

END 'of main program


'***********************************************************************
'*   This procedure solves a symmetric linear system by Gauss method   *
'*                     in simple precision                             *
'* ------------------------------------------------------------------- *
'* The input matrix includes the right hand column, only the upper     *
'* triangle terms are given.                                           *
'* INPUTS:                                                             *
'*           eps:    desired precision (real)                          *
'*             n:    size of system (integer)                          *
'*             A:    extended matrix of system of type MAT             *
'* OUTPUTS:                                                            *
'*            it:    error indicator (0=singular system, 1=OK)         *
'*             X:    system solution (x1,...,xn) of type VEC           *
'***********************************************************************
1000 'Routine RSLSG(eps:real; n:integer; A:MAT; VAR it:integer; VAR X:VEC)
  it=1 : k=1

1100 if abs(A(k,k)) < eps then
    it=0
  else
    for i=k+1 to n
      q0=A(k,i)/A(k,k)
      for j=i to n+1
        A(i,j)=A(i,j)-q0*A(k,j)
      next j
    next i
    k=k+1
  end if
  if it<>0 and k<n then goto 1100

  if it=1 AND abs(A(n,n))>=eps then
    X(n)=A(n,n+1)/A(n,n)
    for i=n-1 to 1 step -1
      s=0#
      for j=i+1 to n
        s=s+A(i,j)*X(j)
      next j
      X(i)=(A(i,n+1)-s)/A(i,i)
    next i
  else
    it=0
  end if
return

'End of file syslin.bas