C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
C     File:        zgaussj.f (Fortran 77 source for the zgaussj routine)
C     Author:      Fredrik Jonsson <fj@phys.soton.ac.uk>
C     Date:        December 26, 2005
C     Last change: December 29, 2006
C     Description: The zgaussj(a,n,np,b,m,mp) routine solves systems of
C                  linear equations by Gauss-Jordan elimination, in similar
C                  to as in Eq. (2.1.1) of W.H. Press et al, "Numerical
C                  Recipes in Fortran 77, The Art of Scientific Computing",
C                  2nd Edn. (Cambridge University Press, Cambridge, 1992),
C                  but here modified so as to solve complex-valued systems.
C                     Here a(1:n,1:n) is an input matrix stored in an array
C                  of physical dimensions np by np. b(1:n,1:m) is an input
C                  matrix containing the m right-hand side vectors, stored
C                  in an array of physical dimensions np by mp. On output,
C                  a(1:n,1:n) is replaced by its matrix inverse, and
C                  b(1:n,1:m) is replaced by the corresponding set of solu-
C                  tion vectors. The parameter nmax contains the largest
C                  anticipated value of n.
C
C     Copyright (C) 2006, Fredrik Jonsson <fj@phys.soton.ac.uk>
C...+....1....+....2....+....3....+....4....+....5....+....6....+....7....+
      subroutine zgaussj(a,n,np,b,m,mp)
      integer m,mp,n,np,nmax
      complex*8 a(np,np),b(np,mp)
      parameter (nmax=50)
      integer i,icol,irow,j,k,l,ll,indxc(nmax),indxr(nmax),ipiv(nmax)
      complex*8 big,dum,pivinv
      do 11 j=1,n
         ipiv(j)=0
 11   continue
      do 22 i=1,n
         big=0.
         do 13 j=1,n
            if(ipiv(j).ne.1)then
               do 12 k=1,n
                  if (ipiv(k).eq.0) then
                     if (abs(a(j,k)).ge.abs(big))then
                        big=abs(a(j,k))
                        irow=j
                        icol=k
                     endif
                  else if (ipiv(k).gt.1) then
                     pause 'Singular matrix ecountered in zgaussj()!'
                  endif
 12            continue
            endif
 13      continue
         ipiv(icol)=ipiv(icol)+1
         if (irow.ne.icol) then
            do 14 l=1,n
               dum=a(irow,l)
               a(irow,l)=a(icol,l)
               a(icol,l)=dum
 14         continue
            do 15 l=1,m
               dum=b(irow,l)
               b(irow,l)=b(icol,l)
               b(icol,l)=dum
 15         continue
         endif
         indxr(i)=irow
         indxc(i)=icol
         if (abs(a(icol,icol)).eq.(0.0))
            pause 'Singular matrix ecountered in zgaussj()! [Case 2]'
         pivinv=1.0/a(icol,icol)
         a(icol,icol)=1.
         do 16 l=1,n
            a(icol,l)=a(icol,l)*pivinv
 16      continue
         do 17 l=1,m
            b(icol,l)=b(icol,l)*pivinv
 17      continue
         do 21 ll=1,n
            if(ll.ne.icol)then
               dum=a(ll,icol)
               a(ll,icol)=0.
               do 18 l=1,n
                  a(ll,l)=a(ll,l)-a(icol,l)*dum
 18            continue
               do 19 l=1,m
                  b(ll,l)=b(ll,l)-b(icol,l)*dum
 19            continue
            endif
 21      continue
 22   continue
      do 24 l=n,1,-1
         if(indxr(l).ne.indxc(l))then
            do 23 k=1,n
               dum=a(k,indxr(l))
               a(k,indxr(l))=a(k,indxc(l))
               a(k,indxc(l))=dum
 23         continue
         endif
 24   continue
      return
      end