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