Reverse Cuthill-MacKee subroutine

cuthill_mckee.f

subroutine cuthill_mckee(matrix,LDA,n,iorder)
  implicit none

!  Uses the Cuthill-McKee algorithm to reorder the columns and rows of the matrix

  integer, intent(in)            :: n,LDA
  integer, intent(out)           :: iorder(n)
  double precision, intent(in)   :: matrix(LDA,n)

  integer                        :: i,j,k
  integer, allocatable           :: queue(:), degree(:), connexions(:,:)
  double precision, parameter    :: threshold = 1.d-20
  logical, allocatable           :: available(:)
  integer                        :: ip, ir, ic, iq, iq0

  allocate(queue(n*n),degree(n),connexions(n,n),available(n))

  iorder = 0
  do j=1,n
    degree(j) = 0
    do i=1,n
      if ( (i /= j).and.(dabs(matrix(i,j)) > threshold) ) then
        degree(j) = degree(j)+1
        connexions(degree(j),j) = i
      endif
    enddo
  enddo

  available(:) = .True.

  ir = n
  iq = 1
  iq0 = 1

  do while (ir > 0)
    ! Select node with lowest degree
    ip = minloc(degree,1)
    if (.not.available(ip)) then
      degree(ip) = huge(1)
      cycle
    endif

    ! Add ip to the result
    iorder(ir) = ip
    ir -= 1

    ! Mark ip as done
    available(ip) = .False.

    ! Add to the queue all the nodes connected to ip
    do j=1,degree(ip)
      queue(iq) = connexions(j,ip)
      iq = iq+1
    enddo

    do while (iq0 < iq)
      ic = queue(iq0)
      iq0 = iq0+1
      if (available(ic)) then
        available(ic) = .False.
        iorder(ir) = ic
        ir = ir-1
        do j=1,degree(ic)
          queue(iq) = connexions(j,ic)
          iq = iq+1
        enddo
      endif
    enddo
  enddo

  deallocate(queue,degree,connexions,available)
end