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