!*******************************************************
!* LU decomposition routines used by test_clu.f90 *
!* *
!* F90 version by J-P Moreau, Paris *
!* (with extension to complex domain) *
!* (www.jpmoreau.fr) *
!* --------------------------------------------------- *
!* Reference: *
!* *
!* "Numerical Recipes by W.H. Press, B. P. Flannery, *
!* S.A. Teukolsky and W.T. Vetterling, Cambridge *
!* University Press, 1986". *
!* *
!*******************************************************
MODULE CLU
CONTAINS
!*****************************************************************
!* Given an N x N complex matrix A, this routine replaces it by *
!* the LU decomposition of a rowwise permutation of itself. A *
!* and N are input. INDX is an output vector which records the *
!* row permutation effected by the partial pivoting; D is output *
!* as -1 or 1, depending on whether the number of row inter- *
!* changes was even or odd, respectively. This routine is used in*
!* combination with LUBKSB to solve linear equations or to invert*
!* a matrix. Return code is 1, if matrix is singular. *
!*****************************************************************
Subroutine CLUDCMP(A,N,INDX,D,CODE)
REAL*8, PARAMETER :: TINY=2.2D-16
COMPLEX AMAX,DUM, SUM, A(N,N)
INTEGER CODE, D, INDX(N)
COMPLEX, pointer :: VV(:) !complex vector (n)
allocate(VV(N),stat=ialloc)
D=1; CODE=0
DO I=1,N
AMAX=CMPLX(0.,0.)
DO J=1,N
IF (ABS(A(I,J)).GT.ABS(AMAX)) AMAX=A(I,J)
END DO ! j loop
IF(ABS(AMAX).LT.TINY) THEN
CODE = 1
RETURN
END IF
VV(I) = 1. / AMAX
END DO ! i loop
DO J=1,N
DO I=1,J-1
SUM = A(I,J)
DO K=1,I-1
SUM = SUM - A(I,K)*A(K,J)
END DO ! k loop
A(I,J) = SUM
END DO ! i loop
AMAX = CMPLX(0.,0.)
DO I=J,N
SUM = A(I,J)
DO K=1,J-1
SUM = SUM - A(I,K)*A(K,J)
END DO ! k loop
A(I,J) = SUM
DUM = VV(I)*ABS(SUM)
IF(CABS(DUM).GE.CABS(AMAX)) THEN
IMAX = I
AMAX = DUM
END IF
END DO ! i loop
IF(J.NE.IMAX) THEN
DO K=1,N
DUM = A(IMAX,K)
A(IMAX,K) = A(J,K)
A(J,K) = DUM
END DO ! k loop
D = -D
VV(IMAX) = VV(J)
END IF
INDX(J) = IMAX
IF(ABS(A(J,J)) < TINY) A(J,J) = CMPLX(TINY,0.)
IF(J.NE.N) THEN
DUM = 1. / A(J,J)
DO I=J+1,N
A(I,J) = A(I,J)*DUM
END DO ! i loop
END IF
END DO ! j loop
RETURN
END subroutine CLUDCMP
!********************************************************************
!* Solves the set of N complex linear equations A . X = B. Here A *
!* is input, not as the matrix A but rather as its LU decomposition,*
!* determined by the routine CLUDCMP. INDX is input as the permuta- *
!* tion vector returned by CLUDCMP. B is input as the right-hand *
!* side complex vector B, and returns with the solution vector X. A,*
!* N and INDX are not modified by this routine and can be used for *
!* successive calls with different right-hand sides. This routine is*
!* also efficient for plain complex matrix inversion. *
!********************************************************************
Subroutine CLUBKSB(A,N,INDX,B)
COMPLEX SUM, A(N,N),B(N)
INTEGER INDX(N)
II = 0
DO I=1,N
LL = INDX(I)
SUM = B(LL)
B(LL) = B(I)
IF(II.NE.0) THEN
DO J=II,I-1
SUM = SUM - A(I,J)*B(J)
END DO ! j loop
ELSE IF(CABS(SUM).NE.0.) THEN
II = I
END IF
B(I) = SUM
END DO ! i loop
DO I=N,1,-1
SUM = B(I)
IF(I < N) THEN
DO J=I+1,N
SUM = SUM - A(I,J)*B(J)
END DO ! j loop
END IF
B(I) = SUM / A(I,I)
END DO ! i loop
RETURN
END subroutine CLUBKSB
!************************************************************
!* Inversion of a complex square matrix by LU decomposition *
!* -------------------------------------------------------- *
!* INPUTS: *
!* A: complex square matrix of size n x n *
!* n: size of matrix A *
!* OUTPUTS: *
!* Y; Inverse of A of size n x n *
!* rc: error code (must be zero) *
!************************************************************
Subroutine CINVERT_LU(A,n,Y,rc)
COMPLEX A(n,n), Y(n,n)
integer rc, D
integer,pointer :: INDX(:) !integer vector (n)
!dynamic allocations
allocate(INDX(n),stat=ialloc)
!call LU decomposition routine
call CLUDCMP(A,n,INDX,D,rc)
!call solver if previous return code is ok
!to obtain inverse of A one column at a time
if (rc.eq.0) then
do j=1, n
call CLUBKSB(A,n,INDX,Y(1,j))
end do
end if
!the inverse matrix is now in matrix Y
!the original matrix A is destroyed
return
END Subroutine CInvert_lu
END MODULE CLU
! end of file clu.f90