! Program to calculate the determinant of an n x n matrix
!  by using L-U decomposition
!
!
! Dec 8, 2011

!-----------------------------------------------------
! Real Kinds Module:
!------------------------------------------------------

module real_kinds

integer, parameter :: spt = kind(1.0  )
integer, parameter :: dpt = kind(1.0d0)

end module real_kinds

!-----------------------------------------------------
! Variables Module:
!------------------------------------------------------

module variables

use real_kinds

integer        :: i, j, k, dimensions, count1

character(len=32) :: formatting

real(kind=dpt) :: factor, det, temp

real(kind=dpt),allocatable :: matrix(:,:)



end module variables

!-----------------------------------------------------
! Main Program:
!------------------------------------------------------

program determinant

use real_kinds
use variables

implicit none

! get dimension of matrix:
write(*,*)
write(*,*) 'Enter the matrix dimension:'
read(*,*) dimensions

! allocate matrix:
allocate(matrix(dimensions,dimensions))

! get each row of matrix:
write(*,*)
write(*,*) 'Enter the matrix elements row by row:'
write(*,*)

do i=1,dimensions,1
   write(*,'("=> ")',advance='no')
   read(*,*) (matrix(i,j), j=1,dimensions,1)
enddo

! format for output of matrices
100 format ('(',i4,'(f10.5,1x))')
write(formatting,100) dimensions

write(*,*)
write(*,*)'Entered matrix:'
! display entered matrix:
do i=1,dimensions,1

  do j=1,dimensions,1
    write(*,formatting,advance='no') matrix(i,j)

  enddo
  write(*,*)

enddo

write(*,*)

! transform matrix into upper triangular:
count1 = 0
do i=1,dimensions-1,1
  if (matrix(i,i) == 0.0d0) then
    do k=i,dimensions-1,1
      if (matrix(k,i) /= 0.0d0) then
        do j=1,dimensions,1
          temp = matrix(i,j)
          matrix(i,j) = matrix(k,j)
          matrix(k,j) = temp
        enddo
      exit
      endif
    enddo
    count1 = count1 + 1  
  else
    do k=i+1,dimensions,1
      factor = -1.0d0 * matrix(k,i) / matrix(i,i)
      do j=i,dimensions,1
        matrix(k,j) = matrix(k,j) + factor * matrix(i,j)
      enddo
    enddo
  endif
enddo

!display upper triangular:
write(*,*)'Upper triangular:'
do i=1,dimensions,1

  do j=1,dimensions,1
    write(*,formatting,advance='no') matrix(i,j)

  enddo
  write(*,*)

enddo

! calculate determinant:

det = 1.0d0
do i=1,dimensions,1
  det = det * matrix(i,i)
enddo

! modify determinant:

if (mod(count1,2) == 0) then
  write(*,*)
  write(*,*)'Determinant: ',det
else
  write(*,*)
  write(*,*)'Determinant: ', -1.0d0*det
endif

write(*,*)
write(*,*)
write(*,*)'---Complete---'
write(*,*)

stop
end program determinant