'===========================================================================
' Subject: SOLVE LINEAR SYSTEM USING GAUSS   Date: 09-23-99 (08:21)       
' Author:  Robert L. Roach                   Code: QB, QBasic, PDS        
' Origin:  rroach@cvd-pro.com              Packet: ALGOR.ABC
'===========================================================================
'Below is an example.  It is a simple algorithm for
'solving a linear system by Gauss elimination with no pivoting (since
'most of my linear systems are, fortunately, diagonally dominant). The
'actual matrix solver is in the subroutine, the main program is thus
'not necessary and only is given as an example of the use of the matrix
'solver.

REM---------------------- PROGRAM MSOLV ----------------------------
REM
REM    This program demonstrates the solution of a full linear
REM    system using Gauss elimination with no pivoting.  Please
REM    note the order of the indices of the augmented matrix.
REM    col, row format is used as I found it faster.  Thus the system
REM    is written as:
REM
REM          a11*x1 + a21*x2 + a31*x3 + ... = a(n+1,1)
REM          a12*x1 + a22*x2 + a32*x3 + ... = a(n+1,2)
REM                             etc.
REM          a1n*x1 + a2n*x2 + a3n*x3 + ... = a(n+1,n)
REM
REM    and the a(n+1,n) is the augmented matrix.
REM
REM    The main program only loads up the matrix.  The real
REM    work of solving the thing is in the subroutine 2000.  Note
REM    that the solution is returned in the column a(n+1,i).
REM
REM-----------------------------------------------------------------------------------
REM---------- DEFAULT NAMES (FORTRAN CONVENTION)
        DEFDBL A-H, O-Z
        DEFINT I-N

REM----------- MAXIMUM ARRAY SIZE -----------
        NMAX = 51
        DIM AA(NMAX + 1, NMAX)

REM------------ USEFUL STUFF -------------
        PI = 3.141592654

REM----------- LOAD MATRIX
        N = 5

        FOR J = 1 TO N
          FOR I = 1 TO N
            AA(I, J) = 1 / (1 + ABS(I - J))
          NEXT
          AA(N + 1, J) = SIN(2 * PI * (J - 1) / (N - 1))
        NEXT
       
       GOSUB 2000

REM------------ PRINT SOLUTION FOR THE FIRST 4 UNKNOWNS
        PRINT USING "   AA(5,1) = #.####^^^^"; AA(5, 1)
        PRINT USING "   AA(5,2) = #.####^^^^"; AA(5, 2)
        PRINT USING "   AA(5,3) = #.####^^^^"; AA(5, 3)
        PRINT USING "   AA(5,4) = #.####^^^^"; AA(5, 4)

REM------------------------------------------------------------
        END



2000 REM-------------------------------------------------------
REM        SUBROUTINE MSOLV(A, N)
REM------------------------------------------------------------
REM----------- ELIMINATE LOWER TRIANGULAR ELEMENTS
        FOR I = 1 TO N - 1
          FOR J = I + 1 TO N
            R = -AA(I, J) / AA(I, I)
            FOR K = I + 1 TO N + 1
              AA(K, J) = AA(K, J) + R * AA(K, I)
            NEXT
          NEXT
        NEXT
     
REM----------- BACK SUBSTITUTION
        AA(N + 1, N) = AA(N + 1, N) / AA(N, N)
       
        FOR J = N - 1 TO 1 STEP -1
          R = AA(N + 1, J)
          FOR I = J + 1 TO N
            R = R - AA(I, J) * AA(N + 1, I)
          NEXT
          AA(N + 1, J) = R / AA(J, J)
        NEXT

REM------------------------------------------------------------
        RETURN