!******************************************************* !* 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