invert_gauss_jordan Subroutine

  subroutine invert_gauss_jordan(a,n)
    implicit none
    integer,intent(in):: n
    !Should we make this assumed shape?
    !<DD>Tagged
    complex,dimension(n,n),intent(inout) :: a

    !Is it better to make these allocatable?
    !<DD>Tagged
    complex :: tmp,fac
    integer i, k

    do i=1,n
       fac=1/a(i,i) !This would become inverse if done on blocks
       a(i,i)=1
       a(:,i)=a(:,i)*fac
       !NOTE:These loops don't shrink as i increases so load should
       !be balanced?
       do k=1,i-1
          tmp=a(i,k)
          a(i,k)=0
          a(:,k)=a(:,k)-a(:,i)*tmp
       enddo
       do k=i+1,n
          tmp=a(i,k)
          a(i,k)=0
          a(:,k)=a(:,k)-a(:,i)*tmp
       enddo
    enddo
  end subroutine invert_gauss_jordan