Choral
graph_mod.F90
Go to the documentation of this file.
1 
32 
33 module graph_mod
34 
36  use basic_tools
37  use algebra_set
38  !$ use OMP_LIB
39 
40  implicit none
41  private
42 
43 
44  ! %----------------------------------------%
45  ! | |
46  ! | PUBLIC DATA |
47  ! | |
48  ! %----------------------------------------%
49  !! TYPE
50  public :: graph
51 
52  !! GENERIC
53  public :: clear ! TESTED
54  public :: valid ! TESTED
55  public :: copy ! TESTED
56  public :: equal ! TESTED
57  public :: sortrows ! TESTED
58  public :: getrow ! TESTED
59  public :: setrow ! TESTED
60  public :: transpose ! TESTED
61  public :: print
62 
63  !! NON GENERIC
64  public :: prod_tab ! TESTED
65  public :: graph_extract ! TESTED
66  public :: union ! TESTED
67  public :: graph_prod ! TESTED
68 
69 
70  ! %----------------------------------------%
71  ! | |
72  ! | DERIVED TYPE |
73  ! | |
74  ! %----------------------------------------%
78  type graph
80  integer :: nl = 0
82  integer :: nc = 0
84  integer :: nnz = 0
87  integer, dimension(:), allocatable :: row
92  integer, dimension(:), allocatable :: col
93 
95  integer :: width = 0
96 
98  logical :: sorted = .false.
99 
100  contains
101 
103  final :: graph_clear
104 
105  end type graph
106 
107  ! %----------------------------------------%
108  ! | |
109  ! | GENERIc SUBROUTINES |
110  ! | |
111  ! %----------------------------------------%
112  interface clear
113  module procedure graph_clear ! TESTED
114  end interface clear
115 
116  interface graph
117  module procedure graph_create ! TESTED
118  end interface graph
119 
120  interface valid
121  module procedure graph_valid ! TESTED
122  end interface valid
123 
124  interface copy
125  module procedure graph_copy ! TESTED
126  end interface copy
127 
129  interface print
130  module procedure graph_print
131  end interface print
132 
134  interface equal
135  module procedure graph_equal ! TESTED
136  end interface equal
137 
139  interface sortrows
140  module procedure graph_sortrows ! TESTED
141  end interface sortrows
142 
144  interface setrow
145  module procedure graph_setrow ! TESTED
146  end interface setrow
147 
149  interface getrow
150  module procedure graph_getrow_1 ! TESTED
151  module procedure graph_getrow_2 ! TESTED
152  end interface getrow
153 
155  interface transpose
156  module procedure graph_transpose ! TESTED
157  end interface transpose
158 
159 contains
160 
163  subroutine graph_clear(g)
164  type(graph), intent(inout) :: g
165 
166  call freemem(g%row)
167  call freemem(g%col)
168  g%nl = 0
169  g%nc = 0
170  g%nnz = 0
171  g%width = 0
172  g%sorted = .false.
173 
174  end subroutine graph_clear
175 
176 
187  function graph_create(nnz, nl, nc) result(g)
188  type(graph) :: g
189  integer, dimension(nl), intent(in) :: nnz
190  integer , intent(in) :: nl, nc
191 
192  integer :: nz, ii, jj, wdt, length
193 
194  if ( nc < 0 ) call quit( &
195  & "graph_mod: graph_create: number of columns nc < 0")
196 
197  if ( nl <= 0 ) call quit( &
198  & "graph_mod: graph_create: number of lines nl <= 0")
199 
200  if ( size(nnz,1) /= nl ) call quit( &
201  & "graph_mod: graph_create: 'nnz'wrong size")
202 
203  if ( .NOT. all( nnz >=0 ) ) call quit( &
204  & "graph_mod: graph_create: 'nnz(ii)' must be >=0 ")
205 
206  if ( sum( nnz ) == 0 ) call quit( &
207  & "graph_mod: graph_create: one line at least must be non-empty'")
208 
209  call clear(g)
210  g%nl = nl
211  g%nc = nc
212  nz = sum(nnz)
213  g%nnz = nz
214  call allocmem(g%row, g%nl+1)
215  call allocmem(g%col, nz )
216 
217  g%row = 0
218  g%col = 0
219  wdt = 0
220 
221  ! Construction of g%row
222  jj = 1
223  do ii = 1, g%nl
224  g%row(ii) = jj
225  length = nnz(ii)
226  if (length>wdt) wdt = length
227  jj = jj + length
228  end do
229  g%row(g%nl+1) = jj
230 
231  g%width = wdt
232 
233  end function graph_create
234 
235 
242  function graph_valid(g) result(info)
243  type(graph), intent(in) :: g
244  integer :: info
245 
246  integer :: ii, wdt, length
247 
248  info = 1
249 
250  if ( .NOT. (allocated(g%row)) ) info = -1
251  if ( .NOT. (allocated(g%col)) ) info = -2
252 
253  if (info > 0) then
254  if ( g%nl <= 0 ) info = -3
255  if ( g%nc <= 0 ) info = -4
256  if ( g%nnz <= 0 ) info = -5
257  if ( minval(g%row) <= 0 ) info = -6
258  if ( minval(g%col) < 0 ) info = -7
259  end if
260 
261  if (info > 0) then
262  if ( maxval(g%row) /= g%nnz+1 ) info = -8
263  if ( maxval(g%col) > g%nc ) info = -9
264  if ( size(g%row,1) /= g%nl+1 ) info = -10
265  if ( size(g%col,1) /= g%nnz ) info = -11
266  end if
267 
268  if (info > 0) then
269  if ( g%row(g%nl+1) /= size(g%col)+1 ) info = -12
270  end if
271 
272  if (info > 0) then
273  wdt = 0
274  do ii = 1, g%nl
275  length = g%row(ii+1) - g%row(ii)
276  if (length>wdt) wdt = length
277  end do
278  if ( wdt /= g%width ) info = -13
279  end if
280 
281  if (choral_verb>2) then
282  write(*,*)"graph_mod : valid, info =", int(info, 1)
283  end if
284 
285  end function graph_valid
286 
287 
290  subroutine graph_print(g)
291  type(graph), intent(in) :: g
292 
293  write(*,*) 'graph_mod : print'
294 
295  write(*,*) ' Valid =', int( valid(g), 1)
296  write(*,*) ' Nb lines =', g%nl
297  write(*,*) ' Nb columns =', g%nc
298  write(*,*) ' Nb non-zero entries =', g%nnz
299  write(*,*) ' Maximal line width =', g%width
300  write(*,*) ' Sorted lines = ', g%sorted
301 
302  end subroutine graph_print
303 
306  subroutine graph_copy(g, g0)
307  type(graph) , intent(inout) :: g
308  type(graph), intent(in) :: g0
309 
310  call clear(g)
311  g%nl = g0%nl
312  g%nc = g0%nc
313  g%nnz = g0%nnz
314  g%width = g0%width
315  g%sorted = g0%sorted
316 
317  call allocmem(g%row, size(g0%row,1) )
318  call allocmem(g%col, size(g0%col,1) )
319  g%row = g0%row
320  g%col = g0%col
321 
322  end subroutine graph_copy
323 
324 
330  function graph_equal(g2, g1) result(b)
331  logical :: b
332  type(graph), intent(in):: g1, g2
333 
334  integer, dimension(:), allocatable :: p1, p2
335  integer :: ii, wdt, sz
336 
337  b = ( valid(g1)>0 ) .AND. ( valid(g2)>0 )
338  if(.NOT.b) return
339 
340  b = b .AND. ( g1%nl == g2%nl ) .AND. ( g1%nnz == g2%nnz )
341  b = b .AND. ( g1%nc == g2%nc )
342  b = b .AND. ( g1%width == g2%width )
343  if(.NOT.b) return
344 
345  b = b .AND. (all( g1%row == g2%row ) )
346  if(.NOT.b) return
347 
348  wdt = g1%width
349  call allocmem( p1, wdt)
350  call allocmem( p2, wdt )
351 
352  do ii=1, g1%nl
353  call getrow(sz, p1, wdt, g1, ii)
354  if( sz == 0 ) cycle
355 
356  call getrow(sz, p2, wdt, g2, ii)
357 
358  call sort( p1(1:sz), sz )
359  call sort( p2(1:sz), sz )
360 
361  b = b .AND. (all( p1(1:sz) == p2(1:sz) ) )
362 
363  end do
364 
365  end function graph_equal
366 
367 
368 
379  subroutine graph_getrow_1(sz, p, n, g, ii)
380  integer , dimension(n), intent(out):: p
381  integer , intent(out):: sz
382  type(graph) , intent(in) :: g
383  integer , intent(in) :: ii, n
384 
385  integer :: j1, j2
386 
387 #if DBG>1
388  if ( valid(g)<0 ) call quit( &
389  & "graph_mod: graph_getRow_1: graph 'g' not valid")
390 
391  if ( (11<=0) .OR. (ii>g%nl) ) call quit( &
392  & "graph_mod: graph_getRow_1: wrong line index 'ii'")
393 #endif
394 
395  j1 = g%row(ii)
396  j2 = g%row(ii+1)
397  sz = j2 - j1
398 
399 #if DBG>1
400  if (sz>n ) call quit( &
401  & "graph_mod: graph_getRow_1: unsufficient size&
402  & for the array 'p'")
403 #endif
404 
405  p(1:sz) = g%col(j1:j2-1)
406 
407  end subroutine graph_getrow_1
408 
409 
420  subroutine graph_getrow_2(p, n, g, ii)
421  integer , dimension(n), intent(out):: p
422  type(graph) , intent(in) :: g
423  integer , intent(in) :: ii, n
424 
425  integer :: j1, j2
426 
427 #if DBG>1
428  if ( valid(g)<0 ) call quit( &
429  & "graph_mod: graph_getRow_2: graph 'g' not valid")
430 
431  if ( (11<=0) .OR. (ii>g%nl) ) call quit( &
432  & "graph_mod: graph_getRow_2: wrong line index 'ii'")
433 #endif
434 
435  j1 = g%row(ii)
436  j2 = g%row(ii+1)
437 
438 #if DBG>1
439  if ( (j2-j1) /= n ) call quit( &
440  & 'graph_mod: graph_getRow: 3')
441 #endif
442 
443  p = g%col(j1:j2-1)
444 
445  end subroutine graph_getrow_2
446 
447 
448 
449  subroutine graph_setrow(g, col, p, ii)
450  type(graph) , intent(inout):: g
451  integer , intent(in) :: p, ii
452  integer , dimension(p), intent(in) :: col
453 
454  integer :: j1
455 
456 
457 #if DBG>1
458  j1 = valid(g)
459 
460  if ( j1<0 ) then
461  write(*,*) "graph_mod: graph_setRow: graph 'g' not valid, info = ", j1
462  call quit("")
463  end if
464  if ( (11<=0) .OR. (ii>g%nl) ) call quit( &
465  & "graph_mod: graph_setRow: wrong line index 'ii'")
466 
467  j1 = g%row(ii+1) - g%row(ii)
468  if (j1 /= p ) call quit( 'graph_mod: graph_setRow: 3')
469 #endif
470 
471  j1 = g%row(ii)
472 
473  g%col(j1:j1+p-1) = col
474 
475  end subroutine graph_setrow
476 
477 
478 
484  subroutine graph_transpose(g2, g1)
485  type(graph), intent(inout) :: g2
486  type(graph), intent(in) :: g1
487 
488  integer, dimension(:), allocatable :: nnz
489  integer :: jj, ii, j1, j2, col, kk, nc
490 
491  if ( valid(g1)<0 ) call quit( &
492  & "graph_mod: graph_transpose: graph 'g' not valid")
493 
494  call clear(g2)
495 
496  !! Number of columns for g1
497  !!
498  nc = g1%nc
499  ! if (nc==0) then
500  ! nc = maxVal(g1%col)
501  ! end if
502 
503 
504  call allocmem(nnz, g1%nc)
505  nnz = 0
506 
507  do ii=1, g1%nnz
508  jj = g1%col(ii)
509  nnz(jj) = nnz(jj) + 1
510  end do
511  g2 = graph(nnz, g1%nc, g1%nl)
512 
513  nnz = 0
514 
515  do ii=1, g1%nl
516 
517  j1 = g1%row(ii)
518  j2 = g1%row(ii+1)-1
519 
520  do jj=j1, j2
521 
522  col = g1%col(jj)
523  kk = g2%row(col) + nnz(col)
524  g2%col(kk) = ii
525  nnz(col) = nnz(col) + 1
526 
527  end do
528 
529  end do
530 
531  ! The graph is sorted by construction
532  !
533  g2%sorted = .true.
534 
535  end subroutine graph_transpose
536 
537 
557  subroutine union(sz, tab, bf1, n, bf2, p, g, rows, m)
558  integer , intent(out):: sz
559  integer , intent(in) :: m, n, p
560  integer, dimension(n), intent(out):: tab, bf1
561  integer, dimension(p), intent(out):: bf2
562  type(graph) , intent(in) :: g
563  integer, dimension(m), intent(in) :: rows
564 
565  integer :: ii, s1, s2
566 
567 #if DBG>1
568  if( maxval(rows) > g%nl ) call quit( 'graph_mod: union: 1')
569  if( minval(rows) <= 0 ) call quit( 'graph_mod: union: 2')
570 #endif
571 
572  call getrow(sz, tab, n, g, rows(1))
573 
574  do ii=2, m
575 
576  call getrow(s2, bf2, p, g, rows(ii))
577 
578  call merge_sorted_set(s1, bf1, n, bf2(1:s2), s2, tab(1:sz), sz)
579 
580 #if DBG>1
581  if( s1 > n ) call quit( 'graph_mod: union: 3')
582 #endif
583 
584  sz = s1
585  tab(1:s1) = bf1(1:s1)
586 
587  end do
588 
589  end subroutine union
590 
591 
592 
593 
597  subroutine graph_extract(g2, g1, line)
598  type(graph) , intent(inout):: g2
599  type(graph) , intent(in) :: g1
600  logical, dimension(:), intent(in) :: line
601 
602  integer, dimension(g1%nl) :: nnz
603  integer :: ii, j1_1, j1_2, j2_1, j2_2
604 
605  if ( valid(g1)<0 ) call quit( &
606  & 'graph_mod: graph_extract: 1')
607 
608  if (size(line,1)/=g1%nl) call quit( &
609  & 'graph_mod: graph_extract: 2')
610 
611  call clear(g2)
612  do ii=1, g1%nl
613  if (line(ii)) then
614  nnz(ii) = g1%row(ii+1)-g1%row(ii)
615  else
616  nnz(ii) = 0
617  end if
618  end do
619  g2 = graph(nnz, g1%nl, g1%nc)
620 
621  !$OMP PARALLEL PRIVATE(j1_1, j1_2, j2_1, j2_2)
622  !$OMP DO
623  do ii=1, g1%nl
624  if (nnz(ii)==0) cycle
625 
626  j1_1 = g1%row(ii)
627  j1_2 = g1%row(ii+1)
628 
629  j2_1 = g2%row(ii)
630  j2_2 = g2%row(ii+1)
631 
632  g2%col(j2_1:j2_2-1) = g1%col(j1_1:j1_2-1)
633 
634  end do
635  !$OMP END DO
636  !$OMP END PARALLEL
637 
638  if (valid(g2)<0) call quit( &
639  & 'graph_mod: graph_extract: 3')
640 
641  end subroutine graph_extract
642 
643 
647  subroutine graph_sortrows(g)
648  type(graph), intent(inout) :: g
649 
650  integer :: ii, j1, j2, ll
651 
652  if ( valid(g)<0 ) call quit( &
653  & 'graph_mod: graph_sortRows: 1')
654 
655  !$OMP PARALLEL PRIVATE(j1, j2, ll)
656  !$OMP DO
657  do ii=1, g%nl
658 
659  j1 = g%row(ii)
660  j2 = g%row(ii+1)
661  ll = j2-j1
662  if (ll<=1) cycle
663 
664  call sort(g%col(j1:j2-1), ll)
665 
666  end do
667  !$OMP END DO
668  !$OMP END PARALLEL
669 
670  g%sorted = .true.
671 
672  end subroutine graph_sortrows
673 
674 
678  subroutine graph_prod(g, A, B)
679  type(graph), intent(inout) :: g
680  type(graph), intent(in) :: A, B
681 
682  type(graph) :: B2
683  integer :: wA, wB, wC, nl
684 
685  if ( valid(a)<0 ) call quit( &
686  & 'graph_mod: graph_prod: 1')
687 
688  if ( valid(b)<0 ) call quit( &
689  & 'graph_mod: graph_prod: 2')
690 
691  if ( maxval(a%col) > b%nl ) call quit( &
692  & 'graph_mod: graph_prod: 3')
693 
694  call clear(g)
695  call copy(b2, b)
696  if (.NOT. b2%sorted) call sortrows(b2)
697 
698  wa = a%width
699  wb = b%width
700  wc = wa * wb
701  nl = a%nl
702 
703  ! graph construction
704  !
705  call loop()
706 
707  ! The graph is sorted by construction
708  !
709  g%sorted = .true.
710 
711  contains
712 
713  subroutine loop()
715  integer, dimension(wc,nl) :: list
716  integer, dimension(nl) :: nnz
717  integer, dimension(wA) :: row_A
718  integer, dimension(wB) :: bf2
719  integer, dimension(wC) :: tab, bf1
720 
721  integer :: ii, jj, szA, sz
722 
723  ! Determines for each row ii :
724  ! its length = nnz(ii)
725  ! its entries = list(1:nnz(ii), ii)
726  !
727  !$OMP PARALLEL PRIVATE(szA, sz, row_A, tab, bf1, bf2)
728  !$OMP DO
729  do ii=1, nl
730  call getrow(sza, row_a, wa, a, ii)
731  if (sza==0) then
732  nnz(ii) = 0
733  cycle
734  end if
735 
736  call union(sz, tab, bf1, wc, bf2, wb, b2, row_a(1:sza), sza)
737 
738  nnz(ii) = sz
739  list(1:sz,ii) = tab(1:sz)
740 
741  end do
742  !$OMP END DO
743  !$OMP END PARALLEL
744 
745  ! clear temporary variables
746  !
747  call clear(b2)
748 
749  ! Graph allocation
750  !
751  g = graph(nnz, a%nl, b%nc)
752 
753  ! Fill in the graph
754  !
755  !$OMP PARALLEL PRIVATE(jj)
756  !$OMP DO
757  do ii=1, nl
758  jj = nnz(ii)
759  call setrow(g, list(1:jj,ii), jj, ii)
760  end do
761  !$OMP END DO
762  !$OMP END PARALLEL
763 
764  end subroutine loop
765 
766  end subroutine graph_prod
767 
768 
773  subroutine prod_tab(g, A, B)
774  type(graph), intent(inout) :: g
775  type(graph), intent(in) :: A, B
776 
777  type(graph) :: C
778 
779  if ( valid(a)<0 ) call quit( &
780  & 'graph_mod: prod_tAB: 1')
781 
782  call transpose(c, a)
783  call graph_prod(g, c, b)
784 
785  end subroutine prod_tab
786 
787 
788 
789 end module graph_mod
790 
set a graph row
Definition: graph_mod.F90:144
subroutine graph_sortrows(g)
Sort all rows of the graph g.
Definition: graph_mod.F90:648
DERIVED TYPE graph: sparse matrices of boolean in CSR format
Definition: graph_mod.F90:33
subroutine loop()
BASIC TOOLS
Definition: basic_tools.F90:8
subroutine graph_copy(g, g0)
g := g0
Definition: graph_mod.F90:307
deallocate memory for real(RP) arrays
Definition: real_type.F90:163
type(graph) function graph_create(nnz, nl, nc)
Constructor for the type graph
Definition: graph_mod.F90:188
subroutine graph_transpose(g2, g1)
Transpose graph.
Definition: graph_mod.F90:485
subroutine graph_print(g)
Print a short description of the graph &#39;g&#39;.
Definition: graph_mod.F90:291
print a short description
Definition: graph_mod.F90:129
subroutine, public graph_prod(g, A, B)
Matrix product g = A*B.
Definition: graph_mod.F90:679
subroutine, public quit(message)
Stop program execution, display an error messahe.
subroutine, public graph_extract(g2, g1, line)
line i of g2 = line i of g1 if line(i) = .TRUE. = void otherwise
Definition: graph_mod.F90:598
subroutine graph_setrow(g, col, p, ii)
Definition: graph_mod.F90:450
subroutine graph_getrow_1(sz, p, n, g, ii)
Extract row ii of the graph g.
Definition: graph_mod.F90:380
transpose graph
Definition: graph_mod.F90:155
subroutine graph_clear(g)
Destructor for the type graph.
Definition: graph_mod.F90:164
integer function graph_valid(g)
Check if the graph structure &#39;g&#39; is valid.
Definition: graph_mod.F90:243
subroutine, public merge_sorted_set(cpt, t, n, t1, n1, t2, n2)
t(1:cpt) = merge of the two arrays t1 and t2 of size n1 and n2
Definition: algebra_set.F90:37
ALGEBRAIC OPERATIONS ON SETS
Definition: algebra_set.F90:14
subroutine, public union(sz, tab, bf1, n, bf2, p, g, rows, m)
graph row union
Definition: graph_mod.F90:558
allocate memory for real(RP) arrays
Definition: real_type.F90:158
subroutine graph_getrow_2(p, n, g, ii)
Extract row ii of the graph g.
Definition: graph_mod.F90:421
logical function graph_equal(g2, g1)
test equality Two graph are equal if the underlying boolean matrices are equal independently of the c...
Definition: graph_mod.F90:331
integer choral_verb
Verbosity level.
subroutine, public sort(t, n)
Sorts the array t of length N in ascending order by the straight insertion method.
DEFINITION OF GLOBAL VARIABLES FOR THE LIBRARY CHORAL
sort graph rows in ascending order
Definition: graph_mod.F90:139
check graph equality
Definition: graph_mod.F90:134
The type graph stores sparse matrices of boolean in CSR format.
Definition: graph_mod.F90:78
subroutine, public prod_tab(g, A, B)
graph product
Definition: graph_mod.F90:774
extract a graph row
Definition: graph_mod.F90:149