Coding the LU decomposition

The F90 routine is given , the coefficient matrix; it returns the LU factorization, stored compactly in and the integer vector of permutations, or pivots, . Note that the content of is overwritten by this routine.

An LU decomposition routine might look like:

SUBROUTINE LU_Fact(A,Pivot,Ierr)
IMPLICIT NONE
!
! The arguments:
!
! The input matrix, at exit
! will hold the LU factorization
REAL, DIMENSION(:,:), INTENT(INOUT) :: A
! Vector of permutations
INTEGER, DIMENSION(:), INTENT(OUT)   :: Pivot
! Singularity indicator, = 0 if A nonsingular,
! and = column number j if the first zero
! pivot was detected in column j
INTEGER, INTENT(OUT)   :: Ierr
!
! The local variables:
!
LOGICAL :: singular        ! Singularity flag 
INTEGER :: i,j,k,n         ! DO loop variables
INTEGER, DIMENSION(:) :: p ! pivot location
REAL, DIMENSION(:) :: Tmp  ! Temporary row
REAL  :: uround            ! rounding unit
!
! Check if the argument is a square matrix
IF( size(A,1).NE.size(A,2) ) THEN
  PRINT*,"Error in Factorize: A must be square"
  RETURN
END IF
n=SIZE(A,1) ! the dimension of the matrix

ALLOCATE(Tmp(n),stat = k)
IF (k.NE.0) THEN
  PRINT*,"Error in Factorize: cannot allocate Tmp"; RETURN
END IF
!
ALLOCATE(p(n),stat = k)
IF (k.NE.0) THEN
  PRINT*,"Error in Factorize: cannot allocate p"; RETURN
END IF
!
Ierr = 0                      ! reset error indicator
singular = .FALSE.            ! reset singularity flag
uround = 1.0E-7               ! unit roundoff, set it to
                              ! 1.0D-14 for double precision
!
DO j=1,n-1  ! columns 1:n-1
!
  p=MAXLOC(ABS(A(j:n,j)))+j-1 ! Look for pivot, A(p(1),j)
  IF (p(1).NE.j) THEN         ! If pivot is non-diagonal
     Tmp(:)    = A(j,:)       ! permute rows j and p(1)
     A(j,:)    = A(p(1),:)
     A(p(1),:) = Tmp(:)
     Pivot(j) = p(1)          ! Save the pivot position
  ELSE
     Pivot(j) = j             ! Save the pivot position
  END IF
!
! If A(j,j)=0 then the matrix is singular
! uround is the rounding unit (machine epsilon)
  IF ( ABS(A(j,j)) .LT. uround  ) THEN
    Ierr = j                  ! Singularity Position
    singular = .TRUE.         ! Singularity Flag
    exit                      ! the `DO j' loop
  END IF
!
  DO i=j+1,n                  ! rows to be processed
     A(i,j) = A(i,j)/A(j,j)   ! Store the multiplier
     A(i,j+1:n) = A(i,j+1:n) - A(i,j)*A(j,j+1:n)
  END DO          ! Row loop
!
END DO            ! Column loop
!
! If A(n,n)=0 then the matrix is singular
! uround is the rounding unit (machine epsilon)
IF ( abs(A(n,n)) .LT. uround  ) THEN
    Ierr = n                  ! Singularity Flag
    singular = .TRUE.
END IF
!  
IF (allocated(Tmp)) DEALLOCATE(Tmp)
!
IF (.NOT. singular) THEN
   Ierr = 0
END IF