How to Solve a System of Linear Equations of n Variables using the Gauss Elimination Method

Gaussian Elimination is an algorithm to solve a system of linear equations of n variables. It uses elementary operations to find the solution to the system of linear equations. It is an efficient method to find the determinant and inverse of a matrix as well. In this program, we will learn how to solve a system of equations using this method. We will also learn how to calculate the determinant of the n x n matrix used in the process.

• We create a module for 3 matrices a(n,m) for the augmented matrix, b(n) for temporary values of x and x(n) for the final solutions. We will need this in the main program as well as the subroutine.
• We read the number of variables, n at the time of execution.
• We read the augmented matrix from an input file and display it on the screen.
• Optional: At this stage, we can create another matrix check(n,m) to store the initial values of elements of a(n,m). This can help us substitute the obtained values of x1,x2,...,xn later and verify them.
• We call a subroutine and display the results only in the square matrix a(n,n) is non-singular. This is also called the coefficient matrix.The subroutine needs only the input m or n. The rest of the values are held in the module. The subroutine generates the arrays and the output sing to tell the status of singularity of the matrix a(n,n).
• Optional: We can evaluate the L.H.S. and R.H.S by using two simple one-column arrays of n rows. R.H.S will simply be the initial  values of a(:,m) and L.H.S. will be a1*x1+a12*x2+...a1n*xn, a21*x1+a22*x2+...a2n*xn,and  so on for all n equations. If L.H.S. and R.H.S. are equal for all rows, we have successfully recreated the equation to verify the results.

Now for the subroutine:

• The first step is forward elimination. Starting with the first row, we make sure that the diagonal element has the largest element for its column in that row and all the subsequent rows. If not, we swap rows with the row from the rows below which has the largest element in that column. This is an elementary operation.
•  Next, we perform another elementary operation. We subtract from each subsequent row     c=a(k,i)/a(i,i) times the row whose diagonal element large has been set now. So, every a(k,:)=a(k,:)-c*a(i,:). These two steps are repeated for all rows till we get a triangulated matrix.
• At this stage we can calculate the determinant by multiplying all the diagonal elements a(i,i) together. If the determinant is non-zero, we proceed with the next step. Else we return to the main program.
•  Finally, we perform back substitution. So, xn=bn/a(n,n). Knowing this, we calculate x(n-1)=[b(n-1)-a(n,n)*x(n)]/a(n-1,n-1)...and so on to find values of xn, x(n-1) through x1.We store this result in the array x(n).

Let us understand forward eliminate and back substitution further for a 3x3 coefficient matrix.

Forward Elimination and Triangularization

Row Swapping Explained further

Back Substitution to Calculate the Final Values of xn through x1 Backwards

Now for the FORTRAN program,

!To use arrays in the program
    module arrays
        real, allocatable :: a(:,:) , x(:), b(:)
    end module arrays

!To solve a system of equations of n variables using Gauss elimination method
program gausepivot
    use  arrays
    implicit none
    integer i,j,n,m,sing,flag
    real, allocatable :: LHS(:),RHS(:), check(:,:) !These arrays are optional. Only to verify  results
   
    write(*,*) "Give the size of the coefficient matrix i.e. the number of variables in the system of linear equations"
    read (*,*) n

   
    write(*,*) "The size of the coefficient matrix is ",n
   
    open(unit=1, file='gause.in')
    m=n+1
    allocate (a(n,m),x(n),b(n))
    write(*,*) "Producing the augmented matrix:"
   
    do i=1,n
      read(unit=1,fmt=*) (a(i,j),j=1,m)
      write(*,*) (a(i,j),j=1,m)
    end do
   
    !To verify results later. Optional
    allocate(check(n,m)) !To save the values of the array a(n,m) before they are altered
        do i=1,n
            do j=1,m
                  check(i,j)=a(i,j)
              end do
        end do

    call gauss (m,sing) !To call the subroutine that calculates the results
        if (sing==1) then
               write(*,*) "the coefficient matrix is singular"
               else ! To display results only if the matrix is non-singular
               write(*,*) "The given equation has the set of solutions in the order given as follows"
               write(*,*) (x(i),i=1,n)
        end if

    allocate (LHS(n),RHS(n)) !To test the values in array x(n) by substituting them in the linear equation. Optional
    do i=1,n
        LHS(i)=0
        do j=1,m-1
            LHS(i)=x(j)*check(i,j)+LHS(i)
        end do
       
        write(*,20) "LHS of row #",i,"=",LHS(i)
          20 format (a,i3,a,f20.8)
          RHS(i)=check(i,m)
          write(*,10) "RHS of row #",i,"=",RHS(i)
          10 format (a,i3,a,f20.8)
    end do
      
    flag=0 !To set the flag=1 if the L.H.S. is not equal to the R.H.S.
    do i=1,n
        if (abs(LHS(i)-RHS(i)).gt.00001) then
              flag=1
              exit
          end if
    end do
   
    if (flag==1) then !To write if results match or not
        write(*,*) "Results incorrect"
          else
          write(*,*) "The program works"
    end if
     
end program gausepivot

    subroutine gauss (m,sing) !Subroutine for the actual calculations
        use arrays
           implicit none
           integer :: n,m,sing,j,imax,i,k
           real:: temp(m),eps,flag,large,aa,det,sumax,c
       
           eps=1.0e-5
           !Forward Elimination
           !Triangulation
           n=m-1
        do i=1,n
               flag=0
             large=abs(a(i,i)) !To define the diagonal element of every row as 'large'
             do j=i+1,n
                   aa=abs(a(j,i)) !To check the subsequent elements in the same column below the row
                   if (aa>large) then !To check if the diagonal element is the largest in the column
                large=aa
                     flag=1
                     imax=j
                   end if
            end do
             !Row swapping
            if (flag==1) then !To swap rows so that the diagonal element of each row is larger than all elements in the same column in rows below it
                   temp=a(i,:)
                   a(i,:)=a(imax,:)
                   a(imax,:)=temp
               end if
            do k=i+1,n !To eliminate a variable from all subsequent rows
                c=a(k,i)/a(i,i) !Note that this value is fixed for a given k. (Elementary operation)
                a(k,:)=a(k,:)-c*a(i,:) !This will become zero when a(k,i)= a(k,i) from the row a(k,:)
            end do
        end do
                
        !Checking for singularity of coefficient matrix
        !Find det(a)
        det=1.0 !At this point the determinant is calculated simply by multiplying all the diagonal elements a(i,i)
        do i=1,n
            det=det*a(i,i)
           end do
                  
        if (abs(det)<eps) then !Calculating the determinant helps us find out if the matrix is singular
            sing=1 !In this case we will not get results
            return
            else
            sing=0
        end if
                      
        !Back substitution
        b=a(:,m)
        do i=n,1,-1 !This time we go from the bottom row to the top
            sumax=0.0
            do j=i+1,n !We look at all elements in the row to the right of the diagonal
                sumax=sumax+a(i,j)*x(j)
               end do
            x(i)=(b(i)-sumax)/a(i,i)
        end do
                          
    end subroutine gauss
                        

Sample Output for a System of 4 Linear Equations with 4  Variables