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
13C
14C History:
15C -------
16C This maximum transversal set of routines are
17C based on the work done by Jacko Koster at CERFACS for
18C his PhD thesis from Institut National Polytechnique de Toulouse
19C at CERFACS (1995-1997) and includes modifications provided
20C by the author as well as work done by Stephane Pralet
21C first at CERFACS during his PhD thesis (2003-2004) then
22C at INPT-IRIT (2004-2005) during his post-doctoral position.
23C
24C The main research publication references for this work are:
25C  [1] I. S. Duff, (1981),
26C      "Algorithm 575. Permutations for a zero-free diagonal",
27C      ACM Trans. Math. Software 7(3), 387-390.
28C  [2] I. S. Duff and J. Koster, (1998),
29C      "The design and use of algorithms for permuting large
30C      entries to the diagonal of sparse matrices",
31C      SIAM J. Matrix Anal. Appl., vol. 20, no. 4, pp. 889-901.
32C  [3] I. S. Duff and J. Koster, (2001),
33C      "On algorithms for permuting large entries to the diagonal
34C      of sparse matrices",
35C      SIAM J. Matrix Anal. Appl., vol. 22, no. 4, pp. 973-996.
36C
37      SUBROUTINE SMUMPS_MTRANSI(ICNTL,CNTL)
38      IMPLICIT NONE
39      INTEGER NICNTL, NCNTL
40      PARAMETER (NICNTL=10, NCNTL=10)
41      INTEGER ICNTL(NICNTL)
42      REAL CNTL(NCNTL)
43      INTEGER I
44      ICNTL(1) =  6
45      ICNTL(2) =  6
46      ICNTL(3) = -1
47      ICNTL(4) = -1
48      ICNTL(5) =  0
49      DO 10 I = 6,NICNTL
50        ICNTL(I) = 0
51   10 CONTINUE
52      CNTL(1) = 0.0E0
53      CNTL(2) = 0.0E0
54      DO 20 I = 3,NCNTL
55        CNTL(I) = 0.0E0
56   20 CONTINUE
57      RETURN
58      END SUBROUTINE SMUMPS_MTRANSI
59      SUBROUTINE SMUMPS_MTRANSB
60     &           (M,N,NE,IP,IRN,A,IPERM,NUM,JPERM,PR,Q,L,D,RINF)
61      IMPLICIT NONE
62      INTEGER ::  M,N,NUM
63      INTEGER(8), INTENT(IN)  :: NE
64      INTEGER :: IRN(NE),IPERM(M),JPERM(N),Q(M),L(M)
65      INTEGER(8), INTENT(IN)  :: IP(N+1)
66      INTEGER(8), INTENT(OUT) :: PR(N)
67      REAL  :: A(NE)
68      REAL  :: D(M), RINF
69      INTEGER    :: I,II,J,JJ,JORD,Q0,QLEN,IDUM,JDUM,ISP,JSP,
70     &              I0,UP,LOW, IK
71      INTEGER(8) :: K,KK,KK1,KK2
72      REAL    CSP,DI,DNEW,DQ0,AI,A0,BV,TBV,RLX
73      REAL    ZERO,MINONE,ONE
74      PARAMETER (ZERO=0.0E0,MINONE=-1.0E0,ONE=1.0E0)
75      INTRINSIC abs,min
76      EXTERNAL SMUMPS_MTRANSD, SMUMPS_MTRANSE,
77     &         SMUMPS_MTRANSF, SMUMPS_MTRANSX
78      RLX = D(1)
79      NUM = 0
80      BV = RINF
81      DO 10 I = 1,N
82        JPERM(I) = 0
83        PR(I) = IP(I)
84   10 CONTINUE
85      DO 12 I = 1,M
86        IPERM(I) = 0
87        D(I) = ZERO
88   12 CONTINUE
89      DO 30 J = 1,N
90        A0 = MINONE
91        DO 20 K = IP(J),IP(J+1)-1_8
92          I = IRN(K)
93          AI = abs(A(K))
94          IF (AI.GT.D(I)) D(I) = AI
95          IF (JPERM(J).NE.0) GO TO 20
96          IF (AI.GE.BV) THEN
97            A0 = BV
98            IF (IPERM(I).NE.0) GO TO 20
99            JPERM(J) = I
100            IPERM(I) = J
101            NUM = NUM + 1
102          ELSE
103            IF (AI.LE.A0) GO TO 20
104            A0 = AI
105            I0 = I
106          ENDIF
107   20   CONTINUE
108        IF (A0.NE.MINONE .AND. A0.LT.BV) THEN
109          BV = A0
110          IF (IPERM(I0).NE.0) GO TO 30
111          IPERM(I0) = J
112          JPERM(J) = I0
113          NUM = NUM + 1
114        ENDIF
115   30 CONTINUE
116      IF (M.EQ.N) THEN
117        DO 35 I = 1,M
118          BV = min(BV,D(I))
119   35   CONTINUE
120      ENDIF
121      IF (NUM.EQ.N) GO TO 1000
122      DO 95 J = 1,N
123        IF (JPERM(J).NE.0) GO TO 95
124        DO 50 K = IP(J),IP(J+1)-1_8
125          I = IRN(K)
126          AI = abs(A(K))
127          IF (AI.LT.BV) GO TO 50
128          IF (IPERM(I).EQ.0) GO TO 90
129          JJ = IPERM(I)
130          KK1 = PR(JJ)
131          KK2 = IP(JJ+1) - 1_8
132          IF (KK1.GT.KK2) GO TO 50
133          DO 70 KK = KK1,KK2
134            II = IRN(KK)
135            IF (IPERM(II).NE.0) GO TO 70
136            IF (abs(A(KK)).GE.BV) GO TO 80
137   70     CONTINUE
138          PR(JJ) = KK2 + 1_8
139   50   CONTINUE
140        GO TO 95
141   80   JPERM(JJ) = II
142        IPERM(II) = JJ
143        PR(JJ) = KK + 1_8
144   90   NUM = NUM + 1
145        JPERM(J) = I
146        IPERM(I) = J
147        PR(J) = K + 1_8
148   95 CONTINUE
149      IF (NUM.EQ.N) GO TO 1000
150      DO 99 I = 1,M
151        D(I) = MINONE
152        L(I) = 0
153   99 CONTINUE
154      TBV = BV * (ONE-RLX)
155      DO 100 JORD = 1,N
156        IF (JPERM(JORD).NE.0) GO TO 100
157        QLEN = 0
158        LOW = M + 1
159        UP = M + 1
160        CSP = MINONE
161        J = JORD
162        PR(J) = -1_8
163        DO 115 K = IP(J),IP(J+1)-1_8
164          I = IRN(K)
165          DNEW = abs(A(K))
166          IF (CSP.GE.DNEW) GO TO 115
167          IF (IPERM(I).EQ.0) THEN
168            CSP = DNEW
169            ISP = I
170            JSP = J
171            IF (CSP.GE.TBV) GO TO 160
172          ELSE
173            D(I) = DNEW
174            IF (DNEW.GE.TBV) THEN
175              LOW = LOW - 1
176              Q(LOW) = I
177            ELSE
178              QLEN = QLEN + 1
179              L(I) = QLEN
180              CALL SMUMPS_MTRANSD(I,M,Q,D,L,1)
181            ENDIF
182            JJ = IPERM(I)
183            PR(JJ) = int(J,8)
184          ENDIF
185  115   CONTINUE
186        DO 150 JDUM = 1,NUM
187          IF (LOW.EQ.UP) THEN
188            IF (QLEN.EQ.0) GO TO 160
189            I = Q(1)
190            IF (CSP.GE.D(I)) GO TO 160
191            BV = D(I)
192            TBV = BV * (ONE-RLX)
193            DO 152 IDUM = 1,M
194              CALL SMUMPS_MTRANSE(QLEN,M,Q,D,L,1)
195              L(I) = 0
196              LOW = LOW - 1
197              Q(LOW) = I
198              IF (QLEN.EQ.0) GO TO 153
199              I = Q(1)
200              IF (D(I).LT.TBV) GO TO 153
201  152       CONTINUE
202          ENDIF
203  153     UP = UP - 1
204          Q0 = Q(UP)
205          DQ0 = D(Q0)
206          L(Q0) = UP
207          J = IPERM(Q0)
208          DO 155 K = IP(J),IP(J+1)-1_8
209            I = IRN(K)
210            IF (L(I).GE.UP) GO TO 155
211            DNEW = min(DQ0,abs(A(K)))
212            IF (CSP.GE.DNEW) GO TO 155
213            IF (IPERM(I).EQ.0) THEN
214              CSP = DNEW
215              ISP = I
216              JSP = J
217              IF (CSP.GE.TBV) GO TO 160
218            ELSE
219              DI = D(I)
220              IF (DI.GE.TBV .OR. DI.GE.DNEW) GO TO 155
221              D(I) = DNEW
222              IF (DNEW.GE.TBV) THEN
223                IF (DI.NE.MINONE) THEN
224                  CALL SMUMPS_MTRANSF(L(I),QLEN,M,Q,D,L,1)
225                ENDIF
226                L(I) = 0
227                LOW = LOW - 1
228                Q(LOW) = I
229              ELSE
230                IF (DI.EQ.MINONE) THEN
231                  QLEN = QLEN + 1
232                  L(I) = QLEN
233                ENDIF
234                CALL SMUMPS_MTRANSD(I,M,Q,D,L,1)
235              ENDIF
236              JJ = IPERM(I)
237              PR(JJ) = int(J,8)
238            ENDIF
239  155     CONTINUE
240  150   CONTINUE
241  160   IF (CSP.EQ.MINONE) GO TO 190
242        BV = min(BV,CSP)
243        TBV = BV * (ONE-RLX)
244        NUM = NUM + 1
245        I = ISP
246        J = JSP
247        DO 170 JDUM = 1,NUM+1
248          I0 = JPERM(J)
249          JPERM(J) = I
250          IPERM(I) = J
251          J = int(PR(J))
252          IF (J.EQ.-1) GO TO 190
253          I = I0
254  170   CONTINUE
255  190   DO 191 IK = UP,M
256          I = Q(IK)
257          D(I) = MINONE
258          L(I) = 0
259  191   CONTINUE
260        DO 192 IK = LOW,UP-1
261          I = Q(IK)
262          D(I) = MINONE
263  192   CONTINUE
264        DO 193 IK = 1,QLEN
265          I = Q(IK)
266          D(I) = MINONE
267          L(I) = 0
268  193   CONTINUE
269  100 CONTINUE
270 1000 IF (M.EQ.N .and. NUM.EQ.N) GO TO 2000
271      CALL SMUMPS_MTRANSX(M,N,IPERM,L,JPERM)
272 2000 RETURN
273      END SUBROUTINE SMUMPS_MTRANSB
274      SUBROUTINE SMUMPS_MTRANSD(I,N,Q,D,L,IWAY)
275      IMPLICIT NONE
276      INTEGER I,N,IWAY
277      INTEGER Q(N),L(N)
278      REAL D(N)
279      INTEGER IDUM,K,POS,POSK,QK
280      PARAMETER (K=2)
281      REAL DI
282      POS = L(I)
283      IF (POS.LE.1) GO TO 20
284      DI = D(I)
285      IF (IWAY.EQ.1) THEN
286        DO 10 IDUM = 1,N
287          POSK = POS/K
288          QK = Q(POSK)
289          IF (DI.LE.D(QK)) GO TO 20
290          Q(POS) = QK
291          L(QK) = POS
292          POS = POSK
293          IF (POS.LE.1) GO TO 20
294   10   CONTINUE
295      ELSE
296        DO 15 IDUM = 1,N
297          POSK = POS/K
298          QK = Q(POSK)
299          IF (DI.GE.D(QK)) GO TO 20
300          Q(POS) = QK
301          L(QK) = POS
302          POS = POSK
303          IF (POS.LE.1) GO TO 20
304   15   CONTINUE
305      ENDIF
306   20 Q(POS) = I
307      L(I) = POS
308      RETURN
309      END SUBROUTINE SMUMPS_MTRANSD
310      SUBROUTINE SMUMPS_MTRANSE(QLEN,N,Q,D,L,IWAY)
311      IMPLICIT NONE
312      INTEGER QLEN,N,IWAY
313      INTEGER Q(N),L(N)
314      REAL D(N)
315      INTEGER I,IDUM,K,POS,POSK
316      PARAMETER (K=2)
317      REAL DK,DR,DI
318      I = Q(QLEN)
319      DI = D(I)
320      QLEN = QLEN - 1
321      POS = 1
322      IF (IWAY.EQ.1) THEN
323        DO 10 IDUM = 1,N
324          POSK = K*POS
325          IF (POSK.GT.QLEN) GO TO 20
326          DK = D(Q(POSK))
327          IF (POSK.LT.QLEN) THEN
328            DR = D(Q(POSK+1))
329            IF (DK.LT.DR) THEN
330              POSK = POSK + 1
331              DK = DR
332            ENDIF
333          ENDIF
334          IF (DI.GE.DK) GO TO 20
335          Q(POS) = Q(POSK)
336          L(Q(POS)) = POS
337          POS = POSK
338   10   CONTINUE
339      ELSE
340        DO 15 IDUM = 1,N
341          POSK = K*POS
342          IF (POSK.GT.QLEN) GO TO 20
343          DK = D(Q(POSK))
344          IF (POSK.LT.QLEN) THEN
345            DR = D(Q(POSK+1))
346            IF (DK.GT.DR) THEN
347              POSK = POSK + 1
348              DK = DR
349            ENDIF
350          ENDIF
351          IF (DI.LE.DK) GO TO 20
352          Q(POS) = Q(POSK)
353          L(Q(POS)) = POS
354          POS = POSK
355   15   CONTINUE
356      ENDIF
357   20 Q(POS) = I
358      L(I) = POS
359      RETURN
360      END SUBROUTINE SMUMPS_MTRANSE
361      SUBROUTINE SMUMPS_MTRANSF(POS0,QLEN,N,Q,D,L,IWAY)
362      IMPLICIT NONE
363      INTEGER POS0,QLEN,N,IWAY
364      INTEGER Q(N),L(N)
365      REAL D(N)
366      INTEGER I,IDUM,K,POS,POSK,QK
367      PARAMETER (K=2)
368      REAL DK,DR,DI
369      IF (QLEN.EQ.POS0) THEN
370        QLEN = QLEN - 1
371        RETURN
372      ENDIF
373      I = Q(QLEN)
374      DI = D(I)
375      QLEN = QLEN - 1
376      POS = POS0
377      IF (IWAY.EQ.1) THEN
378        IF (POS.LE.1) GO TO 20
379        DO 10 IDUM = 1,N
380          POSK = POS/K
381          QK = Q(POSK)
382          IF (DI.LE.D(QK)) GO TO 20
383          Q(POS) = QK
384          L(QK) = POS
385          POS = POSK
386          IF (POS.LE.1) GO TO 20
387   10   CONTINUE
388   20   Q(POS) = I
389        L(I) = POS
390        IF (POS.NE.POS0) RETURN
391        DO 30 IDUM = 1,N
392          POSK = K*POS
393          IF (POSK.GT.QLEN) GO TO 40
394          DK = D(Q(POSK))
395          IF (POSK.LT.QLEN) THEN
396            DR = D(Q(POSK+1))
397            IF (DK.LT.DR) THEN
398              POSK = POSK + 1
399              DK = DR
400            ENDIF
401          ENDIF
402          IF (DI.GE.DK) GO TO 40
403          QK = Q(POSK)
404          Q(POS) = QK
405          L(QK) = POS
406          POS = POSK
407   30   CONTINUE
408      ELSE
409        IF (POS.LE.1) GO TO 34
410        DO 32 IDUM = 1,N
411          POSK = POS/K
412          QK = Q(POSK)
413          IF (DI.GE.D(QK)) GO TO 34
414          Q(POS) = QK
415          L(QK) = POS
416          POS = POSK
417          IF (POS.LE.1) GO TO 34
418   32   CONTINUE
419   34   Q(POS) = I
420        L(I) = POS
421        IF (POS.NE.POS0) RETURN
422        DO 36 IDUM = 1,N
423          POSK = K*POS
424          IF (POSK.GT.QLEN) GO TO 40
425          DK = D(Q(POSK))
426          IF (POSK.LT.QLEN) THEN
427            DR = D(Q(POSK+1))
428            IF (DK.GT.DR) THEN
429              POSK = POSK + 1
430              DK = DR
431            ENDIF
432          ENDIF
433          IF (DI.LE.DK) GO TO 40
434          QK = Q(POSK)
435          Q(POS) = QK
436          L(QK) = POS
437          POS = POSK
438   36   CONTINUE
439      ENDIF
440   40 Q(POS) = I
441      L(I) = POS
442      RETURN
443      END SUBROUTINE SMUMPS_MTRANSF
444      SUBROUTINE SMUMPS_MTRANSQ(IP,LENL,LENH,W,WLEN,A,NVAL,VAL)
445      IMPLICIT NONE
446      INTEGER ::WLEN,NVAL
447      INTEGER :: LENL(*),LENH(*),W(*)
448      INTEGER(8) :: IP(*)
449      REAL :: A(*),VAL
450      INTEGER XX,J,K,S,POS
451      INTEGER(8) :: II
452      PARAMETER (XX=10)
453      REAL SPLIT(XX),HA
454      NVAL = 0
455      DO 10 K = 1,WLEN
456        J = W(K)
457        DO 15 II = IP(J)+int(LENL(J),8),IP(J)+int(LENH(J)-1,8)
458          HA = A(II)
459          IF (NVAL.EQ.0) THEN
460            SPLIT(1) = HA
461            NVAL = 1
462          ELSE
463            DO 20 S = NVAL,1,-1
464              IF (SPLIT(S).EQ.HA) GO TO 15
465              IF (SPLIT(S).GT.HA) THEN
466                POS = S + 1
467                GO TO 21
468              ENDIF
469  20        CONTINUE
470            POS = 1
471  21        DO 22 S = NVAL,POS,-1
472              SPLIT(S+1) = SPLIT(S)
473  22        CONTINUE
474            SPLIT(POS) = HA
475            NVAL = NVAL + 1
476          ENDIF
477          IF (NVAL.EQ.XX) GO TO 11
478  15    CONTINUE
479  10  CONTINUE
480  11  IF (NVAL.GT.0) VAL = SPLIT((NVAL+1)/2)
481      RETURN
482      END SUBROUTINE SMUMPS_MTRANSQ
483      SUBROUTINE SMUMPS_MTRANSR(N,NE,IP,IRN,A)
484      IMPLICIT NONE
485      INTEGER, INTENT(IN)    :: N
486      INTEGER(8), INTENT(IN) :: NE
487      INTEGER(8), INTENT(IN) :: IP(N+1)
488      INTEGER, INTENT(INOUT) :: IRN(NE)
489      REAL, INTENT(INOUT)    :: A(NE)
490      INTEGER  :: THRESH,TDLEN
491      PARAMETER (THRESH=15,TDLEN=50)
492      INTEGER    :: J, LEN, HI
493      INTEGER(8) :: K, IPJ, TD, FIRST, LAST, MID, R, S
494      REAL       :: HA, KEY
495      INTEGER(8) :: TODO(TDLEN)
496      DO 100 J = 1,N
497        LEN = int(IP(J+1) - IP(J))
498        IF (LEN.LE.1) GO TO 100
499        IPJ = IP(J)
500        IF (LEN.LT.THRESH) GO TO 400
501        TODO(1) = IPJ
502        TODO(2) = IPJ +int(LEN,8)
503        TD = 2_8
504  500   CONTINUE
505        FIRST = TODO(TD-1)
506        LAST = TODO(TD)
507        KEY = A((FIRST+LAST)/2)
508        DO 475 K = FIRST,LAST-1
509          HA = A(K)
510          IF (HA.EQ.KEY) GO TO 475
511          IF (HA.GT.KEY) GO TO 470
512          KEY = HA
513          GO TO 470
514  475   CONTINUE
515        TD = TD - 2_8
516        GO TO 425
517  470   MID = FIRST
518        DO 450 K = FIRST,LAST-1
519          IF (A(K).LE.KEY) GO TO 450
520          HA = A(MID)
521          A(MID) = A(K)
522          A(K) = HA
523          HI = IRN(MID)
524          IRN(MID) = IRN(K)
525          IRN(K) = HI
526          MID = MID + 1
527  450   CONTINUE
528        IF (MID-FIRST.GE.LAST-MID) THEN
529          TODO(TD+2) = LAST
530          TODO(TD+1) = MID
531          TODO(TD) = MID
532        ELSE
533          TODO(TD+2) = MID
534          TODO(TD+1) = FIRST
535          TODO(TD) = LAST
536          TODO(TD-1) = MID
537        ENDIF
538        TD = TD + 2_8
539  425   CONTINUE
540        IF (TD.EQ.0_8) GO TO 400
541        IF (TODO(TD)-TODO(TD-1).GE.int(THRESH,8)) GO TO 500
542        TD = TD - 2_8
543        GO TO 425
544  400   DO 200 R = IPJ+1_8,IPJ+int(LEN-1,8)
545          IF (A(R-1) .LT. A(R)) THEN
546            HA = A(R)
547            HI = IRN(R)
548            A(R) = A(R-1_8)
549            IRN(R) = IRN(R-1_8)
550            DO 300 S = R-1,IPJ+1_8,-1_8
551              IF (A(S-1) .LT. HA) THEN
552                A(S) = A(S-1)
553                IRN(S) = IRN(S-1)
554              ELSE
555                A(S) = HA
556                IRN(S) = HI
557                GO TO 200
558              END IF
559  300       CONTINUE
560            A(IPJ) = HA
561            IRN(IPJ) = HI
562          END IF
563  200   CONTINUE
564  100 CONTINUE
565      RETURN
566      END SUBROUTINE SMUMPS_MTRANSR
567      SUBROUTINE SMUMPS_MTRANSS(M,N,NE,IP,IRN,A,IPERM,NUMX,
568     &           W,LEN,LENL,LENH,FC,IW,IW4,RLX,RINF)
569      IMPLICIT NONE
570      INTEGER, INTENT(IN) :: M,N
571      INTEGER(8), INTENT(IN) :: NE
572      INTEGER, INTENT(OUT) :: NUMX
573      INTEGER(8), INTENT(IN) :: IP(N+1)
574      INTEGER  :: IRN(NE),IPERM(N),
575     &        W(N),LEN(N),LENL(N),LENH(N),FC(N),IW(M),IW4(3*N+M)
576      REAL A(NE),RLX,RINF
577      INTEGER :: NUM,NVAL,WLEN,I,J,L,CNT,MOD, IDUM
578      INTEGER(8) :: K, II, KDUM1, KDUM2
579      REAL ::    BVAL,BMIN,BMAX
580      EXTERNAL SMUMPS_MTRANSQ,SMUMPS_MTRANSU,SMUMPS_MTRANSX
581      DO 20 J = 1,N
582        FC(J) = J
583        LEN(J) = int(IP(J+1) - IP(J))
584   20 CONTINUE
585      DO 21 I = 1,M
586        IW(I) = 0
587   21 CONTINUE
588      CNT = 1
589      MOD = 1
590      NUMX = 0
591      CALL SMUMPS_MTRANSU(CNT,MOD,M,N,IRN,NE,IP,LEN,FC,IW,
592     &            NUMX,N,
593     &            IW4(1),IW4(N+1),IW4(2*N+1),IW4(2*N+M+1))
594      NUM = NUMX
595      IF (NUM.NE.N) THEN
596        BMAX = RINF
597      ELSE
598        BMAX = RINF
599        DO 30 J = 1,N
600          BVAL = 0.0E0
601          DO 25 K = IP(J),IP(J+1)-1_8
602            IF (A(K).GT.BVAL) BVAL = A(K)
603   25     CONTINUE
604          IF (BVAL.LT.BMAX) BMAX = BVAL
605   30   CONTINUE
606        BMAX = 1.001E0 * BMAX
607      ENDIF
608      BVAL = 0.0E0
609      BMIN = 0.0E0
610      WLEN = 0
611      DO 48 J = 1,N
612        L = int(IP(J+1) - IP(J))
613        LENH(J) = L
614        LEN(J) = L
615        DO 45 K = IP(J),IP(J+1)-1_8
616          IF (A(K).LT.BMAX) GO TO 46
617   45   CONTINUE
618        K = IP(J+1)
619   46   LENL(J) = int(K - IP(J))
620        IF (LENL(J).EQ.L) GO TO 48
621        WLEN = WLEN + 1
622        W(WLEN) = J
623   48 CONTINUE
624      DO 90 KDUM1 = 1_8,NE
625        IF (NUM.EQ.NUMX) THEN
626          DO 50 I = 1,M
627            IPERM(I) = IW(I)
628   50     CONTINUE
629          DO 80 KDUM2 = 1_8,NE
630            BMIN = BVAL
631            IF (BMAX-BMIN .LE. RLX) GO TO 1000
632            CALL SMUMPS_MTRANSQ(IP,LENL,LEN,W,WLEN,A,NVAL,BVAL)
633            IF (NVAL.LE.1) GO TO 1000
634            K = 1
635            DO 70 IDUM = 1,N
636              IF (K.GT.WLEN) GO TO 71
637              J = W(K)
638              DO 55 II = IP(J)+int(LEN(J)-1,8),
639     &                   IP(J)+int(LENL(J),8),-1_8
640                IF (A(II).GE.BVAL) GO TO 60
641                I = IRN(II)
642                IF (IW(I).NE.J) GO TO 55
643                IW(I) = 0
644                NUM = NUM - 1
645                FC(N-NUM) = J
646   55         CONTINUE
647   60         LENH(J) = LEN(J)
648              LEN(J) = int(II - IP(J) + 1)
649              IF (LENL(J).EQ.LENH(J)) THEN
650                W(K) = W(WLEN)
651                WLEN = WLEN - 1
652              ELSE
653                K = K + 1
654              ENDIF
655   70       CONTINUE
656   71       IF (NUM.LT.NUMX) GO TO 81
657   80     CONTINUE
658   81     MOD = 1
659        ELSE
660          BMAX = BVAL
661          IF (BMAX-BMIN .LE. RLX) GO TO 1000
662          CALL SMUMPS_MTRANSQ(IP,LEN,LENH,W,WLEN,A,NVAL,BVAL)
663          IF (NVAL.EQ.0. OR. BVAL.EQ.BMIN) GO TO 1000
664          K = 1
665          DO 87 IDUM = 1,N
666            IF (K.GT.WLEN) GO TO 88
667            J = W(K)
668            DO 85 II = IP(J)+int(LEN(J),8),IP(J)+int(LENH(J)-1,8)
669              IF (A(II).LT.BVAL) GO TO 86
670   85       CONTINUE
671   86       LENL(J) = LEN(J)
672            LEN(J) = int(II - IP(J))
673            IF (LENL(J).EQ.LENH(J)) THEN
674              W(K) = W(WLEN)
675              WLEN = WLEN - 1
676            ELSE
677              K = K + 1
678            ENDIF
679   87     CONTINUE
680   88     MOD = 0
681        ENDIF
682        CNT = CNT + 1
683        CALL SMUMPS_MTRANSU(CNT,MOD,M,N,IRN,NE,IP,LEN,FC,IW,
684     &              NUM,NUMX,
685     &              IW4(1),IW4(N+1),IW4(2*N+1),IW4(2*N+M+1))
686   90 CONTINUE
687 1000 IF (M.EQ.N .and. NUMX.EQ.N) GO TO 2000
688      CALL SMUMPS_MTRANSX(M,N,IPERM,IW,W)
689 2000 RETURN
690      END SUBROUTINE SMUMPS_MTRANSS
691C
692      SUBROUTINE SMUMPS_MTRANSU
693     &           (ID,MOD,M,N,IRN,LIRN,IP,LENC,FC,IPERM,NUM,NUMX,
694     &           PR,ARP,CV,OUT)
695      IMPLICIT NONE
696      INTEGER  :: ID,MOD,M,N,NUM,NUMX
697      INTEGER(8), INTENT(IN) :: LIRN
698      INTEGER  :: ARP(N),CV(M),IRN(LIRN),
699     &        FC(N),IPERM(M),LENC(N),OUT(N),PR(N)
700      INTEGER(8), INTENT(IN) :: IP(N)
701      INTEGER I,J,J1,JORD,NFC,K,KK,
702     &        NUM0,NUM1,NUM2,ID0,ID1,LAST
703      INTEGER(8) :: IN1, IN2, II
704      IF (ID.EQ.1) THEN
705        DO 5 I = 1,M
706          CV(I) = 0
707    5   CONTINUE
708        DO 6 J = 1,N
709          ARP(J) = 0
710    6   CONTINUE
711        NUM1 = N
712        NUM2 = N
713      ELSE
714        IF (MOD.EQ.1) THEN
715          DO 8 J = 1,N
716            ARP(J) = 0
717    8     CONTINUE
718        ENDIF
719        NUM1 = NUMX
720        NUM2 = N - NUMX
721      ENDIF
722      NUM0 = NUM
723      NFC = 0
724      ID0 = (ID-1)*N
725      DO 100 JORD = NUM0+1,N
726        ID1 = ID0 + JORD
727        J = FC(JORD-NUM0)
728        PR(J) = -1
729        DO 70 K = 1,JORD
730          IF (ARP(J).GE.LENC(J)) GO TO 30
731          IN1 = IP(J) + int(ARP(J),8)
732          IN2 = IP(J) + int(LENC(J) - 1,8)
733          DO 20 II = IN1,IN2
734            I = IRN(II)
735            IF (IPERM(I).EQ.0) GO TO 80
736   20     CONTINUE
737          ARP(J) = LENC(J)
738   30     OUT(J) = LENC(J) - 1
739          DO 60 KK = 1,JORD
740            IN1 = int(OUT(J),8)
741            IF (IN1.LT.0) GO TO 50
742            IN2 = IP(J) + int(LENC(J) - 1,8)
743            IN1 = IN2 - IN1
744            DO 40 II = IN1,IN2
745              I = IRN(II)
746              IF (CV(I).EQ.ID1) GO TO 40
747              J1 = J
748              J = IPERM(I)
749              CV(I) = ID1
750              PR(J) = J1
751              OUT(J1) = int(IN2 - II) - 1
752              GO TO 70
753   40       CONTINUE
754   50       J1 = PR(J)
755            IF (J1.EQ.-1) THEN
756              NFC = NFC + 1
757              FC(NFC) = J
758              IF (NFC.GT.NUM2) THEN
759                LAST = JORD
760                GO TO 101
761              ENDIF
762              GO TO 100
763            ENDIF
764            J = J1
765   60     CONTINUE
766   70   CONTINUE
767   80   IPERM(I) = J
768        ARP(J) = int(II - IP(J)) + 1
769        NUM = NUM + 1
770        DO 90 K = 1,JORD
771          J = PR(J)
772          IF (J.EQ.-1) GO TO 95
773          II = IP(J) + int(LENC(J) - OUT(J) - 2,8)
774          I = IRN(II)
775          IPERM(I) = J
776   90   CONTINUE
777   95   IF (NUM.EQ.NUM1) THEN
778          LAST = JORD
779          GO TO 101
780        ENDIF
781  100 CONTINUE
782      LAST = N
783  101 DO 110 JORD = LAST+1,N
784        NFC = NFC + 1
785        FC(NFC) = FC(JORD-NUM0)
786  110 CONTINUE
787      RETURN
788      END SUBROUTINE SMUMPS_MTRANSU
789C
790      SUBROUTINE SMUMPS_MTRANSW(M,N,NE,IP,IRN,A,IPERM,NUM,
791     &           JPERM,L32,OUT,PR,Q,L,U,D,RINF)
792      IMPLICIT NONE
793      INTEGER :: M,N,NUM
794      INTEGER(8), INTENT(IN) :: NE
795      INTEGER  :: IRN(NE),IPERM(M),Q(M),L32(max(M,N))
796      INTEGER(8)  :: IP(N+1), PR(N), L(M), JPERM(N), OUT(N)
797      REAL A(NE),U(M),D(M),RINF,RINF3
798      INTEGER  :: I,I0,II,J,JJ,JORD,Q0,QLEN,JDUM,JSP,
799     &            UP,LOW,IK
800      INTEGER(8) :: K, KK, KK1, KK2, K0, K1, K2, ISP
801      REAL     :: CSP,DI,DMIN,DNEW,DQ0,VJ,RLX
802      LOGICAL  :: LORD
803      REAL    :: ZERO, ONE
804      PARAMETER (ZERO=0.0E0,ONE=1.0E0)
805      EXTERNAL SMUMPS_MTRANSD, SMUMPS_MTRANSE,
806     &         SMUMPS_MTRANSF, SMUMPS_MTRANSX
807      RLX = U(1)
808      RINF3 = U(2)
809      LORD = (JPERM(1).EQ.6)
810      NUM = 0
811      DO 10 I = 1,N
812        JPERM(I) = 0_8
813        PR(I) = IP(I)
814        D(I) = RINF
815   10 CONTINUE
816      DO 15 I = 1,M
817        U(I) = RINF3
818        IPERM(I) = 0
819        L(I) = 0_8
820   15 CONTINUE
821      DO 30 J = 1,N
822         IF (int(IP(J+1)-IP(J)) .GT. N/10 .AND. N.GT.50) GO TO 30
823        DO 20 K = IP(J),IP(J+1)-1
824          I = IRN(K)
825          IF (A(K).GT.U(I)) GO TO 20
826          U(I) = A(K)
827          IPERM(I) = J
828          L(I) = K
829   20   CONTINUE
830   30 CONTINUE
831      DO 40 I = 1,M
832        J = IPERM(I)
833        IF (J.EQ.0) GO TO 40
834        IF (JPERM(J).EQ.0_8) THEN
835          JPERM(J) = L(I)
836          D(J) = U(I)
837          NUM = NUM + 1
838        ELSEIF (D(J).GT.U(I)) THEN
839          K = JPERM(J)
840          II = IRN(K)
841          IPERM(II) = 0
842          JPERM(J) = L(I)
843          D(J) = U(I)
844        ELSE
845          IPERM(I) = 0
846        ENDIF
847   40 CONTINUE
848      IF (NUM.EQ.N) GO TO 1000
849      DO 45 I = 1,M
850        D(I) = ZERO
851   45 CONTINUE
852      DO 95 J = 1,N
853        IF (JPERM(J).NE.0) GO TO 95
854        K1 = IP(J)
855        K2 = IP(J+1) - 1_8
856        IF (K1.GT.K2) GO TO 95
857        VJ = RINF
858        DO 50 K = K1,K2
859          I = IRN(K)
860          DI = A(K) - U(I)
861          IF (DI.GT.VJ) GO TO 50
862          IF (DI.LT.VJ .OR. DI.EQ.RINF) GO TO 55
863          IF (IPERM(I).NE.0 .OR. IPERM(I0).EQ.0) GO TO 50
864   55     VJ = DI
865          I0 = I
866          K0 = K
867   50   CONTINUE
868        D(J) = VJ
869        K = K0
870        I = I0
871        IF (IPERM(I).EQ.0) GO TO 90
872        DO 60 K = K0,K2
873          I = IRN(K)
874          IF (A(K)-U(I).GT.VJ) GO TO 60
875          JJ = IPERM(I)
876          KK1 = PR(JJ)
877          KK2 = IP(JJ+1) - 1_8
878          IF (KK1.GT.KK2) GO TO 60
879          DO 70 KK = KK1,KK2
880            II = IRN(KK)
881            IF (IPERM(II).GT.0) GO TO 70
882            IF (A(KK)-U(II).LE.D(JJ)) GO TO 80
883   70     CONTINUE
884          PR(JJ) = KK2 + 1_8
885   60   CONTINUE
886        GO TO 95
887   80   JPERM(JJ) = KK
888        IPERM(II) = JJ
889        PR(JJ) = KK + 1_8
890   90   NUM = NUM + 1
891        JPERM(J) = K
892        IPERM(I) = J
893        PR(J) = K + 1_8
894   95 CONTINUE
895      IF (NUM.EQ.N) GO TO 1000
896      DO 99 I = 1,M
897        D(I) = RINF
898        Q(I) = 0
899   99 CONTINUE
900      DO 100 JORD = 1,N
901        IF (JPERM(JORD).NE.0) GO TO 100
902        DMIN = RINF
903        QLEN = 0
904        LOW = M + 1
905        UP = M + 1
906        CSP = RINF
907        J = JORD
908        PR(J) = -1_8
909        DO 115 K = IP(J),IP(J+1)-1_8
910          I = IRN(K)
911          DNEW = A(K) - U(I)
912          IF (DNEW.GE.CSP) GO TO 115
913          IF (IPERM(I).EQ.0) THEN
914            CSP = DNEW
915            ISP = K
916            JSP = J
917          ELSE
918            IF (DNEW.LT.DMIN) DMIN = DNEW
919            D(I) = DNEW
920            QLEN = QLEN + 1
921            L(QLEN) = K
922          ENDIF
923  115   CONTINUE
924        Q0 = QLEN
925        QLEN = 0
926        DO 120 IK = 1,Q0
927          K = L(IK)
928          I = IRN(K)
929          IF (CSP.LE.D(I)) THEN
930            D(I) = RINF
931            GO TO 120
932          ENDIF
933          IF (D(I).LE.DMIN) THEN
934            LOW = LOW - 1
935            L32(LOW) = I
936            Q(I) = LOW
937          ELSE
938            QLEN = QLEN + 1
939            Q(I) = QLEN
940            CALL SMUMPS_MTRANSD(I,M,L32,D,Q,2)
941          ENDIF
942          JJ = IPERM(I)
943          OUT(JJ) = K
944          PR(JJ) = int(J,8)
945  120   CONTINUE
946        DO 150 JDUM = 1,NUM
947          IF (LOW.EQ.UP) THEN
948            IF (QLEN.EQ.0) GO TO 160
949            I = L32(1)
950            IF (D(I).LT.RINF) DMIN = D(I)*(ONE+RLX)
951            IF (DMIN.GE.CSP) GO TO 160
952  152       CALL SMUMPS_MTRANSE(QLEN,M,L32,D,Q,2)
953            LOW = LOW - 1
954            L32(LOW) = I
955            Q(I) = LOW
956            IF (QLEN.EQ.0) GO TO 153
957            I = L32(1)
958            IF (D(I).GT.DMIN) GO TO 153
959            GO TO 152
960          ENDIF
961  153     Q0 = L32(UP-1)
962          DQ0 = D(Q0)
963          IF (DQ0.GE.CSP) GO TO 160
964          IF (DMIN.GE.CSP) GO TO 160
965          UP = UP - 1
966          J = IPERM(Q0)
967          VJ = DQ0 - A(JPERM(J)) + U(Q0)
968          K1 = IP(J+1)-1_8
969          IF (LORD) THEN
970            IF (CSP.NE.RINF) THEN
971              DI = CSP - VJ
972              IF (A(K1).GE.DI) THEN
973                K0 = JPERM(J)
974                IF (K0.GE.K1-6) GO TO 178
975  177           CONTINUE
976                  K = (K0+K1)/2
977                  IF (A(K).GE.DI) THEN
978                    K1 = K
979                  ELSE
980                    K0 = K
981                  ENDIF
982                  IF (K0.GE.K1-6) GO TO 178
983                GO TO 177
984  178           DO 179 K = K0+1,K1
985                  IF (A(K).LT.DI) GO TO 179
986                  K1 = K - 1
987                  GO TO 181
988  179           CONTINUE
989              ENDIF
990            ENDIF
991  181       IF (K1.EQ.JPERM(J)) K1 = K1 - 1
992          ENDIF
993          K0 = IP(J)
994          DI = CSP - VJ
995          DO 155 K = K0,K1
996            I = IRN(K)
997            IF (Q(I).GE.LOW) GO TO 155
998            DNEW = A(K) - U(I)
999            IF (DNEW.GE.DI) GO TO 155
1000            DNEW = DNEW + VJ
1001            IF (DNEW.GT.D(I)) GO TO 155
1002            IF (IPERM(I).EQ.0) THEN
1003              CSP = DNEW
1004              ISP = K
1005              JSP = J
1006              DI = CSP - VJ
1007            ELSE
1008              IF (DNEW.GE.D(I)) GO TO 155
1009              D(I) = DNEW
1010              IF (DNEW.LE.DMIN) THEN
1011                IF (Q(I).NE.0) THEN
1012                  CALL SMUMPS_MTRANSF(Q(I),QLEN,M,L32,D,Q,2)
1013                ENDIF
1014                LOW = LOW - 1
1015                L32(LOW) = I
1016                Q(I) = LOW
1017              ELSE
1018                IF (Q(I).EQ.0) THEN
1019                  QLEN = QLEN + 1
1020                  Q(I) = QLEN
1021                ENDIF
1022                CALL SMUMPS_MTRANSD(I,M,L32,D,Q,2)
1023              ENDIF
1024              JJ = IPERM(I)
1025              OUT(JJ) = K
1026              PR(JJ) = int(J,8)
1027            ENDIF
1028  155     CONTINUE
1029  150   CONTINUE
1030  160   IF (CSP.EQ.RINF) GO TO 190
1031        NUM = NUM + 1
1032        I = IRN(ISP)
1033        J = JSP
1034        IPERM(I) = J
1035        JPERM(J) = ISP
1036        DO 170 JDUM = 1,NUM
1037          JJ = int(PR(J))
1038          IF (JJ.EQ.-1) GO TO 180
1039          K = OUT(J)
1040          I = IRN(K)
1041          IPERM(I) = JJ
1042          JPERM(JJ) = K
1043          J = JJ
1044  170   CONTINUE
1045  180   DO 182 JJ = UP,M
1046          I = L32(JJ)
1047          U(I) = U(I) + D(I) - CSP
1048  182   CONTINUE
1049  190   DO 191 JJ = UP,M
1050          I = L32(JJ)
1051          D(I) = RINF
1052          Q(I) = 0
1053  191   CONTINUE
1054        DO 192 JJ = LOW,UP-1
1055          I = L32(JJ)
1056          D(I) = RINF
1057          Q(I) = 0
1058  192   CONTINUE
1059        DO 193 JJ = 1,QLEN
1060          I = L32(JJ)
1061          D(I) = RINF
1062          Q(I) = 0
1063  193   CONTINUE
1064  100 CONTINUE
1065 1000 CONTINUE
1066      DO 1200 J = 1,N
1067        K = JPERM(J)
1068        IF (K.NE.0) THEN
1069          D(J) = A(K) - U(IRN(K))
1070        ELSE
1071          D(J) = ZERO
1072        ENDIF
1073 1200 CONTINUE
1074      DO 1201 I = 1,M
1075        IF (IPERM(I).EQ.0) U(I) = ZERO
1076 1201 CONTINUE
1077      IF (M.EQ.N .and. NUM.EQ.N) GO TO 2000
1078      CALL SMUMPS_MTRANSX(M,N,IPERM,Q,L32)
1079 2000 RETURN
1080      END SUBROUTINE SMUMPS_MTRANSW
1081      SUBROUTINE SMUMPS_MTRANSZ
1082     &           (M,N,IRN,LIRN,IP,LENC,IPERM,NUM,PR,ARP,CV,OUT)
1083      IMPLICIT NONE
1084      INTEGER :: M,N,NUM
1085      INTEGER(8), INTENT(IN) :: LIRN
1086      INTEGER :: ARP(N),CV(M),IRN(LIRN),IPERM(M),LENC(N),OUT(N),PR(N)
1087      INTEGER(8), INTENT(IN) :: IP(N)
1088C Local variables
1089      INTEGER    :: I,J,J1,JORD,K,KK
1090      INTEGER(8) :: II, IN1, IN2
1091      EXTERNAL SMUMPS_MTRANSX
1092      DO 10 I = 1,M
1093        CV(I) = 0
1094        IPERM(I) = 0
1095   10 CONTINUE
1096      DO 12 J = 1,N
1097        ARP(J) = LENC(J) - 1
1098   12 CONTINUE
1099      NUM = 0
1100      DO 1000 JORD = 1,N
1101        J = JORD
1102        PR(J) = -1
1103        DO 70 K = 1,JORD
1104          IN1 = int(ARP(J),8)
1105          IF (IN1.LT.0_8) GO TO 30
1106          IN2 = IP(J) + int(LENC(J) - 1,8)
1107          IN1 = IN2 - IN1
1108          DO 20 II = IN1,IN2
1109            I = IRN(II)
1110            IF (IPERM(I).EQ.0) GO TO 80
1111   20     CONTINUE
1112          ARP(J) = -1
1113   30     CONTINUE
1114          OUT(J) = LENC(J) - 1
1115          DO 60 KK = 1,JORD
1116            IN1 = int(OUT(J),8)
1117            IF (IN1.LT.0_8) GO TO 50
1118            IN2 = IP(J) + int(LENC(J) - 1,8)
1119            IN1 = IN2 - IN1
1120            DO 40 II = IN1,IN2
1121              I = IRN(II)
1122              IF (CV(I).EQ.JORD) GO TO 40
1123              J1 = J
1124              J = IPERM(I)
1125              CV(I) = JORD
1126              PR(J) = J1
1127              OUT(J1) = int(IN2 - II - 1_8)
1128              GO TO 70
1129   40       CONTINUE
1130   50       CONTINUE
1131            J = PR(J)
1132            IF (J.EQ.-1) GO TO 1000
1133   60     CONTINUE
1134   70   CONTINUE
1135   80   CONTINUE
1136        IPERM(I) = J
1137        ARP(J) = int(IN2 - II - 1_8)
1138        NUM = NUM + 1
1139        DO 90 K = 1,JORD
1140          J = PR(J)
1141          IF (J.EQ.-1) GO TO 1000
1142          II = IP(J) + int(LENC(J) - OUT(J) - 2,8)
1143          I = IRN(II)
1144          IPERM(I) = J
1145   90   CONTINUE
1146 1000 CONTINUE
1147      IF (M.EQ.N .and. NUM.EQ.N) GO TO 2000
1148      CALL SMUMPS_MTRANSX(M,N,IPERM,CV,ARP)
1149 2000 RETURN
1150      END SUBROUTINE SMUMPS_MTRANSZ
1151      SUBROUTINE SMUMPS_MTRANSX(M,N,IPERM,RW,CW)
1152      IMPLICIT NONE
1153      INTEGER M,N
1154      INTEGER RW(M),CW(N),IPERM(M)
1155      INTEGER I,J,K
1156      DO 10 J = 1,N
1157        CW(J) = 0
1158   10 CONTINUE
1159      K = 0
1160      DO 20 I = 1,M
1161        IF (IPERM(I).EQ.0) THEN
1162          K = K + 1
1163          RW(K) = I
1164        ELSE
1165          J = IPERM(I)
1166          CW(J) = I
1167        ENDIF
1168   20 CONTINUE
1169      K = 0
1170      DO 30 J = 1,N
1171        IF (CW(J).NE.0) GO TO 30
1172        K = K + 1
1173        I = RW(K)
1174        IPERM(I) = -J
1175   30 CONTINUE
1176      DO 40 J = N+1,M
1177        K = K + 1
1178        I = RW(K)
1179        IPERM(I) = -J
1180   40 CONTINUE
1181      RETURN
1182      END SUBROUTINE SMUMPS_MTRANSX
1183