SolveLinearSystems Subroutine

  subroutine SolveLinearSystems(n,a,b)
    implicit none
    integer,intent(in) :: n
    real(real64),intent(inout) :: a(1:n,1:n)
    real(real64),intent(inout) :: b(1:n,1:n)
    ! Local
    integer :: i,j,p,ipiv
    real(real64) :: pivmax,factor,swap(1:n)

    do j = 1,n-1

      ! Partial pivoting
      ipiv = j
      pivmax = abs(a(j,j))
      do p = j+1,n
        if(abs(a(p,j)) > pivmax) then
          ipiv = p
          pivmax = abs(a(p,j))
        endif
      enddo
      if(ipiv /= j) then
        swap = a(j,1:n)
        a(j,1:n) = a(ipiv,1:n)
        a(ipiv,1:n) = swap
        swap = b(j,1:n)
        b(j,1:n) = b(ipiv,1:n)
        b(ipiv,1:n) = swap
      endif

      ! Elimination
      do i = j+1,n
        factor = a(i,j)/a(j,j)
        a(i,j:n) = a(i,j:n)-factor*a(j,j:n)
        b(i,1:n) = b(i,1:n)-factor*b(j,1:n)
      enddo

    enddo

    ! Back substitution
    do j = n,1,-1
      do i = j+1,n
        b(j,1:n) = b(j,1:n)-a(j,i)*b(i,1:n)
      enddo
      b(j,1:n) = b(j,1:n)/a(j,j)
    enddo

  endsubroutine SolveLinearSystems