C	ARTURO QUIRANTES SIERRA
C	Department of Applied Physics, Faculty of Sciences
C	University of Granada, 18071 Granada (SPAIN)
C	http://www.ugr.es/local/aquiran/codigos.htm
C	aquiran@ugr.es
C
C	Last update: 20 May 2.003
C
C	Subroutine MOE
C	to perform an LU-based inversion and multiplication 
C	of a square nq*nq matrix
C
C	You need to add a "CALL MOE(A,B,NQ)" in your main routine
C	Input:
C		A: First square matrix
C		B: Second square matrix
C		NQ: size of A, B matrices
C
C	Output
C
C		B: result of the product B*(A^-1) (B times inverse of A)
C
C	To be honest, what MOE does is a product X*(Y^-1)
C	The trick is to use transposition and calculate (At^-1)*Bt
C	which equals (B*(A^-1))t
C	we then transpose the result, and get B*(A^-1)
C
C	If, on the other hand, you do want to calculate (A^-1)*B
c	Just delete at the proper places (see remarks)	
C	This is a taylored adaptation of the ZGESV (=ZGETRF+ZGETRS)
C	subroutine from the LAPACK package 
c	For further info, see http://www.cs.colorado.edu/~lapack/
C	The maximum matrix size is 100*100
C	Of course, you can make it larger
C
	SUBROUTINE MOE(A,B,NQ)
	DOUBLE COMPLEX a(100,100),b(100,100),uno,cero,ztemp
	DOUBLE PRECISION SMAX		  
	INTEGER IPIV(100),NQ,J,K,LL,IMAX,JP,INFO,IP
	uno=(1.0d+0,0.0d+0)
	cero=(0.0d+0,0.0d+0)
	info=0
C	First of all, let's transpose the A, B, matrices
C	If you want to calculate (A^-1)*B, delete next 10 lines
	do j=1,nq
	do i=1,j-1
		ztemp=a(i,j)
		a(i,j)=a(j,i)
		a(j,i)=ztemp
		ztemp=b(i,j)
		b(i,j)=b(j,i)
		b(j,i)=ztemp
	end do
	end do
C	ZGETRF routine to make an LU decomposition
	DO 100 J=1,NQ
C	Search for pivoting and checking for singularity
C	Sub-subroutine IZAMAX to get the value of i such that a(i,j)=max for a given j
		nt=nq-j+1
		imax=1
		if(nt.eq.1) go to 20
                smax=cdabs(a(j,j))
		do 10 i=2,nt
			if(abs(a(j+i-1,j)).le.smax) go to 10
			imax=i
			smax=cdabs(a(j+i-1,j))
10		continue
20		continue
		nt=nq-j
C	End of sub-subroutine IZAMAX
		jp=j-1+imax
		ipiv(j)=jp
		if(a(jp,j).ne.cero) then
C	Sub-subroutine ZSWAP for row interchange
			if(jp.ne.j) then
				do i=1,nq
					ztemp=a(j,i)
					a(j,i)=a(jp,i)
					a(jp,i)=ztemp
				end do
C	End of sub-subroutine ZSWAP
			end if
C	Sub-subroutine ZSCAL to rescale: divide by por a(j,j)
			if(j.lt.nq) then
				ztemp=uno/a(j,j)
				do i=1,nt
					a(j+i,j)=ztemp*a(j+i,j)
				end do
C	End of sub-subroutine ZSCAL
			end if
C	Now if info=1 that means that U(l,l)=0 and invertion of A is not feasible
		else if(info.eq.0) then
			info=j
			write (*,*) 'MATRIZ U SINGULAR EN info=',J
			stop
		end if			 
		if(j.lt.nq) then
C	Sub-subroutine ZGERU para substitute A by A*x*y' (x,y vectors)
C	And therefore update the 'trailing submatrix'
		do k=1,nt
			if(a(j,j+k).ne.cero) then
				ztemp=-uno*a(j,j+k)
				do ll=1,nt
					a(j+ll,j+k)=a(j+ll,j+k)+a(j+ll,j)*ztemp
				end do
			end if
		end do 
C	End of sub-subroutine ZGERU
		end if
C	End of the J loop J
100	continue
C	End of subroutine ZGETRF
C
C	A is now an LU decomposition of a rowwise permutation of A
C	Diagonal unity elements belong tl T.  That is, Lii=1
C
C	You can check it by multiplying them, if you wish.
C
C	Resolving to obtain X = A^-1 * B
c	The procedure is to solve L*X=B, and then U*X=B
C	Sub-subroutine ZLASWP for row interchange in B	
	do i=1,nq
		ip=ipiv(i)
		if(ip.ne.i) then
			do k=1,nq
				ztemp=b(i,k)
				b(i,k)=b(ip,k)
				b(ip,k)=ztemp
			end do
		end if
c	End of subroutine ZLASWP
	end do	
c	Subroutine ZTRSM to solve L*X=B.  The result is overwritten in B
	do j=1,nq
		do k=1,nq
			if(b(k,j).ne.cero) then
				do i=k+1,nq
					b(i,j)=b(i,j)-b(k,j)*a(i,k)
				end do
			end if
		end do
	end do
c	Subroutine ZTRSM to solve U*X=B.  The final result is overwritten in B
	do j=1,nq
		do k=nq,1,-1
			if(b(k,j).ne.cero) then
				b(k,j)=b(k,j)/a(k,k)
				do i=1,k-1
					b(i,j)=b(i,j)-b(k,j)*a(i,k)
				end do
			end if
		end do
	end do
C	End of subroutine ZTRSM
C	And now, let's transpose again
C	The final result, B, equals the product B*A^(-1)
C	If you want to calculate (A^-1)*B, delete next 7 lines
	do j=1,nq
	do i=1,j-1
		ztemp=b(i,j)
		b(i,j)=b(j,i)
		b(j,i)=ztemp
	end do
	end do	 
C	End of subroutine MOE
	return
	end