! $Id: graph_mod.f90 194 2011-02-10 11:19:11Z fhrouet $

module graph

  use utils
  implicit none

  type cell
    integer             :: index  ! row index in column
    integer             :: lfill  ! level of fill 
    real(kind(1.d0))    :: value  ! real value 
    type (cell),pointer :: next   ! next entry in column
  end type cell

contains

  subroutine csc_to_graph(A, g, trace)
    ! Creates the graph of matrix (in CSC format,
    ! lower part only). NB: j is in adjacency list
    ! of i only if j>i (edge j->i is not needed for
    ! our purpose).
    ! g(i)%index is the number of elements in column
    ! i (lower part). g(i)%next is the first element
    ! in column i (lower part); since A is pos. def.
    ! its diagonal elements are non zeros, thus
    ! g(i)%next is A(i,i).
    ! Property: each adjacency list follows the same
    ! order as in the corresponding column in the
    ! CSC structure; thus, if each column is sorted
    ! in the CSC structure, each list is sorted in the
    ! graph.

    type (csc_mat),intent(in) :: A
    type (cell),dimension(:), allocatable, intent(out) :: g
    integer, intent(in) :: trace

    type (cell),pointer::current, prev
    integer :: allocok,i,j

    write(*,*)'  Creating graph...'
    allocate(g(A%n), stat=allocok)
    if (allocok.ne.0) then
      write(*,*)'Problem with allocation of g in create_graph', allocok
      stop
    end if
    do i=1,A%n
      nullify(g(i)%next)
    enddo

    do j=1,A%n
      g(j)%index = 0
      g(j)%lfill = 0
      nullify(prev)
      do i=A%jptr(j+1)-1,A%jptr(j),-1
        allocate(current,stat=allocok)
        if (allocok.ne.0) then
          write(*,*)'Problem with allocation of a cell in create_graph',allocok
          stop
        end if
        current%index=A%irn(i)
        current%value=A%val(i)
        current%next=>prev
        prev=>current
        g(j)%index=g(j)%index+1
      end do
      g(j)%next=>prev
    end do

    if(trace.eq.2) then
      write(*,*)'  Graph:'
      do i=1,A%n
        write(*,'(A,I5,A,I5,A)',advance='no')'     Node ',i,'; index: ',g(i)%index,'; adjacency: '
        current => g(i)%next
        do j=1,g(i)%index
          write(*,'(I5,A)',advance='no')current%index,' '
          current => current%next
        end do
        write(*,*)
      end do
    end if

    return

  end subroutine csc_to_graph

  !------------------------------------------------------------------------------

  subroutine destroy_graph(g,dim)

    integer, intent(in)  :: dim
    type (cell),allocatable, dimension(:), intent(inout) ::g

    type (cell),pointer::next
    integer  :: j

    write(*,*)'  Destroying graph...'
    do j=1,dim
      do while (associated(g(j)%next))
        next => g(j)%next
        g(j)%next => next%next
        deallocate(next)
      end do
    end do
    deallocate(g)

  end subroutine destroy_graph

  !------------------------------------------------------------------------------

  subroutine graph_to_csc(g, dim, L, trace)
  ! Converts a graph to the corresponding CSC structure
  ! (lower part). Values are initialized to 0.

    integer, intent(in) :: dim, trace
    type(cell), dimension(dim), intent(in) :: g
    type(csc_mat), intent(out) :: L

    integer :: i, j, nz, allocok
    type(cell),pointer :: current

    write(*,*)'  Creating L...'

    nz = 0
    do i=1,dim
      nz = nz + g(i)%index
    end do

    L%n = dim
    L%nz = nz

    allocate(L%jptr(L%n+1),L%irn(L%nz),L%val(L%nz), stat=allocok)
    if (allocok.ne.0) then
      write(*,*)'Problem with allocation of tmp arrays in ordering', allocok
      stop
    end if

    L%jptr(1) = 1
    nz = 0
    do i=1, dim
      current => g(i)%next
      do j=1,g(i)%index
        nz = nz+1
        L%irn(nz) = current%index
        L%val(nz) = current%value
        current => current%next
      end do
      L%jptr(i+1) = nz + 1
    end do

    if(trace.eq.2) then
9999   format (a,10i6/(8x, 10i6))
9998   format (a,5f12.5/(8x, 5f12.5))
      write(*,9999) '     L%jptr  =', L%jptr
      write(*,9999) '     L%irn   =', L%irn
      write(*,9998) '     L%val   =', L%val
    end if

  end subroutine graph_to_csc

end module graph
