From: marwk@levels.unisa.edu.au
Newsgroups: unisa.cis.prog.f.02462
Subject: PRACTICAL 8: Solution
Date: 9 Nov 92 11:45:43 +1030
Organization: University of South Australia

      PROGRAM LINSYS
************************************************************************
* Program to solve a linear system  using Gaussian elimination.        *
* Identifiers used are:                                                *
*     AUG    : augmented matrix for the linear system                  *
*     I, J   : indices                                                 *
*     LIMAUG : parameter (LIMIT + 1) for maximum # columns in AUG      *
*     LIMIT  : parameter giving maximum dimensions of matrix           *
*     N      : number of equations and unknowns                        *
*     SINGUL : indicates if matrix is (nearly) singular                *
*     X      : solution vector                                         *
************************************************************************

      INTEGER LIMIT, LIMAUG
      PARAMETER (LIMIT = 10, LIMAUG = 11)
      REAL AUG(LIMIT,LIMAUG), X(LIMIT)
      INTEGER N, I, J
      LOGICAL SINGUL

* Read coefficient matrix and constant vector

      PRINT *, 'ENTER NUMBER OF EQUATIONS'
      READ *, N
      PRINT *, 'ENTER COEFFICIENT MATRIX ROWWISE'
      READ *, ((AUG(I,J), J = 1, N), I = 1, N)
      PRINT *, 'ENTER CONSTANT VECTOR'
      READ *, (AUG(I, N + 1), I = 1, N)

* Display the initial tableau

      CALL DISPLAY_AUG(AUG, LIMIT, N, 0)

* Use subroutine GAUSS to find the solution A,
* and then display X.

      CALL GAUSS(AUG, LIMIT, LIMAUG, N, X, SINGUL)
      IF (.NOT. SINGUL) THEN
          PRINT *, 'SOLUTION VECTOR IS'
          DO 10 I = 1, N
              PRINT 100, I, X(I)
100           FORMAT(1X, 'X(', I2, ') =', F8.3)
10        CONTINUE
      ELSE
          PRINT *, 'MATRIX IS (NEARLY) SINGULAR'
      END IF
      END

**GAUSS*****************************************************************
* Subroutine to find solution of linear system AX = B using Gaussian   *
* elimination, provided a unique solution.  If the matrix is (nearly)  *
* singular, SINGUL is returned as true, and the solution X is          *
* is undefined.   Local identifiers are:                         *
*     I,J,K  : indices                                                 *
*     MULT   : multiplier used to eliminate an unknown                 *
*     ASBPIV : absolute value of pivot element                         *
*     PIVROW : row containing pivot element                            *
*     EPSIL  : a small positive real value ("almost zero")             *
*     TEMP   : used to interchange rows of matrix                      *
************************************************************************

      SUBROUTINE GAUSS(AUG, LIM, LIMAUG, N, X, SINGUL)

      INTEGER LIM, LIMAUG
      REAL AUG(LIM, LIMAUG), X(10), TEMP, MULT, EPSIL
      PARAMETER (EPSIL = 1E-6)
      INTEGER N, PIVROW
      LOGICAL SINGUL

      SINGUL = .FALSE.
      DO 50 I = 1, N

*         Locate pivot element

          ABSPIV = ABS(AUG(I,I))
          PIVROW = I
          DO 10 K = I + 1, N
              IF (ABS(AUG(K,I)) .GT. ABSPIV) THEN
                  ABSPIV = ABS(AUG(K,I))
                  PIVROW = K
              END IF
10        CONTINUE

*         Check if matrix is (nearly) singular

          IF (ABSPIV .LT. EPSIL) THEN
              SINGUL = .TRUE.
              RETURN
          END IF

*         It isn't, so interchange rows PIVROW and I if necessary

          IF (PIVROW .NE. I) THEN
              DO 20 J = 1, N + 1
                  TEMP = AUG(I,J)
                  AUG(I,J) = AUG(PIVROW,J)
                  AUG(PIVROW,J) = TEMP
20            CONTINUE
          END IF

*         Eliminate Ith unknown from equations I + 1, ..., N

          DO 40 J = I + 1, N
              MULT = -AUG(J,I) / AUG(I,I)
              DO 30 K = I, N + 1
                  AUG(J,K) = AUG(J,K) +  MULT * AUG(I,K)
30            CONTINUE
40        CONTINUE

          CALL DISPLAY_AUG(AUG, LIMIT, N, I)

50    CONTINUE

* Find the solutions by back substitution

      X(N) = AUG(N, N + 1) / AUG(N,N)
      DO 70 J = N - 1, 1, -1
          X(J) = AUG(J, N + 1)
          DO 60 K = J + 1, N
              X(J) = X(J) - AUG(J,K) * X(K)
60        CONTINUE
          X(J) = X(J) / AUG(J,J)
70    CONTINUE

      END


      SUBROUTINE DISPLAY_AUG(AUG, LIMIT, N, ITER)
      INTEGER LIMIT, N, ITER
      REAL AUG(10, 11)

      PRINT 100, ITER
  100 FORMAT ('0', 'Iteration: ', I2)

      PRINT 200, ((AUG(I,J), J = 1, N+1), I = 1, N)
  200 FORMAT (' ', <N>F6.2, 4X, F6.2)

      END
--
Raymond Kennington
LECTURER
University of South Australia   | Act in haste and repent at leisure
P.O. Box 1                      | Code too soon and debug forever
Ingle Farm                      | Knobs, knobs everywhere,
South Australia                 |              just vary a knob to think!
#! rnews 2269