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: 16 February 2005	
C
C	Subroutine WILLIE
C	to perform an LU-based inversion and multiplication 
C	of a square nq*nq matrix
C
C	You need to add a "CALL WILLIE(AR,AI,BR,BI,NQ)" in your main routine
C	Input:
C		Ar + i*AI: First two square matrices  (A = Ar + i*AI )
C		Br + i*Bi: Second two square matrices (B = Br + i*Bi )
C		NQ: size of A,B matrices
C
C	Output
C
C		Br + i*BI: result of the product B*(A^-1) (B times inverse of A)
C
C	This is a modification of subroutine MOE, the difference being that the
C	input and output matrices are declared as real and imaginary part
C	This is useful when you are calculating a complex matrix by defining
C	real and imaginary parts
C	Remember than operations with complex number are slower,
C	Of course, you can also use it to multiply+invert real matrices
C	(AI,BI = 0), but you might prefer to use MOE and declare MOE's
C	matrices as real instead.  Your choice.
C
C	To be honest, what WILLIE 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 WILLIE(AR,AI,BR,BI,NQ)
	DOUBLE PRECISION AR(100,100),AI(100,100),BR(100,100),BI(100,100)
	DOUBLE PRECISION SMAX,ZTEMP,ZTEMPX,ZTEMPY,ZTEMP1,ZTEMP2,uno,cero
	INTEGER IPIV(ST2),NQ,J,K,LL,IMAX,JP,INFO,IP
	UNO=1.0D0
	CERO=0.0D0
	info=0
C	First of all, 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=ar(i,j)
		ar(i,j)=ar(j,i)
		ar(j,i)=ztemp
		ztemp=ai(i,j)
		ai(i,j)=ai(j,i)
		ai(j,i)=ztemp
		ztemp=br(i,j)
		br(i,j)=br(j,i)
		br(j,i)=ztemp
		ztemp=bi(i,j)
		bi(i,j)=bi(j,i)
		bi(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=ar(j,j)*ar(j,j)+ai(j,j)*ai(j,j)
		do 10 i=2,nt
			ztemp=ar(j+i-1,j)*ar(j+i-1,j)+ai(j+i-1,j)*ai(j+i-1,j)
			if(ztemp.le.smax) go to 10
			imax=i
			smax=ztemp
10		continue
20		continue
		nt=nq-j
C	End of sub-subroutine IZAMAX
		jp=j-1+imax
		ipiv(j)=jp
		if(ar(jp,j).ne.cero.and.ai(jp,j).ne.cero) then
C	Sub-subroutine ZSWAP for row interchange
			if(jp.ne.j) then
				do i=1,nq
					ztemp=ar(j,i)
					ar(j,i)=ar(jp,i)
					ar(jp,i)=ztemp
					ztemp=ai(j,i)
					ai(j,i)=ai(jp,i)
					ai(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=ar(j,j)*ar(j,j)+ai(j,j)*ai(j,j)
				ztempx=ar(j,j)/ztemp
				ztempy=-ai(j,j)/ztemp
				do i=1,nt
				ztemp1=ar(j+i,j)*ztempx-ai(j+i,j)*ztempy
				ztemp2=ar(j+i,j)*ztempy+ai(j+i,j)*ztempx
				ar(j+i,j)=ztemp1
				ai(j+i,j)=ztemp2
				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 m,info=',m,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(ar(j,j+k).ne.cero.and.ai(j,j+k).ne.cero) then
			ztempx=-ar(j,j+k)
			ztempy=-ai(j,j+k)
			do ll=1,nt
				ztemp1=ar(j+ll,j)*ztempx-ai(j+ll,j)*ztempy
				ztemp2=ar(j+ll,j)*ztempy+ai(j+ll,j)*ztempx
				ar(j+ll,j+k)=ar(j+ll,j+k)+ztemp1
				ai(j+ll,j+k)=ai(j+ll,j+k)+ztemp2
			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	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=br(i,k)
				br(i,k)=br(ip,k)
				br(ip,k)=ztemp
				ztemp=bi(i,k)
				bi(i,k)=bi(ip,k)
				bi(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(br(k,j).ne.cero.and.bi(k,j).ne.cero) then
			do i=k+1,nq
				ztempx=br(k,j)*ar(i,k)-bi(k,j)*ai(i,k)
				ztempy=br(k,j)*ai(i,k)+bi(k,j)*ar(i,k)
				br(i,j)=br(i,j)-ztempx
				bi(i,j)=bi(i,j)-ztempy
			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(br(k,j).ne.cero.and.bi(k,j).ne.cero) then
			ztemp=ar(k,k)*ar(k,k)+ai(k,k)*ai(k,k)
			ztempx=(ar(k,k)*br(k,j)+ai(k,k)*bi(k,j))/ztemp
			ztempy=(ar(k,k)*bi(k,j)-ai(k,k)*br(k,j))/ztemp
			br(k,j)=ztempx
			bi(k,j)=ztempy
			do i=1,k-1
				ztemp1=br(k,j)*ar(i,k)-bi(k,j)*ai(i,k)
				ztemp2=br(k,j)*ai(i,k)+bi(k,j)*ar(i,k)
				br(i,j)=br(i,j)-ztemp1
				bi(i,j)=bi(i,j)-ztemp2
			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)
	do j=1,nq
	do i=1,j-1
		ztemp=br(i,j)
		br(i,j)=br(j,i)
		br(j,i)=ztemp
		ztemp=bi(i,j)
		bi(i,j)=bi(j,i)
		bi(j,i)=ztemp
	end do
	end do	 
C	End of subroutine LULABY
	return
	end