1C
2C  This file is part of MUMPS 4.10.0, built on Tue May 10 12:56:32 UTC 2011
3C
4C
5C  This version of MUMPS is provided to you free of charge. It is public
6C  domain, based on public domain software developed during the Esprit IV
7C  European project PARASOL (1996-1999). Since this first public domain
8C  version in 1999, research and developments have been supported by the
9C  following institutions: CERFACS, CNRS, ENS Lyon, INPT(ENSEEIHT)-IRIT,
10C  INRIA, and University of Bordeaux.
11C
12C  The MUMPS team at the moment of releasing this version includes
13C  Patrick Amestoy, Maurice Bremond, Alfredo Buttari, Abdou Guermouche,
14C  Guillaume Joslin, Jean-Yves L'Excellent, Francois-Henry Rouet, Bora
15C  Ucar and Clement Weisbecker.
16C
17C  We are also grateful to Emmanuel Agullo, Caroline Bousquet, Indranil
18C  Chowdhury, Philippe Combes, Christophe Daniel, Iain Duff, Vincent Espirat,
19C  Aurelia Fevre, Jacko Koster, Stephane Pralet, Chiara Puglisi, Gregoire
20C  Richard, Tzvetomila Slavova, Miroslav Tuma and Christophe Voemel who
21C  have been contributing to this project.
22C
23C  Up-to-date copies of the MUMPS package can be obtained
24C  from the Web pages:
25C  http://mumps.enseeiht.fr/  or  http://graal.ens-lyon.fr/MUMPS
26C
27C
28C   THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
29C   EXPRESSED OR IMPLIED. ANY USE IS AT YOUR OWN RISK.
30C
31C
32C  User documentation of any code that uses this software can
33C  include this complete notice. You can acknowledge (using
34C  references [1] and [2]) the contribution of this package
35C  in any scientific publication dependent upon the use of the
36C  package. You shall use reasonable endeavours to notify
37C  the authors of the package of this publication.
38C
39C   [1] P. R. Amestoy, I. S. Duff, J. Koster and  J.-Y. L'Excellent,
40C   A fully asynchronous multifrontal solver using distributed dynamic
41C   scheduling, SIAM Journal of Matrix Analysis and Applications,
42C   Vol 23, No 1, pp 15-41 (2001).
43C
44C   [2] P. R. Amestoy and A. Guermouche and J.-Y. L'Excellent and
45C   S. Pralet, Hybrid scheduling for the parallel solution of linear
46C   systems. Parallel Computing Vol 32 (2), pp 136-156 (2006).
47C
48       MODULE MUMPS_SOL_ES
49      PRIVATE
50      PUBLIC:: PRUNED_SIZE_LOADED
51      PUBLIC:: MUMPS_797
52      PUBLIC:: MUMPS_802
53      PUBLIC:: MUMPS_798
54      PUBLIC:: MUMPS_803
55      PUBLIC:: MUMPS_804
56      INTEGER(8), POINTER, DIMENSION(:,:) :: SIZE_OF_BLOCK
57      INTEGER(8) :: PRUNED_SIZE_LOADED
58      CONTAINS
59      SUBROUTINE MUMPS_804(SIZE_OF_BLOCK_ARG, KEEP201)
60      IMPLICIT NONE
61      INTEGER, INTENT(IN) :: KEEP201
62      INTEGER(8), POINTER, DIMENSION(:,:) :: SIZE_OF_BLOCK_ARG
63      IF (KEEP201 > 0) THEN
64        SIZE_OF_BLOCK => SIZE_OF_BLOCK_ARG
65      ELSE
66        NULLIFY(SIZE_OF_BLOCK)
67      ENDIF
68      RETURN
69      END SUBROUTINE MUMPS_804
70      SUBROUTINE MUMPS_798(
71     &     fill,
72     &     DAD, NE_STEPS, FRERE, KEEP28,
73     &     FILS, STEP, N,
74     &     nodes_RHS, nb_nodes_RHS,
75     &     TO_PROCESS,
76     &     nb_prun_nodes, nb_prun_roots, nb_prun_leaves,
77     &     Pruned_List, Pruned_Roots, Pruned_Leaves
78     &     )
79      IMPLICIT NONE
80      LOGICAL, INTENT(IN) :: fill
81      INTEGER, INTENT(IN) :: N, KEEP28
82      INTEGER, INTENT(IN) :: DAD(KEEP28),NE_STEPS(KEEP28),FRERE(KEEP28)
83      INTEGER, INTENT(IN) :: FILS(N), STEP(N)
84      INTEGER, INTENT(IN) :: nodes_RHS(KEEP28),  nb_nodes_RHS
85      INTEGER :: nb_prun_nodes
86      INTEGER, OPTIONAL, INTENT(INOUT):: Pruned_List(nb_prun_nodes)
87      INTEGER :: nb_prun_roots
88      INTEGER, OPTIONAL, INTENT(INOUT):: Pruned_Roots(nb_prun_roots)
89      INTEGER :: nb_prun_leaves
90      INTEGER, OPTIONAL, INTENT(INOUT):: Pruned_Leaves(nb_prun_leaves)
91      LOGICAL :: TO_PROCESS(KEEP28)
92      INTEGER :: IN, I, ISTEP, TMP, TMPsave
93      nb_prun_nodes = 0
94      nb_prun_leaves = 0
95      TO_PROCESS(:) = .FALSE.
96      DO I = 1, nb_nodes_RHS
97         TMP = nodes_RHS(I)
98         TMPsave = TMP
99         ISTEP = STEP(TMP)
100         DO WHILE(.NOT.TO_PROCESS(ISTEP))
101            TO_PROCESS(ISTEP) = .TRUE.
102            nb_prun_nodes = nb_prun_nodes + 1
103            IF(fill) THEN
104               Pruned_List(nb_prun_nodes) = TMP
105            END IF
106            IN = FILS(TMP)
107            DO WHILE(IN.GT.0)
108               IN = FILS(IN)
109            END DO
110            IF (IN.LT.0) THEN
111               TMP = -IN
112               ISTEP = STEP(TMP)
113            ELSE
114               nb_prun_leaves = nb_prun_leaves + 1
115               IF(fill) THEN
116                  Pruned_Leaves(nb_prun_leaves) = TMP
117               END IF
118               IF(TMP.NE.TMPsave) THEN
119                  TMP = abs(FRERE(ISTEP))
120                  IF(TMP.NE.0) THEN
121                     ISTEP = STEP(TMP)
122                  END IF
123               END IF
124            END IF
125         END DO
126      END DO
127      nb_prun_roots = 0
128      DO I=1,nb_nodes_RHS
129         TMP = nodes_RHS(I)
130         ISTEP = STEP(TMP)
131         IF(DAD(ISTEP).NE.0) THEN
132            IF(.NOT.TO_PROCESS(STEP(DAD(ISTEP)))) THEN
133               nb_prun_roots = nb_prun_roots + 1
134               IF(fill) THEN
135                  Pruned_Roots(nb_prun_roots) = TMP
136               END IF
137            END IF
138         ELSE
139            nb_prun_roots = nb_prun_roots + 1
140            IF(fill) THEN
141               Pruned_Roots(nb_prun_roots) = TMP
142            END IF
143         END IF
144      END DO
145      RETURN
146      END SUBROUTINE MUMPS_798
147      SUBROUTINE MUMPS_797(
148     &     fill,
149     &     DAD, KEEP28,
150     &     STEP, N,
151     &     nodes_RHS, nb_nodes_RHS,
152     &     Pruned_SONS, TO_PROCESS,
153     &     nb_prun_nodes,nb_prun_roots, nb_prun_leaves,
154     &     Pruned_List, Pruned_Roots, Pruned_Leaves
155     &     )
156      IMPLICIT NONE
157      LOGICAL, INTENT(IN) :: fill
158      INTEGER, INTENT(IN) :: N
159      INTEGER, INTENT(IN) :: STEP(N)
160      INTEGER, INTENT(IN) :: KEEP28
161      INTEGER, INTENT(IN) :: DAD(KEEP28)
162      INTEGER, INTENT(IN) :: nb_nodes_RHS
163      INTEGER, INTENT(IN) :: nodes_RHS(nb_nodes_RHS)
164      INTEGER :: nb_prun_nodes
165      INTEGER, OPTIONAL, INTENT(INOUT):: Pruned_List(nb_prun_nodes)
166      INTEGER :: nb_prun_roots
167      INTEGER, OPTIONAL, INTENT(INOUT):: Pruned_Roots(nb_prun_roots)
168      INTEGER :: nb_prun_leaves
169      INTEGER, OPTIONAL, INTENT(INOUT):: Pruned_Leaves(nb_prun_leaves)
170      INTEGER :: Pruned_SONS(KEEP28)
171      LOGICAL :: TO_PROCESS(KEEP28)
172      INTEGER :: IN, I, ISTEP, TMP
173      nb_prun_nodes = 0
174      nb_prun_roots = 0
175      TO_PROCESS(:) = .FALSE.
176      Pruned_SONS(:) = -1
177      DO I = 1, nb_nodes_RHS
178         TMP = nodes_RHS(I)
179         ISTEP = STEP(TMP)
180         TO_PROCESS(ISTEP) = .TRUE.
181         IF (Pruned_SONS(ISTEP) .eq. -1) THEN
182            Pruned_SONS(ISTEP) = 0
183            nb_prun_nodes = nb_prun_nodes + 1
184            IF(fill) THEN
185               Pruned_List(nb_prun_nodes) = nodes_RHS(I)
186            END IF
187            IN = nodes_RHS(I)
188            IN = DAD(STEP(IN))
189            DO WHILE (IN.NE.0)
190               TO_PROCESS(STEP(IN)) = .TRUE.
191               IF (Pruned_SONS(STEP(IN)).eq.-1) THEN
192                  nb_prun_nodes = nb_prun_nodes + 1
193                  IF(fill) THEN
194                     Pruned_List(nb_prun_nodes) = IN
195                  END IF
196                  Pruned_SONS(STEP(IN)) = 1
197                  TMP = IN
198                  IN = DAD(STEP(IN))
199               ELSE
200                  Pruned_SONS(STEP(IN)) = Pruned_SONS(STEP(IN)) + 1
201                  GOTO 201
202               ENDIF
203            ENDDO
204            nb_prun_roots = nb_prun_roots +1
205            IF(fill) THEN
206               Pruned_Roots(nb_prun_roots) = TMP
207            END IF
208         ENDIF
209  201    CONTINUE
210      ENDDO
211      nb_prun_leaves = 0
212      DO I = 1, nb_nodes_RHS
213         TMP = nodes_RHS(I)
214         ISTEP = STEP(TMP)
215         IF (Pruned_SONS(ISTEP).EQ.0) THEN
216            nb_prun_leaves = nb_prun_leaves +1
217            IF(fill) THEN
218              Pruned_Leaves(nb_prun_leaves) = TMP
219            END IF
220         END IF
221      ENDDO
222      RETURN
223      END SUBROUTINE MUMPS_797
224      SUBROUTINE MUMPS_803(MYID, N, KEEP28, KEEP201,
225     &           KEEP8_31,
226     &           STEP, Pruned_List, nb_prun_nodes, OOC_FCT_TYPE_LOC)
227      INTEGER, intent(in) :: KEEP28, KEEP201, OOC_FCT_TYPE_LOC, MYID, N
228      INTEGER(8), intent(in) :: KEEP8_31
229      INTEGER, intent(in) :: nb_prun_nodes
230      INTEGER, intent(in) :: Pruned_List(nb_prun_nodes)
231      INTEGER, intent(in) :: STEP(N)
232      INTEGER I, ISTEP
233      INTEGER(8) :: Pruned_Size
234#if defined(Mila_Print)
235      write(*,*) ' in Pruned List nodes:',nb_prun_nodes
236      write(*,*) Pruned_List(1:nb_prun_nodes)
237#endif
238      IF (KEEP201 .GT. 0) THEN
239        Pruned_Size = 0_8
240        DO I = 1, nb_prun_nodes
241          ISTEP = STEP(Pruned_List(I))
242          Pruned_Size = Pruned_Size + SIZE_OF_BLOCK
243     &                  (ISTEP, OOC_FCT_TYPE_LOC)
244        ENDDO
245        PRUNED_SIZE_LOADED = PRUNED_SIZE_LOADED +Pruned_Size
246#if defined(Mila_Print)
247        write(*,*) 'Pruned_Size Total_Size:',
248     &        Pruned_Size, KEEP8_31
249        write(*,*) MYID,'Gain (%) = ', dble(100)
250     &        - (dble(Pruned_Size)*dble(100)) /dble(KEEP8_31)
251        IF (Pruned_Size .EQ. 0) THEN
252            WRITE(*,*) "NOT NORMAL BEHAVIOUR !!"
253            DO I = 1, nb_nodes_RHS
254               WRITE(*,*) "starting_node node_size",
255     &              nodes_RHS(I),
256     &              SIZE_OF_BLOCK(STEP(nodes_RHS(I)),OOC_FCT_TYPE_LOC)
257            ENDDO
258        ENDIF
259        write(*,*) '============================='
260#endif
261      ENDIF
262      RETURN
263      END SUBROUTINE MUMPS_803
264      SUBROUTINE MUMPS_802
265     &                (MYID, N, KEEP28, KEEP201, KEEP8_31,
266     &                STEP, Pruned_List, nb_prun_nodes, OOC_FCT_TYPE_LOC
267     & )
268      IMPLICIT NONE
269      INTEGER, intent(in) :: KEEP28, KEEP201, OOC_FCT_TYPE_LOC, N
270      INTEGER(8), intent(in) :: KEEP8_31
271      INTEGER, intent(in) :: nb_prun_nodes, MYID
272      INTEGER, intent(in) :: Pruned_List(nb_prun_nodes)
273      INTEGER, intent(in) :: STEP(N)
274      INCLUDE 'mpif.h'
275      INTEGER I, ISTEP
276      INTEGER(8) :: Pruned_Size
277      Pruned_Size = 0_8
278      DO I = 1, nb_prun_nodes
279        ISTEP = STEP(Pruned_List(I))
280        IF (KEEP201 .GT. 0) THEN
281            Pruned_Size = Pruned_Size + SIZE_OF_BLOCK
282     &                    (ISTEP, OOC_FCT_TYPE_LOC)
283        ENDIF
284      ENDDO
285      IF (KEEP201.GT.0) THEN
286#       if defined(Mila_Print)
287          write(*,*) MYID,'PR leaves NODES',nb_prun_leaves,
288     &     Pruned_Leaves(1:nb_prun_leaves)
289          write(*,*) MYID,'PR  NODES',Pos_List,
290     &     Pruned_List(1:Pos_List)
291          write(*,*) 'PR root NODES',
292     &     Pruned_Roots(nb_prun_roots)
293#       endif
294        IF (KEEP8_31 .NE. 0_8) THEN
295#         if defined(Mila_Print)
296            write(*,*) MYID,'PRUNED and TOTAL Size:',
297     &        Pruned_Size, KEEP8_31
298            write(*,*) MYID,'Gain (%) = ', dble(100)
299     &        - ((dble(Pruned_Size)*dble(100))/dble(KEEP8_31))
300            IF (MYID.EQ.0)
301     &        write(*,*) '============================='
302#         endif
303          PRUNED_SIZE_LOADED = PRUNED_SIZE_LOADED +Pruned_Size
304        ENDIF
305      ENDIF
306      RETURN
307      END SUBROUTINE MUMPS_802
308      END MODULE MUMPS_SOL_ES
309      SUBROUTINE MUMPS_780
310     &          (PERM_STRAT, SYM_PERM,
311     &           IRHS_PTR, NHRS,
312     &           PERM_RHS, SIZEPERM, IERR
313     &         )
314      IMPLICIT NONE
315      INTEGER, INTENT(IN) :: PERM_STRAT, NHRS, SIZEPERM
316      INTEGER, INTENT(IN) :: SYM_PERM(SIZEPERM)
317      INTEGER, INTENT(IN) :: IRHS_PTR(NHRS)
318      INTEGER, INTENT(OUT):: IERR
319      INTEGER, INTENT(OUT):: PERM_RHS(SIZEPERM)
320      DOUBLE PRECISION :: RAND_NUM
321      INTEGER  I, J, STRAT
322      IERR = 0
323      STRAT = PERM_STRAT
324      IF( (STRAT.NE.-3).AND.
325     &    (STRAT.NE.-2).AND.
326     &    (STRAT.NE.-1).AND.
327     &    (STRAT.NE. 1).AND.
328     &    (STRAT.NE. 2).AND.
329     &    (STRAT.NE. 6) ) THEN
330        WRITE(*,*)"Warning: incorrect value for the RHS permutation; ",
331     &            "defaulting to post-order"
332        STRAT = 1
333      END IF
334      IF (STRAT .EQ. -3) THEN
335         WRITE(*,*) "Processing the RHS in random order"
336         PERM_RHS(1:SIZEPERM)=0
337         DO I=1, SIZEPERM
338           CALL random_number(RAND_NUM)
339           RAND_NUM = RAND_NUM*dble(SIZEPERM)
340           J = ceiling(RAND_NUM)
341           DO WHILE (PERM_RHS(J).NE.0)
342             CALL random_number(RAND_NUM)
343             RAND_NUM = RAND_NUM*dble(SIZEPERM)
344             J = ceiling(RAND_NUM)
345           ENDDO
346           PERM_RHS(J)=I
347         ENDDO
348      ELSEIF (STRAT .EQ. -2) THEN
349         WRITE(*,*) "Processing the RHS in inverse order"
350         DO I=1, SIZEPERM
351            PERM_RHS(SIZEPERM -I +1) = I
352         ENDDO
353      ELSEIF (STRAT .EQ. -1) THEN
354         WRITE(*,*) "Processing the RHS in natural order"
355         DO I=1, SIZEPERM
356            PERM_RHS(I) = I
357         ENDDO
358      ELSEIF (STRAT .EQ.  1) THEN
359         WRITE(*,*) "Processing the RHS in post-order"
360         DO I=1, SIZEPERM
361            PERM_RHS(SYM_PERM(I)) = I
362         ENDDO
363      ELSEIF (STRAT .EQ.  2) THEN
364         WRITE(*,*) "Processing the RHS in pre-order"
365         DO I=1, SIZEPERM
366            PERM_RHS(SIZEPERM-SYM_PERM(I)+1) = I
367         ENDDO
368      ENDIF
369      END SUBROUTINE MUMPS_780
370      SUBROUTINE MUMPS_772
371     &    (PERM_RHS, SIZEPERM, N, KEEP_28,
372     &     PROCNODE, STEP_S, Nprocs, Step2node,
373     &     IERR)
374      IMPLICIT NONE
375      INTEGER, INTENT(IN) :: SIZEPERM
376      INTEGER, intent(in) :: N, KEEP_28, Nprocs
377      INTEGER, intent(in) :: PROCNODE(KEEP_28), STEP_S(N)
378      INTEGER, intent(in) :: Step2node(KEEP_28)
379      INTEGER, INTENT(OUT):: IERR
380      INTEGER, INTENT(INOUT):: PERM_RHS(SIZEPERM)
381      INTEGER I, TMP_RHS, TMP2, proc_num
382      INTEGER , ALLOCATABLE :: TEMP_LOC_ARRAY(:)
383      INTEGER PTR(0:Nprocs-1)
384      INTEGER  MUMPS_275, MUMPS_330
385      EXTERNAL MUMPS_275, MUMPS_330
386      IERR = 0
387      ALLOCATE(TEMP_LOC_ARRAY(SIZEPERM), stat=IERR)
388      IF (IERR.GT.0) THEN
389       WRITE(6,*) " Not enough memory to allocate working ",
390     &            " arrays in MUMPS_772 "
391       CALL MUMPS_ABORT()
392      ENDIF
393      proc_num = 0
394      PTR(:) = 1
395      DO I = 1, SIZEPERM
396 555     CONTINUE
397         IF ( PTR(proc_num).LE.SIZEPERM ) THEN
398            TMP_RHS = PERM_RHS(PTR(proc_num))
399            TMP2 = Step2node(abs (STEP_S(TMP_RHS) ))
400              IF (proc_num .EQ. MUMPS_275
401     &           (PROCNODE(STEP_S(TMP2)),Nprocs)) THEN
402               TEMP_LOC_ARRAY(I) = TMP_RHS
403               PTR(proc_num) = PTR(proc_num)+1
404               IF ( (MUMPS_330(PROCNODE(STEP_S(TMP2)),
405     &              Nprocs).EQ.1)
406     &             ) THEN
407                  proc_num = mod(proc_num+1,Nprocs)
408                  proc_num = mod(proc_num+1,Nprocs)
409               ENDIF
410            ELSE
411               PTR(proc_num) = PTR(proc_num)+1
412               GOTO 555
413            ENDIF
414         ELSE
415            proc_num = mod(proc_num+1,Nprocs)
416            GOTO 555
417         ENDIF
418      ENDDO
419      WRITE(*,*) "Used interleaving of the RHS"
420      DO I = 1, SIZEPERM
421         PERM_RHS(I) = TEMP_LOC_ARRAY(I)
422      ENDDO
423      IF (allocated(TEMP_LOC_ARRAY)) DEALLOCATE (TEMP_LOC_ARRAY)
424      RETURN
425      END SUBROUTINE MUMPS_772
426