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