graph_mod.F90

Go to the documentation of this file.

module graph_mod

  use choral_variables
  use basic_tools
  use algebra_set
  !$ use OMP_LIB
  
  implicit none
  private


  !       %----------------------------------------%
  !       |                                        |
  !       |          PUBLIC DATA                   |
  !       |                                        |
  !       %----------------------------------------%
  !! TYPE
  public :: graph

  !! GENERIC
  public :: clear          ! TESTED
  public :: valid          ! TESTED
  public :: copy           ! TESTED
  public :: equal          ! TESTED
  public :: sortrows       ! TESTED
  public :: getrow         ! TESTED
  public :: setrow         ! TESTED
  public :: transpose      ! TESTED
  public :: print
  
  !! NON GENERIC
  public :: prod_tab       ! TESTED
  public :: graph_extract  ! TESTED
  public :: union          ! TESTED
  public :: graph_prod     ! TESTED


  !       %----------------------------------------%
  !       |                                        |
  !       |          DERIVED TYPE                  |
  !       |                                        |
  !       %----------------------------------------%
  type graph
     integer :: nl     =  0
     integer :: nc     =  0
     integer :: nnz    =  0 
     integer, dimension(:), allocatable :: row
     integer, dimension(:), allocatable :: col

     integer :: width  = 0

     logical :: sorted = .false.

   contains
     
     final :: graph_clear

  end type graph

  !       %----------------------------------------%
  !       |                                        |
  !       |       GENERIc SUBROUTINES              |
  !       |                                        |
  !       %----------------------------------------%
  interface clear
     module procedure graph_clear     ! TESTED
  end interface clear

  interface graph
     module procedure graph_create    ! TESTED
  end interface graph

  interface valid
     module procedure graph_valid   ! TESTED
  end interface valid

  interface copy
     module procedure graph_copy        ! TESTED
  end interface copy

  interface print
     module procedure graph_print
  end interface print

  interface equal 
     module procedure  graph_equal      ! TESTED
  end interface equal

  interface sortrows
     module procedure graph_sortrows    ! TESTED
  end interface sortrows

  interface setrow
     module procedure graph_setrow      ! TESTED
  end interface setrow

  interface getrow
     module procedure graph_getrow_1    ! TESTED
     module procedure graph_getrow_2    ! TESTED
  end interface getrow

  interface transpose
     module procedure graph_transpose   ! TESTED
  end interface transpose

contains

  subroutine graph_clear(g)
    type(graph), intent(inout) :: g

    call freemem(g%row)
    call freemem(g%col)
    g%nl     = 0
    g%nc     = 0
    g%nnz    = 0
    g%width  = 0
    g%sorted = .false. 

  end subroutine graph_clear


  function graph_create(nnz, nl, nc) result(g)
    type(graph)                           :: g
    integer, dimension(nl), intent(in)    :: nnz
    integer               , intent(in)    :: nl, nc

    integer :: nz, ii, jj, wdt, length

    if ( nc < 0      ) call quit( &
         & "graph_mod: graph_create: number of columns nc < 0")

    if ( nl <= 0      ) call quit( &
         & "graph_mod: graph_create: number of lines nl <= 0")

    if ( size(nnz,1) /= nl     ) call quit( &
         & "graph_mod: graph_create: 'nnz'wrong size")

    if ( .NOT. all( nnz >=0 ) ) call quit( &
         & "graph_mod: graph_create: 'nnz(ii)' must be >=0 ")

    if ( sum( nnz ) == 0      ) call quit( &
         & "graph_mod: graph_create: one line at least must be non-empty'")

    call clear(g)
    g%nl  = nl
    g%nc  = nc
    nz    = sum(nnz)
    g%nnz = nz
    call allocmem(g%row, g%nl+1)
    call allocmem(g%col, nz    )

    g%row = 0
    g%col = 0
    wdt   = 0

    ! Construction of g%row
    jj = 1
    do ii = 1, g%nl
       g%row(ii) = jj
       length = nnz(ii)
       if (length>wdt) wdt = length
       jj = jj + length
    end do
    g%row(g%nl+1) = jj

    g%width = wdt

  end function graph_create


  function graph_valid(g) result(info)
    type(graph), intent(in) :: g
    integer :: info

    integer :: ii, wdt, length

    info = 1

    if ( .NOT. (allocated(g%row)) ) info = -1
    if ( .NOT. (allocated(g%col)) ) info = -2

    if (info > 0) then
       if ( g%nl          <= 0 ) info = -3 
       if ( g%nc          <= 0 ) info = -4
       if ( g%nnz         <= 0 ) info = -5
       if ( minval(g%row) <= 0 ) info = -6
       if ( minval(g%col) <  0 ) info = -7
    end if

    if (info > 0) then
       if ( maxval(g%row) /= g%nnz+1 ) info = -8
       if ( maxval(g%col)  > g%nc    ) info = -9
       if ( size(g%row,1) /= g%nl+1  ) info = -10
       if ( size(g%col,1) /= g%nnz   ) info = -11
    end if

    if (info > 0) then
       if ( g%row(g%nl+1) /= size(g%col)+1 ) info = -12 
    end if

    if (info > 0) then
       wdt = 0
       do ii = 1, g%nl
          length = g%row(ii+1) - g%row(ii)
          if (length>wdt) wdt = length
       end do
       if ( wdt /= g%width ) info = -13 
    end if

    if (choral_verb>2) then
       write(*,*)"graph_mod       : valid,    info =", int(info, 1)
    end if

  end function graph_valid


  subroutine graph_print(g)
    type(graph), intent(in) :: g

    write(*,*) 'graph_mod  : print'
 
    write(*,*) '  Valid                 =', int( valid(g), 1)
    write(*,*) '  Nb lines              =', g%nl
    write(*,*) '  Nb columns            =', g%nc
    write(*,*) '  Nb non-zero entries   =', g%nnz
    write(*,*) '  Maximal line width    =', g%width
    write(*,*) '  Sorted lines          = ', g%sorted
    
  end subroutine graph_print
  
  subroutine graph_copy(g, g0)
    type(graph) , intent(inout) :: g
    type(graph), intent(in)    :: g0

    call clear(g)
    g%nl     = g0%nl
    g%nc     = g0%nc
    g%nnz    = g0%nnz
    g%width  = g0%width
    g%sorted = g0%sorted
    
    call allocmem(g%row, size(g0%row,1) )
    call allocmem(g%col, size(g0%col,1) )
    g%row = g0%row
    g%col = g0%col
    
  end subroutine graph_copy
    

  function graph_equal(g2, g1) result(b)
    logical                :: b
    type(graph), intent(in):: g1, g2

    integer, dimension(:), allocatable :: p1, p2
    integer :: ii, wdt, sz

    b = ( valid(g1)>0 ) .AND. ( valid(g2)>0 )
    if(.NOT.b) return

    b = b .AND. ( g1%nl == g2%nl ) .AND. ( g1%nnz == g2%nnz )
    b = b .AND. ( g1%nc == g2%nc )
    b = b .AND. ( g1%width == g2%width )
    if(.NOT.b) return

    b = b .AND. (all( g1%row == g2%row ) )
    if(.NOT.b) return

    wdt = g1%width
    call allocmem( p1, wdt)
    call allocmem( p2, wdt )
    
    do ii=1, g1%nl
       call getrow(sz, p1, wdt, g1, ii)
       if( sz == 0 ) cycle

       call getrow(sz, p2, wdt, g2, ii)

       call sort( p1(1:sz), sz )
       call sort( p2(1:sz), sz )
       
       b = b .AND. (all( p1(1:sz) == p2(1:sz) ) )
       
    end do

  end function graph_equal



  subroutine graph_getrow_1(sz, p, n, g, ii)
    integer    , dimension(n), intent(out):: p
    integer                  , intent(out):: sz
    type(graph)              , intent(in) :: g
    integer                  , intent(in) :: ii, n

    integer :: j1, j2

#if DBG>1
    if ( valid(g)<0 )     call quit( &
         & "graph_mod: graph_getRow_1: graph 'g' not valid")

    if ( (11<=0) .OR. (ii>g%nl) ) call quit( &
         & "graph_mod: graph_getRow_1: wrong line index 'ii'")
#endif

    j1 = g%row(ii)
    j2 = g%row(ii+1)
    sz = j2 - j1

#if DBG>1
    if (sz>n )  call quit( &
         & "graph_mod: graph_getRow_1: unsufficient size&
         &  for the array 'p'")
#endif

    p(1:sz) = g%col(j1:j2-1)
    
  end subroutine graph_getrow_1


  subroutine graph_getrow_2(p, n, g, ii)
    integer    , dimension(n), intent(out):: p
    type(graph)              , intent(in) :: g
    integer                  , intent(in) :: ii, n

    integer :: j1, j2

#if DBG>1
    if ( valid(g)<0 )     call quit( &
         & "graph_mod: graph_getRow_2: graph 'g' not valid")

    if ( (11<=0) .OR. (ii>g%nl) ) call quit( &
         & "graph_mod: graph_getRow_2: wrong line index 'ii'")
#endif

    j1 = g%row(ii)
    j2 = g%row(ii+1)

#if DBG>1
    if ( (j2-j1) /= n )  call quit( &
         & 'graph_mod: graph_getRow: 3')
#endif

    p  = g%col(j1:j2-1)
    
  end subroutine graph_getrow_2



  subroutine graph_setrow(g, col, p, ii)
    type(graph)            , intent(inout):: g
    integer                , intent(in)   :: p, ii
    integer  , dimension(p), intent(in)   :: col

    integer :: j1


#if DBG>1
    j1 = valid(g)

    if ( j1<0 ) then
       write(*,*) "graph_mod: graph_setRow: graph 'g' not valid, info = ", j1
       call quit("")
    end if
    if ( (11<=0) .OR. (ii>g%nl) ) call quit( &
         & "graph_mod: graph_setRow: wrong line index 'ii'")

    j1 = g%row(ii+1) - g%row(ii)
    if (j1 /= p ) call quit( 'graph_mod: graph_setRow: 3')
#endif

    j1 = g%row(ii)

    g%col(j1:j1+p-1) = col

  end subroutine graph_setrow



  subroutine graph_transpose(g2, g1)
    type(graph), intent(inout)  :: g2
    type(graph), intent(in)     :: g1

    integer, dimension(:), allocatable :: nnz
    integer :: jj, ii, j1, j2, col, kk, nc

    if ( valid(g1)<0 ) call quit( &
         & "graph_mod: graph_transpose: graph 'g' not valid")

    call clear(g2)

    !! Number of columns for g1
    !!
    nc = g1%nc
    ! if (nc==0) then
    !    nc = maxVal(g1%col)
    ! end if


    call allocmem(nnz, g1%nc)
    nnz = 0

    do ii=1, g1%nnz
       jj = g1%col(ii)
       nnz(jj) = nnz(jj) + 1
    end do
    g2 = graph(nnz, g1%nc, g1%nl)

    nnz = 0

    do ii=1, g1%nl

       j1 = g1%row(ii)
       j2 = g1%row(ii+1)-1

       do jj=j1, j2

          col = g1%col(jj)
          kk  = g2%row(col)  + nnz(col)
          g2%col(kk) = ii
          nnz(col) = nnz(col) + 1

       end do

    end do

    ! The graph is sorted by construction
    !
    g2%sorted = .true.

  end subroutine graph_transpose


  subroutine union(sz, tab, bf1, n, bf2, p, g, rows, m)
    integer              , intent(out):: sz
    integer              , intent(in) :: m, n, p
    integer, dimension(n), intent(out):: tab, bf1
    integer, dimension(p), intent(out):: bf2
    type(graph)          , intent(in) :: g
    integer, dimension(m), intent(in) :: rows
    
    integer :: ii, s1, s2

#if DBG>1
       if( maxval(rows)  > g%nl ) call quit( 'graph_mod: union: 1')
       if( minval(rows) <= 0    ) call quit( 'graph_mod: union: 2')
#endif

    call getrow(sz, tab, n, g, rows(1))

    do ii=2, m

       call getrow(s2, bf2, p, g, rows(ii))

       call merge_sorted_set(s1, bf1, n, bf2(1:s2), s2, tab(1:sz), sz)

#if DBG>1
       if( s1 > n ) call quit( 'graph_mod: union: 3')
#endif

       sz        = s1
       tab(1:s1) = bf1(1:s1)

    end do

  end subroutine union




  subroutine graph_extract(g2, g1, line)
    type(graph)          , intent(inout):: g2
    type(graph)          , intent(in)   :: g1
    logical, dimension(:), intent(in)   :: line

    integer, dimension(g1%nl) :: nnz
    integer :: ii, j1_1, j1_2, j2_1, j2_2

    if ( valid(g1)<0 ) call quit( &
         & 'graph_mod: graph_extract: 1')

    if (size(line,1)/=g1%nl) call quit( &
         & 'graph_mod: graph_extract: 2')

    call clear(g2)
    do ii=1, g1%nl
       if (line(ii)) then
          nnz(ii) = g1%row(ii+1)-g1%row(ii)
       else
          nnz(ii) = 0
       end if
    end do
    g2 = graph(nnz, g1%nl, g1%nc)

    !$OMP PARALLEL PRIVATE(j1_1, j1_2, j2_1, j2_2) 
    !$OMP DO
    do ii=1, g1%nl
       if (nnz(ii)==0) cycle

       j1_1 = g1%row(ii)
       j1_2 = g1%row(ii+1)

       j2_1 = g2%row(ii)
       j2_2 = g2%row(ii+1)

       g2%col(j2_1:j2_2-1) = g1%col(j1_1:j1_2-1)
       
    end do
    !$OMP END DO
    !$OMP END PARALLEL
    
    if (valid(g2)<0) call quit( &
         & 'graph_mod: graph_extract: 3')

  end subroutine graph_extract

  
  subroutine graph_sortrows(g)
    type(graph), intent(inout)  :: g

    integer :: ii, j1, j2, ll

    if ( valid(g)<0 ) call quit( &
         & 'graph_mod: graph_sortRows: 1')

    !$OMP PARALLEL PRIVATE(j1, j2, ll) 
    !$OMP DO
    do ii=1, g%nl
       
       j1 = g%row(ii)
       j2 = g%row(ii+1)
       ll = j2-j1
       if (ll<=1) cycle

       call sort(g%col(j1:j2-1), ll)

    end do
    !$OMP END DO
    !$OMP END PARALLEL

    g%sorted = .true.
    
  end subroutine graph_sortrows


  subroutine graph_prod(g, A, B)
    type(graph), intent(inout) :: g
    type(graph), intent(in)    :: A, B

    type(graph) :: B2
    integer     :: wA, wB, wC, nl

    if ( valid(a)<0 )   call quit( &
         & 'graph_mod: graph_prod: 1')

    if ( valid(b)<0 )   call quit( &
         & 'graph_mod: graph_prod: 2')

    if ( maxval(a%col) > b%nl ) call quit( &
         & 'graph_mod: graph_prod: 3')

    call clear(g)
    call copy(b2, b)
    if (.NOT. b2%sorted) call sortrows(b2)
    
    wa = a%width
    wb = b%width
    wc = wa * wb
    nl = a%nl

    ! graph construction
    !
    call loop()
    
    ! The graph is sorted by construction
    !
    g%sorted = .true.

  contains

    subroutine loop()

      integer, dimension(wc,nl) :: list
      integer, dimension(nl)    :: nnz
      integer, dimension(wA)    :: row_A
      integer, dimension(wB)    :: bf2
      integer, dimension(wC)    :: tab, bf1

      integer :: ii, jj, szA, sz

      ! Determines for each row ii :
      !            its length  = nnz(ii)
      !            its entries = list(1:nnz(ii), ii)
      !
      !$OMP PARALLEL PRIVATE(szA, sz, row_A, tab, bf1, bf2) 
      !$OMP DO
      do ii=1, nl
         call getrow(sza, row_a, wa, a, ii)
         if (sza==0) then
            nnz(ii) = 0
            cycle
         end if

         call union(sz, tab, bf1, wc, bf2, wb, b2, row_a(1:sza), sza)
        
         nnz(ii) = sz
         list(1:sz,ii) = tab(1:sz)
         
      end do
      !$OMP END DO
      !$OMP END PARALLEL
      
      ! clear temporary variables
      !
      call clear(b2)
      
      ! Graph allocation
      !
      g = graph(nnz, a%nl, b%nc)
      
      ! Fill in the graph
      !
      !$OMP PARALLEL PRIVATE(jj) 
      !$OMP DO
      do ii=1, nl
         jj = nnz(ii)
         call setrow(g, list(1:jj,ii), jj, ii)
      end do
      !$OMP END DO
      !$OMP END PARALLEL

    end subroutine loop
    
  end subroutine graph_prod


  subroutine prod_tab(g, A, B)
    type(graph), intent(inout) :: g
    type(graph), intent(in)    :: A, B

    type(graph) :: C

    if ( valid(a)<0 ) call quit( &
         & 'graph_mod: prod_tAB: 1')

    call transpose(c, a)
    call graph_prod(g, c, b)

  end subroutine prod_tab


  
end module graph_mod