Code: GAUSS

 

!  The below subroutine GAUSS_ELIM solves the input linear equations
!  via gaussian elimination withpartial pivoting.

       program gauss

       interface
         subroutine gauss_elim(c,n)
           integer n
           real, dimension(1:n, 1:n+1) :: c
         end subroutine
       end interface

       parameter(n =256)
       integer i
       real, dimension(n*(n+1))::vector
       real, dimension(n, n+1)::augmented

!HPF$ DISTRIBUTE augmented (BLOCK,*)

       open (10, file="datafile", status = "OLD")
       do i = 1,n
          read(10,*) augmented(i,1:n+1)
       enddo


       call gauss_elim(augmented,n)
       
       print *, '   '

       end

!**********************************************************************
       subroutine gauss_elim(mat,n)
       integer n
       real, dimension (1:n, 1:n+1)  :: mat
       ! local var
       integer i,j,pivot_index
       integer, dimension (1) :: pivot_index_array
       real, dimension (1:n+1) :: temprow

!HPF$  INHERIT MAT

       !  main loop for gauss elimination

       do i = 1,n
          ! find pivot index
          pivot_index_array = maxloc(abs(mat(i:n,i)))
          pivot_index = pivot_index_array(1) + i - 1
          ! if pivot row not already in place, make it so
          if( pivot_index .ne. i) then
            temprow = mat(i,1:n+1)
             mat(i,1:n+1) = mat(pivot_index,1:n+1)
             mat(pivot_index,1:n+1) = temprow
          endif
          ! check for singular matrix
          if(abs(mat(i,i)) .lt. 1.0e-6) then
             write(*,*) 'error: matrix is singular to machine precision'
             stop
          endif
          mat(i,i:n+1) = mat(i,i:n+1)/mat(i,i)
          do j = i+1,n
             mat(j,i:n+1) = mat(j,i:n+1) - mat(j,i) * mat(i,i:n+1)
          enddo
       enddo

       !  back substitution

       do i = n-1,1,-1
          mat(i,n+1) = mat(i,n+1) - sum(mat(i,i+1:n) * mat(i+1:n,n+1))
          mat(i,i+1:n) = 0.0
       enddo

       end