1C
2C  This file is part of MUMPS 5.1.2, released
3C  on Mon Oct  2 07:37:01 UTC 2017
4C
5C
6C  Copyright 1991-2017 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
7C  University of Bordeaux.
8C
9C  This version of MUMPS is provided to you free of charge. It is
10C  released under the CeCILL-C license:
11C  http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
12C
13      MODULE DMUMPS_ANA_LR
14      USE DMUMPS_LR_CORE
15      USE DMUMPS_LR_STATS
16      USE MUMPS_LR_COMMON
17      USE MUMPS_ANA_ORD_WRAPPERS
18      IMPLICIT NONE
19      CONTAINS
20      SUBROUTINE GET_CUT(IWR, NASS, NCB, LRGROUPS, NPARTSCB,
21     &   NPARTSASS, CUT)
22      INTEGER, INTENT(IN) :: NASS, NCB
23      INTEGER, INTENT(IN) :: IWR(*)
24      INTEGER, INTENT(IN), DIMENSION(:) :: LRGROUPS
25      INTEGER, INTENT(OUT) :: NPARTSCB, NPARTSASS
26      INTEGER, POINTER, DIMENSION(:) :: CUT
27      INTEGER :: I, CURRENT_PART, CUTBUILDER
28      INTEGER, ALLOCATABLE, DIMENSION(:) :: BIG_CUT
29      ALLOCATE(BIG_CUT(max(NASS,1)+NCB+1))
30      CURRENT_PART = LRGROUPS(IWR(1))
31      BIG_CUT(1) = 1
32      BIG_CUT(2) = 2
33      CUTBUILDER = 2
34      NPARTSASS = 0
35      NPARTSCB  = 0
36      DO I = 2,NASS + NCB
37        IF (LRGROUPS(IWR(I)) == CURRENT_PART) THEN
38          BIG_CUT(CUTBUILDER) = BIG_CUT(CUTBUILDER) + 1
39        ELSE
40          CUTBUILDER = CUTBUILDER + 1
41          BIG_CUT(CUTBUILDER) = BIG_CUT(CUTBUILDER-1) + 1
42          CURRENT_PART = LRGROUPS(IWR(I))
43        END IF
44        IF (I == NASS) NPARTSASS = CUTBUILDER - 1
45      END DO
46      IF (NASS.EQ.1) NPARTSASS= 1
47      NPARTSCB = CUTBUILDER - 1 - NPARTSASS
48      ALLOCATE(CUT(max(NPARTSASS,1)+NPARTSCB+1))
49      IF (NPARTSASS.EQ.0) THEN
50        CUT(1) = 1
51        CUT(2:2+NPARTSCB) = BIG_CUT(1:1+NPARTSCB)
52      ELSE
53        CUT = BIG_CUT(1:NPARTSASS+NPARTSCB+1)
54      ENDIF
55      if(allocated(BIG_CUT)) DEALLOCATE(BIG_CUT)
56      END SUBROUTINE
57      SUBROUTINE SEP_GROUPING(NV, VLIST, N, NZ, LRGROUPS, NBGROUPS, IW,
58     &   LW, IPE, LEN, GROUP_SIZE, HALO_DEPTH, TRACE, WORKH, NODE,
59     &   GEN2HALO, K482, K472, K469, SEP_SIZE,
60     &   KEEP10, LP, LPOK, IFLAG, IERROR)
61        INTEGER(8), INTENT(IN) :: NZ, LW
62        INTEGER, INTENT(IN)    :: NV, N, GROUP_SIZE, HALO_DEPTH
63        INTEGER, INTENT(IN)    :: IW(LW), LEN(N), NODE, K482
64        INTEGER(8), INTENT(IN) :: IPE(N+1)
65        INTEGER, INTENT(IN)    :: K472, K469, SEP_SIZE, KEEP10, LP
66        LOGICAL                :: LPOK
67        INTEGER, INTENT(INOUT) :: NBGROUPS, WORKH(N)
68        INTEGER, INTENT(INOUT) :: LRGROUPS(N), VLIST(NV), TRACE(N)
69        INTEGER, INTENT(INOUT) :: GEN2HALO(N)
70        INTEGER, INTENT(INOUT) :: IFLAG, IERROR
71        INTEGER(8), ALLOCATABLE, DIMENSION(:) :: IPTRHALO
72        INTEGER,    ALLOCATABLE, DIMENSION(:) :: PARTS, JCNHALO
73        INTEGER(8) :: HALOEDGENBR
74        INTEGER    :: NHALO,
75     &       NBGROUPS_KWAY, I, GROUP_SIZE2, LRGROUPS_SIGN, IERR
76#if defined (metis) || defined (parmetis) || defined (metis4) || defined (parmetis3)
77        INTEGER :: METIS_IDX_SIZE
78#endif
79#if defined (scotch) || defined (ptscotch)
80        INTEGER :: SCOTCH_IDX_SIZE
81#endif
82        CALL COMPUTE_BLR_VCS(K472, GROUP_SIZE2, GROUP_SIZE, NV)
83        NBGROUPS_KWAY = MAX(NINT(dble(NV)/dble(GROUP_SIZE2)),1)
84        IF (NV .GE. SEP_SIZE) THEN
85          LRGROUPS_SIGN = 1
86        ELSE
87          LRGROUPS_SIGN = -1
88        ENDIF
89        IF (NBGROUPS_KWAY > 1) THEN
90          IF (K469.EQ.3) THEN
91!$OMP CRITICAL(gethalo_cri)
92            CALL GETHALONODES(N, IW, LW, IPE, VLIST, NV, HALO_DEPTH,
93     &                 NHALO, TRACE, WORKH, NODE, LEN, HALOEDGENBR,
94     &                 GEN2HALO)
95            ALLOCATE(PARTS(NHALO), IPTRHALO(NHALO+1),
96     &      JCNHALO(HALOEDGENBR), STAT=IERR)
97            IF (IERR.GT.0) THEN
98              IF (LPOK) WRITE(LP,*)
99     &        " Error allocate integer array of size: ",
100     &         int(NHALO+(KEEP10*(NHALO+1)),8) + HALOEDGENBR
101              IFLAG  = -7
102              CALL MUMPS_SET_IERROR
103     &         (int(NHALO+(KEEP10*(NHALO+1)),8) + HALOEDGENBR,
104     &          IERROR)
105            ENDIF
106            CALL GETHALOGRAPH(WORKH, NHALO, N, IW, LW, IPE, IPTRHALO,
107     &     JCNHALO, HALOEDGENBR,TRACE,NODE, GEN2HALO)
108!$OMP END CRITICAL(gethalo_cri)
109            IF (IFLAG.LT.0) RETURN
110          ELSE
111            CALL GETHALONODES(N, IW, LW, IPE, VLIST, NV, HALO_DEPTH,
112     &                 NHALO, TRACE, WORKH, NODE, LEN, HALOEDGENBR,
113     &                 GEN2HALO)
114            ALLOCATE(PARTS(NHALO), IPTRHALO(NHALO+1),
115     &      JCNHALO(HALOEDGENBR), STAT=IERR)
116            IF (IERR.GT.0) THEN
117              IF (LPOK) WRITE(LP,*)
118     &        " Error allocate integer array of size: ",
119     &         int(NHALO+(KEEP10*(NHALO+1)),8) + HALOEDGENBR
120              IFLAG  = -7
121              CALL MUMPS_SET_IERROR
122     &         (int(NHALO+(KEEP10*(NHALO+1)),8) + HALOEDGENBR,
123     &          IERROR)
124              RETURN
125            ENDIF
126            CALL GETHALOGRAPH(WORKH, NHALO, N, IW, LW, IPE, IPTRHALO,
127     &     JCNHALO, HALOEDGENBR,TRACE,NODE, GEN2HALO)
128          ENDIF
129          IF (K482.EQ.1) THEN
130#if defined (metis) || defined (parmetis) || defined (metis4) || defined (parmetis3)
131             CALL MUMPS_METIS_IDXSIZE(METIS_IDX_SIZE)
132             IF (METIS_IDX_SIZE .EQ. 64) THEN
133                CALL MUMPS_METIS_KWAY_MIXEDto64(NHALO, HALOEDGENBR,
134     &               IPTRHALO,
135     &               JCNHALO,
136     &               NBGROUPS_KWAY, PARTS, LP, LPOK, KEEP10,
137     &               IFLAG, IERROR)
138             ELSE
139               IF (KEEP10.EQ.1) THEN
140                IFLAG  = -52
141                IERROR = 1
142               ELSE
143                CALL MUMPS_METIS_KWAY_MIXEDto32(NHALO, HALOEDGENBR,
144     &               IPTRHALO,
145     &               JCNHALO,
146     &               NBGROUPS_KWAY, PARTS, LP, LPOK, KEEP10,
147     &               IFLAG, IERROR)
148               ENDIF
149             ENDIF
150#endif
151          ELSE IF (K482.EQ.2) THEN
152#if defined (scotch) || defined (ptscotch)
153             CALL MUMPS_SCOTCH_INTSIZE(SCOTCH_IDX_SIZE)
154             IF (SCOTCH_IDX_SIZE .EQ. 32) THEN
155               IF (KEEP10.EQ.1) THEN
156                IFLAG  = -52
157                IERROR = 2
158               ELSE
159                CALL MUMPS_SCOTCH_KWAY_MIXEDto32(
160     &               NHALO, HALOEDGENBR, IPTRHALO, JCNHALO,
161     &               NBGROUPS_KWAY, PARTS, LP, LPOK, KEEP10,
162     &               IFLAG, IERROR)
163               ENDIF
164             ELSE
165                CALL MUMPS_SCOTCH_KWAY_MIXEDto64(
166     &               NHALO, HALOEDGENBR, IPTRHALO, JCNHALO,
167     &               NBGROUPS_KWAY, PARTS, LP, LPOK, KEEP10,
168     &               IFLAG, IERROR)
169             END IF
170#endif
171          ELSE
172            WRITE(6,*) " Internal ERROR K482=", K482
173            CALL MUMPS_ABORT()
174          END IF
175          IF (IFLAG.LT.0) GOTO 500
176          CALL GET_GLOBAL_GROUPS(PARTS,VLIST,NV,
177     &     NBGROUPS_KWAY, LRGROUPS, N, NBGROUPS, LRGROUPS_SIGN)
178        ELSE
179!$OMP CRITICAL(lrgrouping_cri)
180          DO I=1,NV
181            LRGROUPS(VLIST(I)) = LRGROUPS_SIGN*(NBGROUPS + 1)
182          END DO
183            NBGROUPS = NBGROUPS + 1
184!$OMP END CRITICAL(lrgrouping_cri)
185        END IF
186  500   IF (allocated(IPTRHALO)) then
187           DEALLOCATE(IPTRHALO)
188        ENDIF
189        IF (allocated(PARTS)) then
190           DEALLOCATE(PARTS)
191        ENDIF
192        IF (allocated(JCNHALO)) then
193           DEALLOCATE(JCNHALO   )
194        ENDIF
195        RETURN
196      END SUBROUTINE SEP_GROUPING
197      SUBROUTINE GET_GLOBAL_GROUPS(PARTS, SEP, NSEP, NPARTS,
198     &                    LRGROUPS, N, NBGROUPS, LRGROUPS_SIGN)
199      INTEGER,INTENT(IN) :: NSEP, N, LRGROUPS_SIGN
200      INTEGER :: PARTS(:)
201      INTEGER,DIMENSION(:),INTENT(INOUT) :: SEP
202      INTEGER, INTENT(INOUT) :: NPARTS
203      INTEGER, INTENT(INOUT) :: NBGROUPS
204      INTEGER, INTENT(INOUT) :: LRGROUPS(N)
205      INTEGER::I,CNT,NB_PARTS_WITHOUT_SEP_NODE
206      INTEGER,DIMENSION(:),ALLOCATABLE::SIZES, RIGHTPART
207      INTEGER,DIMENSION(:),ALLOCATABLE::PARTPTR
208      INTEGER,DIMENSION(:),ALLOCATABLE :: NEWSEP
209      ALLOCATE( NEWSEP(NSEP),
210     &          SIZES(NPARTS),
211     &          RIGHTPART(NPARTS),
212     &          PARTPTR(NPARTS+1))
213      NB_PARTS_WITHOUT_SEP_NODE = 0
214      RIGHTPART = 0
215      SIZES = 0
216      DO I=1,NSEP
217        SIZES(PARTS(I)) = SIZES(PARTS(I)) + 1
218      END DO
219      PARTPTR(1)=1
220      CNT = 0
221      DO I=2,NPARTS+1
222        PARTPTR(I) = PARTPTR(I-1) + SIZES(I-1)
223        IF (SIZES(I-1)==0) THEN
224          NB_PARTS_WITHOUT_SEP_NODE = NB_PARTS_WITHOUT_SEP_NODE + 1
225        ELSE
226          CNT = CNT + 1
227          RIGHTPART(I-1) = CNT
228        END IF
229      END DO
230      NPARTS = NPARTS - NB_PARTS_WITHOUT_SEP_NODE
231!$OMP CRITICAL(lrgrouping_cri)
232      DO I=1,NSEP
233        NEWSEP(PARTPTR(PARTS(I))) = SEP(I)
234        LRGROUPS(SEP(I)) = LRGROUPS_SIGN*(RIGHTPART(PARTS(I))
235     &    + NBGROUPS)
236        PARTPTR(PARTS(I)) =
237     &    PARTPTR(PARTS(I)) + 1
238      END DO
239      NBGROUPS = NBGROUPS + NPARTS
240!$OMP END CRITICAL(lrgrouping_cri)
241      SEP = NEWSEP
242      DEALLOCATE(NEWSEP,SIZES,RIGHTPART,PARTPTR)
243      END SUBROUTINE GET_GLOBAL_GROUPS
244      SUBROUTINE GETHALONODES(N, IW, LW, IPE, IND, NIND, PMAX,
245     &                        NHALO, TRACE, WORKH, NODE, LEN, CNT,
246     &                        GEN2HALO)
247        INTEGER,DIMENSION(:),INTENT(IN) :: IND
248        INTEGER(8), INTENT(IN)          :: LW
249        INTEGER, INTENT(IN)             :: N, NODE
250        INTEGER, INTENT(IN)             :: IW(LW), LEN(N)
251        INTEGER(8), INTENT(IN)          :: IPE(N+1)
252        INTEGER, INTENT(IN)             :: PMAX,NIND
253        INTEGER, INTENT(OUT)            :: NHALO
254        INTEGER, INTENT(INOUT)          :: TRACE(N), WORKH(N)
255        INTEGER                         :: GEN2HALO(N)
256        INTEGER(8), INTENT(OUT)         :: CNT
257        INTEGER                         :: DEPTH, I, LAST_LVL_START
258        INTEGER                         :: HALOI
259        INTEGER(8)                      :: J
260        WORKH(1:NIND) = IND
261        LAST_LVL_START = 1
262        NHALO = NIND
263        CNT = 0
264        DO I=1,NIND
265          HALOI = WORKH(I)
266          GEN2HALO(HALOI) = I
267          IF (TRACE(HALOI) .NE. NODE) THEN
268             TRACE(HALOI) = NODE
269          END IF
270          DO J=IPE(HALOI),IPE(HALOI+1)-1
271            IF (TRACE(IW(J)).EQ.NODE) THEN
272              CNT = CNT + 2
273            END IF
274          END DO
275        END DO
276        DO DEPTH=1,PMAX
277          CALL NEIGHBORHOOD(WORKH, NHALO, N, IW, LW, IPE,
278     &                      TRACE, NODE, LEN, CNT, LAST_LVL_START,
279     &                      DEPTH, PMAX, GEN2HALO)
280        END DO
281      END SUBROUTINE
282      SUBROUTINE NEIGHBORHOOD(HALO, NHALO, N, IW, LW, IPE,
283     &                        TRACE, NODE, LEN, CNT, LAST_LVL_START,
284     &                        DEPTH, PMAX, GEN2HALO)
285        INTEGER, INTENT(IN)                 :: N, NODE, DEPTH, PMAX
286        INTEGER,INTENT(INOUT)               :: NHALO, GEN2HALO(N)
287        INTEGER, INTENT(INOUT)              :: LAST_LVL_START
288        INTEGER(8), INTENT(INOUT)           :: CNT
289        INTEGER,DIMENSION(:),INTENT(INOUT)  :: HALO
290        INTEGER(8), INTENT(IN)              :: LW
291        INTEGER(8), INTENT(IN)              :: IPE(N+1)
292        INTEGER, TARGET, INTENT(IN)         :: IW(LW)
293        INTEGER, INTENT(IN)                 :: LEN(N)
294        INTEGER,DIMENSION(:)                :: TRACE
295        INTEGER :: AvgDens, THRESH
296        INTEGER :: I,INEI,NADJI,NEWNHALO, NEIGH
297        INTEGER, DIMENSION(:), POINTER :: ADJI
298        INTEGER(8) :: J
299        NEWNHALO = 0
300        AvgDens = nint(dble(IPE(N+1)-1_8)/dble(N))
301        THRESH  = AvgDens*10
302        DO I=LAST_LVL_START,NHALO
303          NADJI = LEN(HALO(I))
304          IF (NADJI.GT.THRESH) CYCLE
305          ADJI => IW(IPE(HALO(I)):IPE(HALO(I)+1)-1)
306          DO INEI=1,NADJI
307            IF (TRACE(ADJI(INEI)) .NE. NODE) THEN
308              NEIGH = ADJI(INEI)
309              IF (LEN(NEIGH).GT.THRESH) CYCLE
310              TRACE(NEIGH) = NODE
311              NEWNHALO = NEWNHALO + 1
312              HALO(NHALO+NEWNHALO) = NEIGH
313              GEN2HALO(NEIGH) = NHALO + NEWNHALO
314              DO J=IPE(NEIGH),IPE(NEIGH+1)-1
315                IF (TRACE(IW(J)).EQ.NODE) THEN
316                  CNT = CNT + 2
317                END IF
318              END DO
319            END IF
320          END DO
321        END DO
322        LAST_LVL_START = NHALO + 1
323        NHALO = NHALO + NEWNHALO
324      END SUBROUTINE
325      SUBROUTINE GETHALOGRAPH(HALO,NHALO,N,IW,LW,IPE,IPTRHALO,JCNHALO,
326     &     HALOEDGENBR,TRACE,NODE, GEN2HALO)
327      INTEGER, INTENT(IN) :: N
328      INTEGER,INTENT(IN):: NHALO, NODE
329      INTEGER,INTENT(IN):: GEN2HALO(N)
330      INTEGER,DIMENSION(NHALO),INTENT(IN) :: HALO
331      INTEGER(8), INTENT(IN)              :: LW
332      INTEGER(8), INTENT(IN)              :: IPE(N+1)
333      INTEGER, INTENT(IN)     :: IW(LW), TRACE(N)
334      INTEGER(8),INTENT(IN)   :: HALOEDGENBR
335      INTEGER(8), INTENT(OUT) :: IPTRHALO(NHALO+1)
336      INTEGER, INTENT(OUT)    :: JCNHALO(HALOEDGENBR)
337      INTEGER::I,IPTR_CNT,JCN_CNT,HALOI
338      INTEGER(8) :: J, CNT
339      CNT = 0
340      IPTR_CNT = 2
341      JCN_CNT = 1
342      IPTRHALO(1) = 1
343      DO I=1,NHALO
344         HALOI = HALO(I)
345         DO J=IPE(HALOI),IPE(HALOI+1)-1
346            IF (TRACE(IW(J))==NODE) THEN
347               CNT = CNT + 1
348               JCNHALO(JCN_CNT) = GEN2HALO(IW(J))
349               JCN_CNT = JCN_CNT + 1
350            END IF
351         END DO
352         IPTRHALO(IPTR_CNT) = CNT + 1
353         IPTR_CNT = IPTR_CNT + 1
354      END DO
355      END SUBROUTINE
356      SUBROUTINE GET_GROUPS(NHALO,PARTS,SEP,NSEP,NPARTS,
357     &     CUT,NEWSEP,PERM,IPERM)
358      INTEGER,INTENT(IN) :: NHALO,NSEP
359      INTEGER,DIMENSION(:),INTENT(IN) :: SEP
360      INTEGER,POINTER,DIMENSION(:)::PARTS
361      INTEGER,POINTER,DIMENSION(:)::CUT,NEWSEP,PERM,
362     &   IPERM
363      INTEGER,INTENT(INOUT) :: NPARTS
364      INTEGER::I,CNT,NB_PARTS_WITHOUT_SEP_NODE
365      INTEGER,DIMENSION(:),ALLOCATABLE::SIZES
366      INTEGER,DIMENSION(:),ALLOCATABLE::PARTPTR
367      ALLOCATE(NEWSEP(NSEP))
368      ALLOCATE(PERM(NSEP))
369      ALLOCATE(IPERM(NSEP))
370      ALLOCATE(SIZES(NPARTS))
371      ALLOCATE(PARTPTR(NPARTS+1))
372      NB_PARTS_WITHOUT_SEP_NODE = 0
373      SIZES = 0
374      DO I=1,NSEP
375        SIZES(PARTS(I)) =
376     &     SIZES(PARTS(I))+1
377      END DO
378      PARTPTR(1)=1
379      DO I=2,NPARTS+1
380        PARTPTR(I) = PARTPTR(I-1) + SIZES(I-1)
381        IF (SIZES(I-1)==0) THEN
382          NB_PARTS_WITHOUT_SEP_NODE = NB_PARTS_WITHOUT_SEP_NODE + 1
383        END IF
384      END DO
385      ALLOCATE(CUT(NPARTS-NB_PARTS_WITHOUT_SEP_NODE+1))
386      CUT(1) = 1
387      CNT = 2
388      DO I=2,NPARTS+1
389        IF (SIZES(I-1).NE.0) THEN
390          CUT(CNT) = PARTPTR(I)
391          CNT = CNT + 1
392        END IF
393      END DO
394      NPARTS = NPARTS - NB_PARTS_WITHOUT_SEP_NODE
395      CUT(NPARTS+1) = NSEP+1
396      DO I=1,NSEP
397        NEWSEP(PARTPTR(PARTS(I))) = SEP(I)
398        PERM(PARTPTR(PARTS(I))) = I
399        IPERM(I) = PARTPTR(PARTS(I))
400        PARTPTR(PARTS(I)) =
401     &     PARTPTR(PARTS(I)) + 1
402      END DO
403      DEALLOCATE(SIZES,PARTPTR)
404      END SUBROUTINE
405      END MODULE DMUMPS_ANA_LR
406