1!-------------------------------------------------------------------------------
2! Copyright (c) 2020 FrontISTR Commons
3! This software is released under the MIT License, see LICENSE.txt
4!-------------------------------------------------------------------------------
5
6!----------------------------------------------------------------------
7!> @brief HECMW_ORDERING_QMD is a program for the minimum degree
8!         ordering using quotient graphs
9!----------------------------------------------------------------------
10module hecmw_ordering_qmd
11  use hecmw_util
12  implicit none
13
14  private
15  public :: hecmw_ordering_genqmd
16
17contains
18
19  !======================================================================!
20  !> @brief hecmw_ordering_GENQMD
21  !======================================================================!
22  subroutine hecmw_ordering_GENQMD(Neqns,Nttbr,Xadj,Adj0,Perm,Invp)
23    implicit none
24    !------
25    integer(kind=kint), intent(in):: Neqns
26    integer(kind=kint), intent(in):: Nttbr
27    integer(kind=kint), intent(in):: Adj0(:)
28    integer(kind=kint), intent(in):: Xadj(:)
29    integer(kind=kint), intent(out):: Perm(:)
30    integer(kind=kint), intent(out):: Invp(:)
31    !------
32    integer(kind=kint), allocatable:: Deg(:)
33    integer(kind=kint), allocatable:: Marker(:)
34    integer(kind=kint), allocatable:: Rchset(:)
35    integer(kind=kint), allocatable:: Nbrhd(:)
36    integer(kind=kint), allocatable:: Qsize(:)
37    integer(kind=kint), allocatable:: Qlink(:)
38    integer(kind=kint), allocatable:: Adjncy(:)
39    integer(kind=kint):: inode
40    integer(kind=kint):: ip
41    integer(kind=kint):: irch
42    integer(kind=kint):: j
43    integer(kind=kint):: mindeg
44    integer(kind=kint):: ndeg
45    integer(kind=kint):: nhdsze
46    integer(kind=kint):: node
47    integer(kind=kint):: np
48    integer(kind=kint):: num
49    integer(kind=kint):: nump1
50    integer(kind=kint):: nxnode
51    integer(kind=kint):: rchsze
52    integer(kind=kint):: search
53    integer(kind=kint):: thresh
54    integer(kind=kint):: ierror
55    logical:: found
56
57    allocate (DEG(NEQns+1),stat=IERror)
58    if ( IERror/=0 ) stop "ALLOCATION ERROR, deg: SUB. genqmd"
59    allocate (MARker(NEQns+1),stat=IERror)
60    if ( IERror/=0 ) stop "ALLOCATION ERROR, marker: SUB. genqmd"
61    allocate (RCHset(NEQns+2),stat=IERror)
62    if ( IERror/=0 ) stop "ALLOCATION ERROR, rchset: SUB. genqmd"
63    allocate (NBRhd(NEQns+1),stat=IERror)
64    if ( IERror/=0 ) stop "ALLOCATION ERROR, nbrhd: SUB. genqmd"
65    allocate (QSIze(NEQns+1),stat=IERror)
66    if ( IERror/=0 ) stop "ALLOCATION ERROR, qsize: SUB. genqmd"
67    allocate (QLInk(NEQns+1),stat=IERror)
68    if ( IERror/=0 ) stop "ALLOCATION ERROR, qlink: SUB. genqmd"
69    allocate (ADJncy(2*NTTbr),stat=IERror)
70    if ( IERror/=0 ) stop "ALLOCATION ERROR, adjncy: SUB. genqmd"
71
72    mindeg = Neqns
73    Adjncy(1:Xadj(Neqns+1) - 1) = Adj0(1:Xadj(Neqns+1) - 1)
74    do node = 1, Neqns
75      Perm(node) = node
76      Invp(node) = node
77      Marker(node) = 0
78      Qsize(node) = 1
79      Qlink(node) = 0
80      ndeg = Xadj(node+1) - Xadj(node)
81      Deg(node) = ndeg
82      if ( ndeg<mindeg ) mindeg = ndeg
83    enddo
84
85    num = 0
86    loop1: do
87      search = 1
88      thresh = mindeg
89      mindeg = Neqns
90      loop2: do
91        nump1 = num + 1
92        if ( nump1>search ) search = nump1
93        found = .false.
94        do j = search, Neqns
95          node = Perm(j)
96          if ( Marker(node)>=0 ) then
97            ndeg = Deg(node)
98            if ( ndeg<=thresh ) then
99              found = .true.
100              exit
101            endif
102            if ( ndeg<mindeg ) mindeg = ndeg
103          endif
104        enddo
105        if (.not. found) cycle loop1
106
107        search = j
108        Marker(node) = 1
109        call QMDRCH(node,Xadj,Adjncy,Deg,Marker,rchsze,Rchset,nhdsze,Nbrhd)
110        nxnode = node
111        do
112          num = num + 1
113          np = Invp(nxnode)
114          ip = Perm(num)
115          Perm(np) = ip
116          Invp(ip) = np
117          Perm(num) = nxnode
118          Invp(nxnode) = num
119          Deg(nxnode) = -1
120          nxnode = Qlink(nxnode)
121          if ( nxnode<=0 ) then
122            if ( rchsze>0 ) then
123              call QMDUPD(Xadj,Adjncy,rchsze,Rchset,Deg,Qsize,Qlink,Marker,Rchset(rchsze+1:),Nbrhd(nhdsze+1:))
124              Marker(node) = 0
125              do irch = 1, rchsze
126                inode = Rchset(irch)
127                if ( Marker(inode)>=0 ) then
128                  Marker(inode) = 0
129                  ndeg = Deg(inode)
130                  if ( ndeg<mindeg ) mindeg = ndeg
131                  if ( ndeg<=thresh ) then
132                    mindeg = thresh
133                    thresh = ndeg
134                    search = Invp(inode)
135                  endif
136                endif
137              enddo
138              if ( nhdsze>0 ) call QMDOT(node,Xadj,Adjncy,Marker,rchsze,Rchset,Nbrhd)
139            endif
140            if ( num>=Neqns ) exit
141            cycle loop2
142          endif
143        enddo
144        exit
145      enddo loop2
146      exit
147    enddo loop1
148
149    deallocate (DEG)
150    deallocate (MARker)
151    deallocate (RCHset)
152    deallocate (NBRhd)
153    deallocate (QSIze)
154    deallocate (QLInk)
155    deallocate (ADJncy)
156  end subroutine HECMW_ORDERING_GENQMD
157
158  !======================================================================!
159  !> @brief QMDRCH
160  !======================================================================!
161  subroutine QMDRCH(Root,Xadj,Adjncy,Deg,Marker,Rchsze,Rchset,Nhdsze,Nbrhd)
162    implicit none
163    !------
164    integer(kind=kint), intent(in):: Root
165    integer(kind=kint), intent(in):: Adjncy(:)
166    integer(kind=kint), intent(in):: Deg(:)
167    integer(kind=kint), intent(in):: Xadj(:)
168    integer(kind=kint), intent(out):: Nhdsze
169    integer(kind=kint), intent(out):: Rchsze
170    integer(kind=kint), intent(inout):: Marker(:)
171    integer(kind=kint), intent(out):: Rchset(:)
172    integer(kind=kint), intent(out):: Nbrhd(:)
173    !------
174    integer(kind=kint):: i
175    integer(kind=kint):: istrt
176    integer(kind=kint):: istop
177    integer(kind=kint):: j
178    integer(kind=kint):: jstrt
179    integer(kind=kint):: jstop
180    integer(kind=kint):: nabor
181    integer(kind=kint):: node
182
183    Nhdsze = 0
184    Rchsze = 0
185    istrt = Xadj(Root)
186    istop = Xadj(Root+1) - 1
187    if ( istop<istrt ) return
188    do i = istrt, istop
189      nabor = Adjncy(i)
190      if ( nabor==0 ) return
191      if ( Marker(nabor)==0 ) then
192        if ( Deg(nabor)<0 ) then
193          Marker(nabor) = -1
194          Nhdsze = Nhdsze + 1
195          Nbrhd(Nhdsze) = nabor
196          loop1: do
197            jstrt = Xadj(nabor)
198            jstop = Xadj(nabor+1) - 1
199            do j = jstrt, jstop
200              node = Adjncy(j)
201              nabor = -node
202              if ( node<0 ) cycle loop1
203              if ( node==0 ) exit
204              if ( Marker(node)==0 ) then
205                Rchsze = Rchsze + 1
206                Rchset(Rchsze) = node
207                Marker(node) = 1
208              endif
209            enddo
210            exit
211          enddo loop1
212        else
213          Rchsze = Rchsze + 1
214          Rchset(Rchsze) = nabor
215          Marker(nabor) = 1
216        endif
217      endif
218    enddo
219  end subroutine QMDRCH
220
221  !======================================================================!
222  !> @brief QMDUPD
223  !======================================================================!
224  subroutine QMDUPD(Xadj,Adjncy,Nlist,List,Deg,Qsize,Qlink,Marker,Rchset,Nbrhd)
225    implicit none
226    !------
227    integer(kind=kint), intent(in):: Nlist
228    integer(kind=kint), intent(in):: Adjncy(:)
229    integer(kind=kint), intent(in):: List(:)
230    integer(kind=kint), intent(in):: Xadj(:)
231    integer(kind=kint), intent(inout):: Deg(:)
232    integer(kind=kint), intent(inout):: Marker(:)
233    integer(kind=kint), intent(out):: Rchset(:)
234    integer(kind=kint), intent(out):: Nbrhd(:)
235    integer(kind=kint), intent(inout):: Qsize(:)
236    integer(kind=kint), intent(inout):: Qlink(:)
237    !------
238    integer(kind=kint):: deg0
239    integer(kind=kint):: deg1
240    integer(kind=kint):: il
241    integer(kind=kint):: inhd
242    integer(kind=kint):: inode
243    integer(kind=kint):: irch
244    integer(kind=kint):: j
245    integer(kind=kint):: jstrt
246    integer(kind=kint):: jstop
247    integer(kind=kint):: mark
248    integer(kind=kint):: nabor
249    integer(kind=kint):: nhdsze
250    integer(kind=kint):: node
251    integer(kind=kint):: rchsze
252
253    if ( Nlist<=0 ) return
254    deg0 = 0
255    nhdsze = 0
256    do il = 1, Nlist
257      node = List(il)
258      deg0 = deg0 + Qsize(node)
259      jstrt = Xadj(node)
260      jstop = Xadj(node+1) - 1
261      do j = jstrt, jstop
262        nabor = Adjncy(j)
263        if ( Marker(nabor)==0 .and. Deg(nabor)<0 ) then
264          Marker(nabor) = -1
265          nhdsze = nhdsze + 1
266          Nbrhd(nhdsze) = nabor
267        endif
268      enddo
269    enddo
270
271    if ( nhdsze>0 ) call QMDMRG(Xadj,Adjncy,Deg,Qsize,Qlink,Marker,deg0,nhdsze,Nbrhd,Rchset,Nbrhd(nhdsze+1:))
272    do il = 1, Nlist
273      node = List(il)
274      mark = Marker(node)
275      if ( mark<=1 .and. mark>=0 ) then
276        call QMDRCH(node,Xadj,Adjncy,Deg,Marker,rchsze,Rchset,nhdsze,Nbrhd)
277        deg1 = deg0
278        if ( rchsze>0 ) then
279          do irch = 1, rchsze
280            inode = Rchset(irch)
281            deg1 = deg1 + Qsize(inode)
282            Marker(inode) = 0
283          enddo
284        endif
285        Deg(node) = deg1 - 1
286        if ( nhdsze>0 ) then
287          do inhd = 1, nhdsze
288            inode = Nbrhd(inhd)
289            Marker(inode) = 0
290          enddo
291        endif
292      endif
293    enddo
294  end subroutine QMDUPD
295
296  !======================================================================!
297  !> @brief QMDMRG
298  !======================================================================!
299  subroutine QMDMRG(Xadj,Adjncy,Deg,Qsize,Qlink,Marker,Deg0,Nhdsze,Nbrhd,Rchset,Ovrlp)
300    implicit none
301    !------
302    integer(kind=kint), intent(in):: Deg0
303    integer(kind=kint), intent(in):: Nhdsze
304    integer(kind=kint), intent(in):: Adjncy(:)
305    integer(kind=kint), intent(in):: Nbrhd(:)
306    integer(kind=kint), intent(in):: Xadj(:)
307    integer(kind=kint), intent(inout):: Deg(:)
308    integer(kind=kint), intent(inout):: Qsize(:)
309    integer(kind=kint), intent(inout):: Qlink(:)
310    integer(kind=kint), intent(inout):: Marker(:)
311    integer(kind=kint), intent(out):: Rchset(:)
312    integer(kind=kint), intent(out):: Ovrlp(:)
313    !------
314    integer(kind=kint):: deg1
315    integer(kind=kint):: head
316    integer(kind=kint):: inhd
317    integer(kind=kint):: iov
318    integer(kind=kint):: irch
319    integer(kind=kint):: j
320    integer(kind=kint):: jstrt
321    integer(kind=kint):: jstop
322    integer(kind=kint):: link
323    integer(kind=kint):: lnode
324    integer(kind=kint):: mark
325    integer(kind=kint):: mrgsze
326    integer(kind=kint):: nabor
327    integer(kind=kint):: node
328    integer(kind=kint):: novrlp
329    integer(kind=kint):: rchsze
330    integer(kind=kint):: root
331
332    if ( Nhdsze<=0 ) return
333    do inhd = 1, Nhdsze
334      root = Nbrhd(inhd)
335      Marker(root) = 0
336    enddo
337    do inhd = 1, Nhdsze
338      root = Nbrhd(inhd)
339      Marker(root) = -1
340      rchsze = 0
341      novrlp = 0
342      deg1 = 0
343      loop1: do
344        jstrt = Xadj(root)
345        jstop = Xadj(root+1) - 1
346        do j = jstrt, jstop
347          nabor = Adjncy(j)
348          root = -nabor
349          if ( nabor<0 ) cycle loop1
350          if ( nabor==0 ) exit
351          mark = Marker(nabor)
352
353          if ( mark>=0 ) then
354            if ( mark<=0 ) then
355              rchsze = rchsze + 1
356              Rchset(rchsze) = nabor
357              deg1 = deg1 + Qsize(nabor)
358              Marker(nabor) = 1
359            elseif ( mark<=1 ) then
360              novrlp = novrlp + 1
361              Ovrlp(novrlp) = nabor
362              Marker(nabor) = 2
363            endif
364          endif
365        enddo
366        exit
367      enddo loop1
368      head = 0
369      mrgsze = 0
370      loop2: do iov = 1, novrlp
371        node = Ovrlp(iov)
372        jstrt = Xadj(node)
373        jstop = Xadj(node+1) - 1
374        do j = jstrt, jstop
375          nabor = Adjncy(j)
376          if ( Marker(nabor)==0 ) then
377            Marker(node) = 1
378            cycle loop2
379          endif
380        enddo
381        mrgsze = mrgsze + Qsize(node)
382        Marker(node) = -1
383        lnode = node
384        do
385          link = Qlink(lnode)
386          if ( link<=0 ) then
387            Qlink(lnode) = head
388            head = node
389            exit
390          else
391            lnode = link
392          endif
393        enddo
394      enddo loop2
395      if ( head>0 ) then
396        Qsize(head) = mrgsze
397        Deg(head) = Deg0 + deg1 - 1
398        Marker(head) = 2
399      endif
400      root = Nbrhd(inhd)
401      Marker(root) = 0
402      if ( rchsze>0 ) then
403        do irch = 1, rchsze
404          node = Rchset(irch)
405          Marker(node) = 0
406        enddo
407      endif
408    enddo
409  end subroutine QMDMRG
410
411  !======================================================================!
412  !> @brief QMDOT
413  !======================================================================!
414  subroutine QMDOT(Root,Xadj,Adjncy,Marker,Rchsze,Rchset,Nbrhd)
415    implicit none
416    !------
417    integer(kind=kint), intent(in):: Rchsze
418    integer(kind=kint), intent(in):: Root
419    integer(kind=kint), intent(in):: Marker(:)
420    integer(kind=kint), intent(in):: Rchset(:)
421    integer(kind=kint), intent(in):: Nbrhd(:)
422    integer(kind=kint), intent(in):: Xadj(:)
423    integer(kind=kint), intent(inout):: Adjncy(:)
424    !------
425    integer(kind=kint):: inhd
426    integer(kind=kint):: irch
427    integer(kind=kint):: j
428    integer(kind=kint):: jstrt
429    integer(kind=kint):: jstop
430    integer(kind=kint):: link
431    integer(kind=kint):: nabor
432    integer(kind=kint):: node
433
434    irch = 0
435    inhd = 0
436    node = Root
437    loop1: do
438      jstrt = Xadj(node)
439      jstop = Xadj(node+1) - 2
440      if ( jstop>=jstrt ) then
441        do j = jstrt, jstop
442          irch = irch + 1
443          Adjncy(j) = Rchset(irch)
444          if ( irch>=Rchsze ) exit loop1
445        enddo
446      endif
447      link = Adjncy(jstop+1)
448      node = -link
449      if ( link>=0 ) then
450        inhd = inhd + 1
451        node = Nbrhd(inhd)
452        Adjncy(jstop+1) = -node
453      endif
454    enddo loop1
455
456    Adjncy(j+1) = 0
457    do irch = 1, Rchsze
458      node = Rchset(irch)
459      if ( Marker(node)>=0 ) then
460        jstrt = Xadj(node)
461        jstop = Xadj(node+1) - 1
462        do j = jstrt, jstop
463          nabor = Adjncy(j)
464          if ( Marker(nabor)<0 ) then
465            Adjncy(j) = Root
466            exit
467          endif
468        enddo
469      endif
470    enddo
471  end subroutine QMDOT
472
473end module hecmw_ordering_qmd
474