MODFLOW 6  version 6.6.0.dev0
USGS Modular Hydrologic Model
dag_module.f90
Go to the documentation of this file.
1 !*******************************************************************************
2 !>
3 ! DAG Module.
4 
5  module dag_module
6 
7  implicit none
8 
9  private
10 
11  type :: vertex
12  !! a vertex of a directed acyclic graph (DAG)
13  private
14  integer,dimension(:),allocatable :: edges !! these are the vertices that this vertex depends on
15  integer :: ivertex = 0 !! vertex number
16  logical :: checking = .false. !! used for toposort
17  logical :: marked = .false. !! used for toposort
18  character(len=:),allocatable :: label !! used for diagraph
19  character(len=:),allocatable :: attributes !! used for diagraph
20  contains
21  generic :: set_edges => set_edge_vector, add_edge
22  procedure :: set_edge_vector
23  procedure :: add_edge
24  end type vertex
25 
26  type,public :: dag
27  !! a directed acyclic graph (DAG)
28  private
29  integer :: n = 0 !! number of `vertices`
30  type(vertex),dimension(:),allocatable :: vertices !! the vertices in the DAG.
31  contains
32  procedure,public :: set_vertices => dag_set_vertices
33  procedure,public :: set_edges => dag_set_edges
34  procedure,public :: set_vertex_info => dag_set_vertex_info
35  procedure,public :: toposort => dag_toposort
36  procedure,public :: generate_digraph => dag_generate_digraph
37  procedure,public :: generate_dependency_matrix => dag_generate_dependency_matrix
38  procedure,public :: save_digraph => dag_save_digraph
39  procedure,public :: get_edges => dag_get_edges
40  procedure,public :: get_dependencies => dag_get_dependencies
41  procedure,public :: destroy => dag_destroy
42  end type dag
43 
44  contains
45 !*******************************************************************************
46 
47 !*******************************************************************************
48 !>
49 ! Destroy the `dag`.
50 
51  subroutine dag_destroy(me)
52 
53  class(dag),intent(inout) :: me
54 
55  me%n = 0
56  if (allocated(me%vertices)) deallocate(me%vertices)
57 
58  end subroutine dag_destroy
59 !*******************************************************************************
60 
61 !*******************************************************************************
62 !>
63 ! specify the edge indices for this vertex
64 
65  subroutine set_edge_vector(me,edges)
66 
67  class(vertex),intent(inout) :: me
68  integer,dimension(:),intent(in) :: edges
69 
70  integer :: i !! counter
71 
72  if (allocated(me%edges)) then
73  do i=1,size(edges)
74  call me%add_edge(edges(i))
75  end do
76  else
77  allocate(me%edges(size(edges))) ! note: not checking for uniqueness here.
78  me%edges = edges
79  end if
80 
81  end subroutine set_edge_vector
82 !*******************************************************************************
83 
84 !*******************************************************************************
85 !>
86 ! add an edge index for this vertex
87 
88  subroutine add_edge(me,edge)
89 
90  class(vertex),intent(inout) :: me
91  integer,intent(in) :: edge
92 
93  if (allocated(me%edges)) then
94  if (.not. any(edge==me%edges)) then
95  me%edges = [me%edges, edge] ! auto lhs reallocation
96  end if
97  else
98  allocate(me%edges(1))
99  me%edges = [edge]
100  end if
101 
102  end subroutine add_edge
103 !*******************************************************************************
104 
105 !*******************************************************************************
106 !>
107 ! get the edges for the vertex (all the the vertices
108 ! that this vertex depends on).
109 
110  pure function dag_get_edges(me,ivertex) result(edges)
111 
112  implicit none
113 
114  class(dag),intent(in) :: me
115  integer,intent(in) :: ivertex
116  integer,dimension(:),allocatable :: edges
117 
118  if (ivertex>0 .and. ivertex <= me%n) then
119  edges = me%vertices(ivertex)%edges ! auto LHS allocation
120  end if
121 
122  end function dag_get_edges
123 !*******************************************************************************
124 
125 !*******************************************************************************
126 !>
127 ! get all the vertices that depend on this vertex.
128 
129  pure function dag_get_dependencies(me,ivertex) result(dep)
130 
131  implicit none
132 
133  class(dag),intent(in) :: me
134  integer,intent(in) :: ivertex
135  integer,dimension(:),allocatable :: dep !! the set of all vertices
136  !! than depend on `ivertex`
137 
138  integer :: i !! vertex counter
139 
140  if (ivertex>0 .and. ivertex <= me%n) then
141 
142  ! have to check all the vertices:
143  do i=1, me%n
144  if (allocated(me%vertices(i)%edges)) then
145  if (any(me%vertices(i)%edges == ivertex)) then
146  if (allocated(dep)) then
147  dep = [dep, i] ! auto LHS allocation
148  else
149  dep = [i] ! auto LHS allocation
150  end if
151  end if
152  end if
153  end do
154 
155  end if
156 
157  end function dag_get_dependencies
158 !*******************************************************************************
159 
160 !*******************************************************************************
161 !>
162 ! set the number of vertices in the dag
163 
164  subroutine dag_set_vertices(me,nvertices)
165 
166  class(dag),intent(inout) :: me
167  integer,intent(in) :: nvertices !! number of vertices
168 
169  integer :: i
170 
171  me%n = nvertices
172  allocate(me%vertices(nvertices))
173  me%vertices%ivertex = [(i,i=1,nvertices)]
174 
175  end subroutine dag_set_vertices
176 !*******************************************************************************
177 
178 !*******************************************************************************
179 !>
180 ! set info about a vertex in a dag.
181 
182  subroutine dag_set_vertex_info(me,ivertex,label,attributes)
183 
184  class(dag),intent(inout) :: me
185  integer,intent(in) :: ivertex !! vertex number
186  character(len=*),intent(in),optional :: label !! if a label is not set,
187  !! then the integer vertex
188  !! number is used.
189  character(len=*),intent(in),optional :: attributes !! other attributes when
190  !! saving as a diagraph.
191 
192  if (present(label)) then
193  me%vertices(ivertex)%label = label
194  else
195  ! just use the vertex number
196  me%vertices(ivertex)%label = integer_to_string(ivertex)
197  end if
198 
199  if (present(attributes)) then
200  me%vertices(ivertex)%attributes = attributes
201  end if
202 
203  end subroutine dag_set_vertex_info
204 !*******************************************************************************
205 
206 !*******************************************************************************
207 !>
208 ! set the edges for a vertex in a dag
209 
210  subroutine dag_set_edges(me,ivertex,edges)
211 
212  class(dag),intent(inout) :: me
213  integer,intent(in) :: ivertex !! vertex number
214  integer,dimension(:),intent(in) :: edges
215 
216  call me%vertices(ivertex)%set_edges(edges)
217 
218  end subroutine dag_set_edges
219 !*******************************************************************************
220 
221 !*******************************************************************************
222 !>
223 ! Main toposort routine
224 
225  subroutine dag_toposort(me,order,istat)
226 
227  class(dag),intent(inout) :: me
228  integer,dimension(:),allocatable,intent(out) :: order !! the toposort order
229  integer,intent(out) :: istat !! Status flag:
230  !!
231  !! * 0 if no errors
232  !! * -1 if circular dependency
233  !! (in this case, `order` will not be allocated)
234 
235  integer :: i,iorder
236 
237  if (me%n==0) return
238 
239  allocate(order(me%n))
240 
241  iorder = 0 ! index in order array
242  istat = 0 ! no errors so far
243  do i=1,me%n
244  if (.not. me%vertices(i)%marked) call dfs(me%vertices(i))
245  if (istat==-1) exit
246  end do
247 
248  if (istat==-1) deallocate(order)
249 
250  contains
251 
252  recursive subroutine dfs(v)
253 
254  !! depth-first graph traversal
255 
256  type(vertex),intent(inout) :: v
257  integer :: j
258 
259  if (istat==-1) return
260 
261  if (v%checking) then
262  ! error: circular dependency
263  istat = -1
264  else
265  if (.not. v%marked) then
266  v%checking = .true.
267  if (allocated(v%edges)) then
268  do j=1,size(v%edges)
269  call dfs(me%vertices(v%edges(j)))
270  if (istat==-1) return
271  end do
272  end if
273  v%checking = .false.
274  v%marked = .true.
275  iorder = iorder + 1
276  order(iorder) = v%ivertex
277  end if
278  end if
279 
280  end subroutine dfs
281 
282  end subroutine dag_toposort
283 !*******************************************************************************
284 
285 !*******************************************************************************
286 !>
287 ! Generate a Graphviz digraph structure for the DAG.
288 !
289 !### Example
290 ! * To convert this to a PDF using `dot`: `dot -Tpdf -o test.pdf test.dot`,
291 ! where `test.dot` is `str` written to a file.
292 
293  function dag_generate_digraph(me,rankdir,dpi,edge) result(str)
294 
295  implicit none
296 
297  class(dag),intent(in) :: me
298  character(len=:),allocatable :: str
299  character(len=*),intent(in),optional :: rankdir !! right to left orientation (e.g. 'RL')
300  integer,intent(in),optional :: dpi !! resolution (e.g. 300)
301  character(len=*),intent(in),optional :: edge !! right to left orientation (e.g. 'forward' or 'back)
302 
303  integer :: i,j !! counter
304  integer :: n_edges !! number of edges
305  character(len=:),allocatable :: attributes,label
306  logical :: has_label, has_attributes
307 
308  character(len=*),parameter :: tab = ' ' !! for indenting
309  character(len=*),parameter :: newline = new_line(' ') !! line break character
310 
311  if (me%n == 0) return
312 
313  str = 'digraph G {'//newline//newline
314  if (present(rankdir)) &
315  str = str//tab//'rankdir='//rankdir//newline//newline
316  if (present(dpi)) &
317  str = str//tab//'graph [ dpi = '//integer_to_string(dpi)//' ]'//newline//newline
318  if (present(edge)) &
319  str = str//tab//'edge [ dir = "'//trim(adjustl(edge))//'" ]'//newline//newline
320 
321  ! define the vertices:
322  do i=1,me%n
323  has_label = allocated(me%vertices(i)%label)
324  has_attributes = allocated(me%vertices(i)%attributes)
325  if (has_label) label = 'label="'//trim(adjustl(me%vertices(i)%label))//'"'
326  if (has_label .and. has_attributes) then
327  attributes = '['//trim(adjustl(me%vertices(i)%attributes))//','//label//']'
328  elseif (has_label .and. .not. has_attributes) then
329  attributes = '['//label//']'
330  elseif (.not. has_label .and. has_attributes) then
331  attributes = '['//trim(adjustl(me%vertices(i)%attributes))//']'
332  else ! neither
333  attributes = ''
334  end if
335  str = str//tab//integer_to_string(i)//' '//attributes//newline
336  if (i==me%n) str = str//newline
337  end do
338 
339  ! define the dependencies:
340  do i=1,me%n
341  if (allocated(me%vertices(i)%edges)) then
342  n_edges = size(me%vertices(i)%edges)
343  str = str//tab//integer_to_string(i)//' -> '
344  do j=1,n_edges
345  ! comma-separated list:
346  str = str//integer_to_string(me%vertices(i)%edges(j))
347  if (n_edges>1 .and. j<n_edges) str = str//','
348  end do
349  str = str//';'//newline
350  end if
351  end do
352 
353  str = str//newline//'}'
354 
355  end function dag_generate_digraph
356 !*******************************************************************************
357 
358 !*******************************************************************************
359 !>
360 ! Generate the dependency matrix for the DAG.
361 !
362 ! This is an \‍(n \times n \‍) matrix with elements \‍(A_{ij}\‍),
363 ! such that \‍(A_{ij}\‍) is true if vertex \‍(i\‍) depends on vertex \‍(j\‍).
364 
366 
367  implicit none
368 
369  class(dag),intent(in) :: me
370  logical,dimension(:,:),intent(out),allocatable :: mat !! dependency matrix
371 
372  integer :: i !! vertex counter
373 
374  if (me%n > 0) then
375 
376  allocate(mat(me%n,me%n))
377  mat = .false.
378 
379  do i=1,me%n
380  if (allocated(me%vertices(i)%edges)) then
381  mat(i,me%vertices(i)%edges) = .true.
382  end if
383  end do
384 
385  end if
386 
387  end subroutine dag_generate_dependency_matrix
388 !*******************************************************************************
389 
390 !*******************************************************************************
391 !>
392 ! Generate a Graphviz digraph structure for the DAG and write it to a file.
393 
394  subroutine dag_save_digraph(me,filename,rankdir,dpi,edge)
395 
396  implicit none
397 
398  class(dag),intent(in) :: me
399  character(len=*),intent(in),optional :: filename !! file name for diagraph
400  character(len=*),intent(in),optional :: rankdir !! right to left orientation (e.g. 'RL')
401  integer,intent(in),optional :: dpi !! resolution (e.g. 300)
402  character(len=*),intent(in),optional :: edge !! right to left orientation (e.g. 'forward' or 'back)
403 
404  integer :: iunit, istat
405  character(len=:),allocatable :: diagraph
406 
407  diagraph = me%generate_digraph(rankdir,dpi,edge)
408 
409  open(newunit=iunit,file=filename,status='REPLACE',iostat=istat)
410 
411  if (istat==0) then
412  write(iunit,fmt='(A)',iostat=istat) diagraph
413  else
414  write(*,*) 'error opening '//trim(filename)
415  end if
416 
417  close(iunit,iostat=istat)
418 
419  end subroutine dag_save_digraph
420 !*******************************************************************************
421 
422 !*******************************************************************************
423 !>
424 ! Integer to allocatable string.
425 
426  pure function integer_to_string(i) result(s)
427 
428  implicit none
429 
430  integer,intent(in) :: i
431  character(len=:),allocatable :: s
432 
433  integer :: istat
434 
435  allocate( character(len=64) :: s ) ! should be big enough
436  write(s,fmt='(ss,I0)',iostat=istat) i
437  if (istat==0) then
438  s = trim(adjustl(s))
439  else
440  s = '***'
441  end if
442 
443  end function integer_to_string
444 !*******************************************************************************
445 
446 !*******************************************************************************
447  end module dag_module
448 !*******************************************************************************
recursive subroutine dfs(v)
Definition: dag_module.f90:253
pure character(len=:) function, allocatable integer_to_string(i)
Definition: dag_module.f90:427
character(len=:) function, allocatable dag_generate_digraph(me, rankdir, dpi, edge)
Definition: dag_module.f90:294
subroutine add_edge(me, edge)
Definition: dag_module.f90:89
pure integer function, dimension(:), allocatable dag_get_dependencies(me, ivertex)
Definition: dag_module.f90:130
subroutine dag_generate_dependency_matrix(me, mat)
Definition: dag_module.f90:366
subroutine dag_set_vertex_info(me, ivertex, label, attributes)
Definition: dag_module.f90:183
subroutine dag_save_digraph(me, filename, rankdir, dpi, edge)
Definition: dag_module.f90:395
subroutine dag_destroy(me)
Definition: dag_module.f90:52
subroutine dag_toposort(me, order, istat)
Definition: dag_module.f90:226
pure integer function, dimension(:), allocatable dag_get_edges(me, ivertex)
Definition: dag_module.f90:111
subroutine set_edge_vector(me, edges)
Definition: dag_module.f90:66
subroutine dag_set_edges(me, ivertex, edges)
Definition: dag_module.f90:211
subroutine dag_set_vertices(me, nvertices)
Definition: dag_module.f90:165