OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
ana_AMDMF.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
14 SUBROUTINE mumps_symqamd_new
15 & ( job, thresh, ndense,
16 & n, iwlen, pe, pfree, len, iw, nv,
17 & elen, last, ncmpa, degree, head, next, w,
18 & perm, complem_list, size_complem_list,
19 & agg6 )
20 IMPLICIT NONE
21 INTEGER, INTENT(IN) :: N, SIZE_COMPLEM_LIST
22 INTEGER(8), INTENT(IN) :: IWLEN
23 INTEGER, INTENT(IN) :: THRESH
24 LOGICAL, INTENT(IN) :: AGG6
25 INTEGER, INTENT (IN) :: COMPLEM_LIST(max(1,SIZE_COMPLEM_LIST))
26 INTEGER, INTENT(INOUT) :: JOB
27 INTEGER, INTENT(INOUT) :: LEN(N), IW(IWLEN)
28 INTEGER(8), INTENT(INOUT) :: PFREE
29 INTEGER(8), INTENT(INOUT) :: PE(N)
30 INTEGER, INTENT(INOUT) :: PERM(N)
31 INTEGER, INTENT(OUT) :: NCMPA
32 INTEGER, INTENT(OUT) :: NV(N), LAST(N)
33 INTEGER, INTENT(INOUT) :: ELEN(N)
34 INTEGER, INTENT(OUT) :: NDENSE(N), DEGREE(N),
35 & HEAD(N), NEXT(N), W(N)
36 INTEGER THRESM, NDME, PERMeqN
37 INTEGER NBD,NBED, NBDM, LASTD, NELME
38 LOGICAL IDENSE
39 INTEGER :: FDEG, ThresMin, ThresPrev, IBEGSchur,
40 & ThresMinINIT
41 LOGICAL :: AGG6_loc
42 INTEGER :: THD_AGG
43 LOGICAL :: SchurON
44 INTEGER :: DEG, DEGME, DEXT, DMAX, E, ELENME, ELN, I,
45 & ILAST, INEXT, J, JLAST, JNEXT, K, KNT1, KNT2, KNT3,
46 & LENJ, LN, ME, MINDEG, NEL,
47 & NLEFT, NVI, NVJ, NVPIV, SLENME, WE, WFLG, WNVI, X
48 INTEGER KNT1_UPDATED, KNT2_UPDATED
49 INTEGER(8) MAXMEM, MEM, NEWMEM
50 INTEGER :: MAXINT_N
51 INTEGER(8) :: HASH, HMOD
52 INTEGER(8) :: P, P1, P2, P3, PDST, PEND, PJ, PME, PME1, PME2,
53 & PN, PSRC, PLN, PELN
54 INTRINSIC max, min, mod
55 IF (n.EQ.1) THEN
56 elen(1) = 1
57 last(1) = 1
58 pe(1) = 0_8
59 nv(1) = 1
60 RETURN
61 ENDIF
62 agg6_loc = agg6
63 thd_agg = max(128, min(n/2048, 1024))
64 IF ( size_complem_list < 0 .OR. size_complem_list > n ) THEN
65 WRITE(*,*) "Internal MUMPS_SYMQAMD_NEW", size_complem_list,n
66 CALL mumps_abort()
67 ENDIF
68 IF (job.EQ.2) THEN
69 schuron = .false.
70 ENDIF
71 thresm = thresh
72 IF (job.NE.2) THEN
73 schuron = (size_complem_list > 0)
74 IF ((job.EQ.1) .AND. (.NOT.schuron) .AND. (n .GT. 0)) THEN
75 ENDIF
76 ibegschur = n-size_complem_list+1
77 IF (thresm.GT.n) thresm = n
78 IF (thresm.LT.0) thresm = 0
79 IF ( schuron ) THEN
80 DO i= 1, n
81 IF ( perm(i) .GE. ibegschur) THEN
82 perm(i) = n + 1
83 IF (len(i) .EQ.0) THEN
84 pe(i) = 0_8
85 ENDIF
86 ENDIF
87 ENDDO
88 ENDIF
89 ENDIF
90 IF (schuron) THEN
91 thresm = n
92 thresmin = n
93 thresprev = n
94 ELSE
95 thresm = max(int(31*n/32),thresm)
96 thresm = max(thresm,1)
97 thresmin = max( 3*thresm / 4, 1)
98 thresprev = thresm
99 ENDIF
100 thresmininit = thresmin/4
101 IF (thresm.GT.0) THEN
102 IF ((thresm.GT.n).OR.(thresm.LT.2)) THEN
103 thresm = n
104 ENDIF
105 ENDIF
106 IF (job.EQ.2) THEN
107 ENDIF
108 permeqn = 0
109 lastd = 0
110 nbd = 0
111 nbed = 0
112 nbdm = 0
113 nel = 0
114 wflg = 2
115 maxint_n=huge(wflg)-n
116 mindeg = 1
117 ncmpa = 0
118 hmod = int(max(1, n-1),kind=8)
119 dmax = 0
120 mem = pfree - 1
121 maxmem = mem
122 DO 10 i = 1, n
123 ndense(i)= 0
124 last(i) = 0
125 head(i) = 0
126 next(i) = 0
127 nv(i) = 1
128 w(i) = 1
129 10 CONTINUE
130 IF (job.EQ.2) THEN
131 DO i = 1,size_complem_list
132 x = complem_list(i)
133 elen(x) = -i
134 nv(x) = len(x)+1
135 dmax = max(dmax, len(x))
136 ENDDO
137 nel = nel + size_complem_list
138 DO i=1,n
139 degree(i) = len(i)
140 ENDDO
141 ELSE
142 DO i=1, n
143 elen(i) = 0
144 degree(i) = len(i)
145 ENDDO
146 ENDIF
147 DO 20 i = 1, n
148 IF (elen(i).LT.0) cycle
149 deg = degree(i)
150 IF (perm(i).EQ.n) THEN
151 permeqn = i
152 perm(i) = n-1
153 ENDIF
154 fdeg = perm(i)
155 IF ( (deg .GT. 0).OR.(perm(i).EQ.n+1) ) THEN
156 IF ( (thresm.GT.0) .AND.
157 & (fdeg .GT.thresm) ) THEN
158 nbd = nbd+1
159 IF (fdeg.NE.n+1) THEN
160 degree(i) = degree(i)+n+2
161 deg = n
162 inext = head(deg)
163 IF (inext .NE. 0) last(inext) = i
164 next(i) = inext
165 head(deg) = i
166 last(i) = 0
167 IF (lastd.EQ.0) lastd=i
168 ELSE
169 nbed = nbed+1
170 degree(i) = n+1
171 deg = n
172 IF (lastd.EQ.0) THEN
173 lastd = i
174 head(deg) = i
175 next(i) = 0
176 last(i) = 0
177 ELSE
178 next(lastd) = i
179 last(i) = lastd
180 lastd = i
181 next(i) = 0
182 ENDIF
183 ENDIF
184 ELSE
185 inext = head(fdeg)
186 IF (inext .NE. 0) last(inext) = i
187 next(i) = inext
188 head(fdeg) = i
189 ENDIF
190 ELSE
191 nel = nel + 1
192 elen(i) = -nel
193 pe(i) = 0_8
194 w(i) = 0
195 ENDIF
196 20 CONTINUE
197 IF ((nbd.EQ.0).AND.(thresm.GT.0)) thresm = n
198 30 IF (nel .LT. n) THEN
199 DO 40 deg = mindeg, n
200 me = head(deg)
201 IF (me .GT. 0) GO TO 50
202 40 CONTINUE
203 50 mindeg = deg
204 IF ( (deg.NE.n) .AND.
205 & (deg.GT.thresm+1) .AND. (nbd.GT.0) ) THEN
206 mindeg = n
207 GOTO 30
208 ENDIF
209 IF (degree(me).LE.n) THEN
210 inext = next(me)
211 IF (inext .NE. 0) last(inext) = 0
212 head(deg) = inext
213 ELSE
214 mindeg = 1
215 nbdm = max(nbdm,nbd)
216 IF (degree(me).GT.n+1) THEN
217 IF (wflg .GT. maxint_n) THEN
218 DO 52 x = 1, n
219 IF (w(x) .NE. 0) w(x) = 1
220 52 CONTINUE
221 wflg = 2
222 ENDIF
223 wflg = wflg + 1
224 51 CONTINUE
225 inext = next(me)
226 IF (inext .NE. 0) THEN
227 last(inext) = 0
228 ELSE
229 lastd = 0
230 ENDIF
231 ndense(me) = 0
232 w(me) = wflg
233 p1 = pe(me)
234 p2 = p1 + int(len(me) -1,8)
235 pln = p1
236 peln = p1
237 DO 55 p=p1,p2
238 e= iw(p)
239 IF (w(e).EQ.wflg) GOTO 55
240 w(e) = wflg
241 IF (pe(e).LT.0_8) THEN
242 x = e
243 53 x = int(-pe(x))
244 IF (w(x) .EQ.wflg) GOTO 55
245 w(x) = wflg
246 IF ( pe(x) .LT. 0 ) GOTO 53
247 e = x
248 ENDIF
249 IF (elen(e).LT.0) THEN
250 ndense(e) = ndense(e) - nv(me)
251 iw(pln) = iw(peln)
252 iw(peln) = e
253 pln = pln + 1_8
254 peln = peln + 1_8
255 pme1 = pe(e)
256 DO 54 pme = pme1, pme1+len(e)-1
257 x = iw(pme)
258 IF ((elen(x).GE.0).AND.(w(x).NE.wflg)) THEN
259 ndense(me) = ndense(me) + nv(x)
260 w(x) = wflg
261 ENDIF
262 54 CONTINUE
263 ELSE
264 ndense(me) = ndense(me) + nv(e)
265 iw(pln)=e
266 pln = pln+1_8
267 ENDIF
268 55 CONTINUE
269 wflg = wflg + 1
270 len(me) = int(pln-p1)
271 elen(me) = int(peln-p1)
272 ndme = ndense(me)+nv(me)
273 IF (ndense(me).EQ.0) ndense(me) =1
274 degree(me) = ndense(me)
275 deg = perm(me)
276 mindeg = min(deg,mindeg)
277 jnext = head(deg)
278 IF (jnext.NE. 0) last(jnext) = me
279 next(me) = jnext
280 head(deg) = me
281 me = inext
282 IF (me.NE.0) THEN
283 IF (degree(me).GT.(n+1) ) GOTO 51
284 ENDIF
285 head(n) = me
286 IF (thresm.LT.n) THEN
287 thresmin = max(thresm+thresmin,thresprev+thresmin/2+1)
288 thresmin = min(thresmin, n)
289 thresprev = thresprev+(n-thresprev)/2+thresmininit
290 thresm = max(
291 & thresm + int(sqrt(dble(thresmin)))+ thresmininit ,
292 & thresprev)
293 thresm = min(thresm,n)
294 thresmin = min(thresm, thresmin)
295 thresprev = thresm
296 ENDIF
297 nbd = nbed
298 GOTO 30
299 ENDIF
300 IF (degree(me).EQ.n+1) THEN
301 IF (nbd.NE.nbed) THEN
302 write(6,*) ' ERROR in MUMPS_SYMQAMD_NEW ',
303 & ' quasi dense rows remains'
304 CALL mumps_abort()
305 ENDIF
306 IF (job.EQ.1) THEN
307 DO i = 1,size_complem_list
308 x = complem_list(i)
309 elen(x) = -(n-size_complem_list+i)
310 nv(x) = 1
311 pe(x) = 0_8
312 ENDDO
313 GOTO 265
314 ENDIF
315 nelme = -(nel+1)
316 DO 59 x=1,n
317 IF ((pe(x).GT.0_8) .AND. (elen(x).LT.0)) THEN
318 pe(x) = int(-complem_list(1),8)
319 ELSEIF (degree(x).EQ.n+1) THEN
320 nel = nel + nv(x)
321 pe(x) = int(-me,8)
322 elen(x) = 0
323 nv(x) = 0
324 ENDIF
325 59 CONTINUE
326 elen(me) = nelme
327 nv(me) = nbd
328 pe(me) = 0_8
329 IF (nel.NE.n) THEN
330 write(6,*) 'Internal ERROR 2 detected in QAMD'
331 write(6,*) ' NEL not equal to N: N, NEL =',n,nel
332 CALL mumps_abort()
333 ENDIF
334 IF (me.NE. complem_list(1)) THEN
335 DO i=1, size_complem_list
336 pe(complem_list(i)) = int(-complem_list(1),8)
337 ENDDO
338 pe(complem_list(1)) = 0_8
339 nv( complem_list(1))= nv(me)
340 nv(me) = 0
341 elen( complem_list(1)) = elen(me)
342 elen(me) = 0
343 ENDIF
344 GOTO 265
345 ENDIF
346 ENDIF
347 elenme = elen(me)
348 elen(me) = - (nel + 1)
349 nvpiv = nv(me)
350 nel = nel + nvpiv
351 ndense(me) = 0
352 nv(me) = -nvpiv
353 degme = 0
354 IF (elenme .EQ. 0) THEN
355 pme1 = pe(me)
356 pme2 = pme1 - 1
357 DO 60 p = pme1, pme1 + int(len(me) - 1,8)
358 i = iw(p)
359 nvi = nv(i)
360 IF (nvi .GT. 0) THEN
361 degme = degme + nvi
362 nv(i) = -nvi
363 pme2 = pme2 + 1
364 iw(pme2) = i
365 IF (degree(i).LE.n) THEN
366 ilast = last(i)
367 inext = next(i)
368 IF (inext .NE. 0) last(inext) = ilast
369 IF (ilast .NE. 0) THEN
370 next(ilast) = inext
371 ELSE
372 head(perm(i)) = inext
373 ENDIF
374 ELSE
375 ndense(me) = ndense(me) + nvi
376 ENDIF
377 ENDIF
378 60 CONTINUE
379 newmem = 0
380 ELSE
381 p = pe(me)
382 pme1 = pfree
383 slenme = len(me) - elenme
384 knt1_updated = 0
385 DO 120 knt1 = 1, elenme + 1
386 knt1_updated = knt1_updated +1
387 IF (knt1 .GT. elenme) THEN
388 e = me
389 pj = p
390 ln = slenme
391 ELSE
392 e = iw(p)
393 p = p + 1
394 pj = pe(e)
395 ln = len(e)
396 ENDIF
397 knt2_updated = 0
398 DO 110 knt2 = 1, ln
399 knt2_updated = knt2_updated+1
400 i = iw(pj)
401 pj = pj + 1
402 nvi = nv(i)
403 IF (nvi .GT. 0) THEN
404 IF (pfree .GT. iwlen) THEN
405 pe(me) = p
406 len(me) = len(me) - knt1_updated
407 knt1_updated = 0
408 IF (len(me) .EQ. 0) pe(me) = 0_8
409 pe(e) = pj
410 len(e) = ln - knt2_updated
411 knt2_updated = 0
412 IF (len(e) .EQ. 0) pe(e) = 0_8
413 ncmpa = ncmpa + 1
414 DO 70 j = 1, n
415 pn = pe(j)
416 IF (pn .GT. 0) THEN
417 pe(j) = int(iw(pn),8)
418 iw(pn) = -j
419 ENDIF
420 70 CONTINUE
421 pdst = 1
422 psrc = 1
423 pend = pme1 - 1
424 80 CONTINUE
425 IF (psrc .LE. pend) THEN
426 j = -iw(psrc)
427 psrc = psrc + 1
428 IF (j .GT. 0) THEN
429 iw(pdst) = int(pe(j))
430 pe(j) = pdst
431 pdst = pdst + 1_8
432 lenj = len(j)
433 DO 90 knt3 = 0, lenj - 2
434 iw(pdst + knt3) = iw(psrc + knt3)
435 90 CONTINUE
436 pdst = pdst + int(lenj - 1,8)
437 psrc = psrc + int(lenj - 1,8)
438 ENDIF
439 GO TO 80
440 ENDIF
441 p1 = pdst
442 DO 100 psrc = pme1, pfree - 1
443 iw(pdst) = iw(psrc)
444 pdst = pdst + 1
445 100 CONTINUE
446 pme1 = p1
447 pfree = pdst
448 pj = pe(e)
449 p = pe(me)
450 ENDIF
451 degme = degme + nvi
452 nv(i) = -nvi
453 iw(pfree) = i
454 pfree = pfree + 1
455 IF (degree(i).LE.n) THEN
456 ilast = last(i)
457 inext = next(i)
458 IF (inext .NE. 0) last(inext) = ilast
459 IF (ilast .NE. 0) THEN
460 next(ilast) = inext
461 ELSE
462 head(perm(i)) = inext
463 ENDIF
464 ELSE
465 ndense(me) = ndense(me) + nvi
466 ENDIF
467 ENDIF
468 110 CONTINUE
469 IF (e .NE. me) THEN
470 pe(e) = int(-me,8)
471 w(e) = 0
472 ENDIF
473 120 CONTINUE
474 pme2 = pfree - 1
475 newmem = pfree - pme1
476 mem = mem + newmem
477 maxmem = max(maxmem, mem)
478 ENDIF
479 degree(me) = degme
480 pe(me) = pme1
481 len(me) = int(pme2 - pme1 + 1_8)
482 IF (wflg .GT. maxint_n) THEN
483 DO 130 x = 1, n
484 IF (w(x) .NE. 0) w(x) = 1
485 130 CONTINUE
486 wflg = 2
487 ENDIF
488 DO 150 pme = pme1, pme2
489 i = iw(pme)
490 IF (degree(i).GT.n) GOTO 150
491 eln = elen(i)
492 IF (eln .GT. 0) THEN
493 nvi = -nv(i)
494 wnvi = wflg - nvi
495 DO 140 p = pe(i), pe(i) + int(eln - 1,8)
496 e = iw(p)
497 we = w(e)
498 IF (we .GE. wflg) THEN
499 we = we - nvi
500 ELSE IF (we .NE. 0) THEN
501 we = degree(e) + wnvi - ndense(e)
502 ENDIF
503 w(e) = we
504 140 CONTINUE
505 ENDIF
506 150 CONTINUE
507 agg6_loc = (agg6 .OR. (degree(me) .LT. thd_agg))
508 DO 180 pme = pme1, pme2
509 i = iw(pme)
510 IF (degree(i).GT.n) GOTO 180
511 p1 = pe(i)
512 p2 = p1 + elen(i) - 1
513 pn = p1
514 hash = 0_8
515 deg = 0
516 DO 160 p = p1, p2
517 e = iw(p)
518 dext = w(e) - wflg
519 IF (dext .GT. 0) THEN
520 deg = deg + dext
521 iw(pn) = e
522 pn = pn + 1
523 hash = hash + int(e,kind=8)
524 ELSE IF (.NOT. agg6_loc .AND. dext .EQ. 0) THEN
525 iw(pn) = e
526 pn = pn + 1
527 hash = hash + int(e,kind=8)
528 ELSE IF (agg6_loc .AND. (dext .EQ. 0) .AND.
529 & ((ndense(me).EQ.nbd).OR.(ndense(e).EQ.0))) THEN
530 pe(e) = int(-me,8)
531 w(e) = 0
532 ELSE IF (agg6_loc .AND. dext.EQ.0) THEN
533 iw(pn) = e
534 pn = pn+1
535 hash = hash + int(e,kind=8)
536 ENDIF
537 160 CONTINUE
538 elen(i) = int(pn - p1 + 1_8)
539 p3 = pn
540 DO 170 p = p2 + 1, p1 + len(i) - 1
541 j = iw(p)
542 nvj = nv(j)
543 IF (nvj .GT. 0) THEN
544 IF (degree(j).LE.n) deg=deg+nvj
545 iw(pn) = j
546 pn = pn + 1
547 hash = hash + int(j,kind=8)
548 ENDIF
549 170 CONTINUE
550 IF (((elen(i).EQ.1).AND.(p3.EQ.pn))
551 & .OR.
552 & (agg6_loc.AND.(deg .EQ. 0).AND.(ndense(me).EQ.nbd))
553 & )
554 & THEN
555 pe(i) = int(-me, 8)
556 nvi = -nv(i)
557 degme = degme - nvi
558 nvpiv = nvpiv + nvi
559 nel = nel + nvi
560 nv(i) = 0
561 elen(i) = 0
562 ELSE
563 degree(i) = min(deg+nbd-ndense(me),
564 & degree(i))
565 iw(pn) = iw(p3)
566 iw(p3) = iw(p1)
567 iw(p1) = me
568 len(i) = int(pn - p1 + 1)
569 hash = mod(hash, hmod) + 1_8
570 j = head(hash)
571 IF (j .LE. 0) THEN
572 next(i) = -j
573 head(hash) = -i
574 ELSE
575 next(i) = last(j)
576 last(j) = i
577 ENDIF
578 last(i) = int(hash,kind=kind(last))
579 ENDIF
580 180 CONTINUE
581 degree(me) = degme
582 dmax = max(dmax, degme)
583 wflg = wflg + dmax
584 IF (wflg .GT. maxint_n) THEN
585 DO 190 x = 1, n
586 IF (w(x) .NE. 0) w(x) = 1
587 190 CONTINUE
588 wflg = 2
589 ENDIF
590 DO 250 pme = pme1, pme2
591 i = iw(pme)
592 IF ( (nv(i).LT.0) .AND. (degree(i).LE.n) ) THEN
593 hash = int(last(i),kind=8)
594 j = head(hash)
595 IF (j .EQ. 0) GO TO 250
596 IF (j .LT. 0) THEN
597 i = -j
598 head(hash) = 0
599 ELSE
600 i = last(j)
601 last(j) = 0
602 ENDIF
603 IF (i .EQ. 0) GO TO 250
604 200 CONTINUE
605 IF (next(i) .NE. 0) THEN
606 x = i
607 ln = len(i)
608 eln = elen(i)
609 DO 210 p = pe(i) + 1, pe(i) + int(ln - 1,8)
610 w(iw(p)) = wflg
611 210 CONTINUE
612 jlast = i
613 j = next(i)
614 220 CONTINUE
615 IF (j .NE. 0) THEN
616 IF (len(j) .NE. ln) GO TO 240
617 IF (elen(j) .NE. eln) GO TO 240
618 DO 230 p = pe(j) + 1, pe(j) + int(ln - 1,8)
619 IF (w(iw(p)) .NE. wflg) GO TO 240
620 230 CONTINUE
621 IF (perm(j).GT.perm(x)) THEN
622 pe(j) = int(-x,8)
623 nv(x) = nv(x) + nv(j)
624 nv(j) = 0
625 elen(j) = 0
626 ELSE
627 pe(x) = int(-j,8)
628 nv(j) = nv(x) + nv(j)
629 nv(x) = 0
630 elen(x) = 0
631 x = j
632 ENDIF
633 j = next(j)
634 next(jlast) = j
635 GO TO 220
636 240 CONTINUE
637 jlast = j
638 j = next(j)
639 GO TO 220
640 ENDIF
641 wflg = wflg + 1
642 i = next(i)
643 IF (i .NE. 0) GO TO 200
644 ENDIF
645 ENDIF
646 250 CONTINUE
647 IF ( (thresm .GT. 0).AND.(thresm.LT.n) ) THEN
648 thresm = max(thresmin, thresm-nvpiv)
649 ENDIF
650 p = pme1
651 nleft = n - nel
652 DO 260 pme = pme1, pme2
653 i = iw(pme)
654 nvi = -nv(i)
655 IF (nvi .GT. 0) THEN
656 nv(i) = nvi
657 IF (degree(i).LE.n) THEN
658 deg = min(degree(i)+ degme - nvi, nleft - nvi)
659 degree(i) = deg
660 idense = .false.
661 IF (thresm.GT.0) THEN
662 IF (perm(i) .GT. thresm) THEN
663 idense = .true.
664 degree(i) = degree(i)+n+2
665 ENDIF
666 IF (idense) THEN
667 p1 = pe(i)
668 p2 = p1 + int(elen(i) - 1, 8)
669 IF (p2.GE.p1) THEN
670 DO 264 pj=p1,p2
671 e= iw(pj)
672 ndense(e) = ndense(e) + nvi
673 264 CONTINUE
674 ENDIF
675 nbd = nbd+nvi
676 fdeg = n
677 deg = n
678 inext = head(deg)
679 IF (inext .NE. 0) last(inext) = i
680 next(i) = inext
681 head(deg) = i
682 last(i) = 0
683 IF (lastd.EQ.0) lastd=i
684 ENDIF
685 ENDIF
686 IF (.NOT.idense) THEN
687 fdeg = perm(i)
688 inext = head(fdeg)
689 IF (inext .NE. 0) last(inext) = i
690 next(i) = inext
691 last(i) = 0
692 head(fdeg) = i
693 ENDIF
694 mindeg = min(mindeg, fdeg)
695 ENDIF
696 iw(p) = i
697 p = p + 1
698 ENDIF
699 260 CONTINUE
700 nv(me) = nvpiv + degme
701 len(me) = int(p - pme1)
702 IF (len(me) .EQ. 0) THEN
703 pe(me) = 0_8
704 w(me) = 0
705 ENDIF
706 IF (newmem .NE. 0) THEN
707 pfree = p
708 mem = mem - newmem + int(len(me),8)
709 ENDIF
710 GO TO 30
711 ENDIF
712 265 CONTINUE
713 DO 290 i = 1, n
714 IF (elen(i) .EQ. 0) THEN
715 j = int(-pe(i))
716 270 CONTINUE
717 IF (elen(j) .GE. 0) THEN
718 j = int(-pe(j))
719 GO TO 270
720 ENDIF
721 e = j
722 k = -elen(e)
723 j = i
724 280 CONTINUE
725 IF (elen(j) .GE. 0) THEN
726 jnext = int(-pe(j))
727 pe(j)= int(-e,8)
728 IF (elen(j) .EQ. 0) THEN
729 elen(j) = k
730 k = k + 1
731 ENDIF
732 j = jnext
733 GO TO 280
734 ENDIF
735 elen(e) = -k
736 ENDIF
737 290 CONTINUE
738 DO 300 i = 1, n
739 k = abs(elen(i))
740 last(k) = i
741 elen(i) = k
742 300 CONTINUE
743 IF (.NOT.schuron) THEN
744 IF (permeqn.GT.0) perm(permeqn) = n
745 ENDIF
746 pfree = maxmem
747 RETURN
748 END SUBROUTINE mumps_symqamd_new
749 SUBROUTINE mumps_wrap_ginp94
750 & ( n, ipe, iw, liw8,
751 & perm, sizeofblocks,
752 & keep60, listvar_schur, size_schur, keep378,
753 & colcount, parent,
754 & porder, iwtmp1, iwtmp2, iwtmp3, iwtmp4,
755 & iwtmp5,
756 & info )
757 IMPLICIT NONE
758 INTEGER, INTENT(IN) :: N
759 INTEGER(8), INTENT(IN) :: LIW8
760 INTEGER(8), INTENT(IN) :: IPE(N+1)
761 INTEGER, INTENT(IN) :: SizeOfBlocks(N)
762 INTEGER, INTENT(INOUT) :: PERM(N)
763 INTEGER, INTENT(IN) :: IW(LIW8)
764 INTEGER, INTENT(OUT) :: COLCOUNT(N)
765 INTEGER, INTENT(OUT) :: PARENT(N)
766 INTEGER, INTENT(IN) :: KEEP60, SIZE_SCHUR
767 INTEGER, INTENT(IN) :: LISTVAR_SCHUR(SIZE_SCHUR)
768 INTEGER, INTENT(IN) :: KEEP378
769 INTEGER, INTENT(INOUT) :: INFO(2)
770 INTEGER, INTENT(OUT):: PORDER(N), IWTMP1(N), IWTMP2(N)
771 INTEGER, INTENT(OUT):: IWTMP3(N), IWTMP4(N), IWTMP5(N)
772 INTEGER :: I, KEEP378_loc
773 LOGICAL :: SizeOfBlocks_provided
774 sizeofblocks_provided = (sizeofblocks(1).NE.-1)
775 IF (keep378.NE.0) keep378_loc=1
776 DO i=1, n
777 iwtmp1(perm(i)) = i
778 END DO
779 CALL mumps_ginp94_elim_tree (
780 & n, ipe, iw, liw8, iwtmp1, perm, parent,
781 & iwtmp2, info)
782 IF (info(1).LT.0) RETURN
783 CALL mumps_ginp94_postorder(parent, n, porder,
784 & iwtmp1, iwtmp2, iwtmp3,
785 & info)
786 IF (info(1).LT.0) RETURN
787 CALL mumps_ginp94_colcounts(
788 & n, liw8, ipe, iw, parent, porder, colcount,
789 & sizeofblocks_provided, sizeofblocks, keep378_loc,
790 & iwtmp1, iwtmp2, iwtmp3, iwtmp4, iwtmp5,
791 & info)
792 IF (info(1).LT.0) RETURN
793 IF (keep60.NE.0) THEN
794 CALL mumps_ginp94_postprocess_schur (
795 & n, parent, colcount, perm,
796 & listvar_schur, size_schur )
797 ENDIF
798 RETURN
799 END SUBROUTINE mumps_wrap_ginp94
800 SUBROUTINE mumps_ginp94_elim_tree (
801 & n, iptr, jcn, ljcn, invperm, perm, parent,
802 & ancestor, info)
803 IMPLICIT NONE
804 INTEGER(8), INTENT(IN) :: ljcn
805 integer :: n
806 INTEGER(8), INTENT(IN) :: iptr(n+1)
807 integer, INTENT(IN) :: jcn(ljcn), invperm(n), perm(n)
808 integer, INTENT(OUT) :: parent(n)
809 integer, INTENT(INOUT) :: INFO(2)
810 integer, INTENT(OUT) :: ancestor(n)
811 integer :: jpos, i, j, k
812 integer(8) :: iidx8
813 ancestor=0
814 parent =0
815 do jpos = 1, n
816 j = invperm(jpos)
817 do iidx8 = iptr(j), iptr(j+1)-1
818 i = jcn(iidx8)
819 if (perm(i).ge.jpos) cycle
820 k = i
821 call add_node(n, ancestor, parent, j, k)
822 end do
823 end do
824 return
825 contains
826 subroutine add_node(n, ancestor, parent, j, i)
827 implicit none
828 integer, intent(in):: n
829 integer :: parent(n)
830 integer :: ancestor(n)
831 integer :: i, j
832 integer :: r, k
833 if(i.eq.0) return
834 k = i
835 do
836 r = ancestor(k)
837 if (r .eq. j) then
838 return
839 end if
840 ancestor(k) = j
841 if(r .eq. 0) then
842 parent(k) = j
843 return
844 end if
845 k = r
846 end do
847 end subroutine add_node
848 END SUBROUTINE mumps_ginp94_elim_tree
849 SUBROUTINE mumps_ginp94_postorder(parent, n, porder,
850 & son, brother, stack,
851 & INFO
852 & )
853 IMPLICIT NONE
854 integer, intent(in) :: n
855 integer, intent(in) :: parent(n)
856 integer, intent(out) :: porder(n)
857 integer, intent(inout):: INFO(2)
858 integer, intent(out) :: son(n), brother(n), stack(n)
859 integer :: i, father, br, head, hp, pp, t
860 son = 0
861 do i=n, 1, -1
862 father = parent(i)
863 if (father .ne. 0) then
864 br = son(father)
865 brother(i) = br
866 son(father) = i
867 end if
868 end do
869 head = 0
870 hp = 0
871 pp = 1
872 do t=1, n
873 if (parent(t) .ne. 0) cycle
874 hp = hp+1
875 stack(hp) = t
876 head = t
877 do
878 if(son(head) .eq. 0) then
879 porder(pp) = head
880 pp = pp+1
881 hp = hp-1
882 if (parent(head) .ne. 0) then
883 son(parent(head)) = brother(head)
884 end if
885 if (hp .eq. 0) then
886 exit
887 end if
888 head = stack(hp)
889 else
890 hp = hp+1
891 stack(hp) = son(head)
892 head = son(head)
893 end if
894 end do
895 end do
896 RETURN
897 END SUBROUTINE mumps_ginp94_postorder
898 SUBROUTINE mumps_ginp94_colcounts(
899 & n, ljcn, iptr, jcn, parent, porder, cc,
900 & SizeOfBlocks_provided, SizeOfBlocks, KEEP378,
901 & fst_desc, iporder, prev_p, prev_nbr, setpath,
902 & INFO)
903 IMPLICIT NONE
904 integer, intent(in) :: n
905 integer(8), intent(in) :: ljcn
906 integer(8), intent(in) :: iptr(n+1)
907 integer, intent(in) :: jcn(ljcn)
908 integer, intent(inout) :: parent(n), porder(n)
909 integer, intent(in) :: SizeOfBlocks(n)
910 logical, intent(in) :: SizeOfBlocks_provided
911 integer, intent(in) :: KEEP378
912 integer, intent(out) :: cc(n)
913 integer, intent(inout):: INFO(2)
914 integer, intent(out) :: fst_desc(n), iporder(n), prev_p(n)
915 integer, intent(out) :: prev_nbr(n), setpath(n)
916 integer :: i, curr, fd, j, jidx, k
917 integer(8) :: iidx8
918 integer :: f, ref, p_leaf, q, jj
919 do j=1, n
920 iporder(porder(j)) = j
921 end do
922 cc = 0
923 fst_desc=-1
924 do i=1, n
925 curr = porder(i)
926 fd = curr
927 if(fst_desc(curr) .eq. -1) then
928 if (sizeofblocks_provided) then
929 cc(curr) = sizeofblocks(curr)
930 else
931 cc(curr) = 1
932 endif
933 end if
934 do
935 if (fst_desc(curr) .gt. 0) exit
936 fst_desc(curr) = fd
937 if (parent(curr) .eq. 0) exit
938 curr = parent(curr)
939 end do
940 end do
941 do j=1, n
942 setpath(j)=j
943 end do
944 prev_p = 0
945 prev_nbr = 0
946 do jidx=1, n
947 j = abs(porder(jidx))
948#if defined(check)
949 if (parent(j).eq.0)
950 & write(6,*) " ========= jidx,j= ", jidx,j," is a rootnode "
951#endif
952 if(parent(j) .ne. 0) then
953 if (keep378.eq.1) then
954 if (cc(parent(j)) .lt. 0) then
955 porder(iporder(parent(j)))= -parent(j)
956 endif
957 endif
958 if (sizeofblocks_provided) then
959 cc(parent(j)) = cc(parent(j)) - sizeofblocks(j)
960 else
961 cc(parent(j)) = cc(parent(j))-1
962 endif
963 endif
964 do iidx8=iptr(j), iptr(j+1)-1
965 i = jcn(iidx8)
966 if (iporder(i).le.jidx) cycle
967 if(prev_nbr(i) .eq. 0) then
968 ref = 0
969 else
970 ref = iporder(prev_nbr(i))
971 end if
972 if(iporder(fst_desc(j)) .gt. ref) then
973 if (keep378.eq.1) then
974 porder(iporder(j))= -j
975 endif
976 if (sizeofblocks_provided) then
977 cc(j) = cc(j) + sizeofblocks(i)
978 else
979 cc(j) = cc(j) + 1
980 endif
981 p_leaf = prev_p(i)
982 if (p_leaf .ne. 0) then
983 q = setfind(setpath, p_leaf)
984 if (sizeofblocks_provided) then
985 cc(q) = cc(q) - sizeofblocks(i)
986 else
987 cc(q) = cc(q) - 1
988 endif
989 end if
990 prev_p(i) = j
991 end if
992 prev_nbr(i) = j
993 end do
994 if (parent(j).ne.0) setpath(j)=parent(j)
995 end do
996 do jj=1, n-1
997 j=abs(porder(jj))
998 if(parent(j) .ne. 0) cc(parent(j)) = cc(parent(j)) + cc(j)
999 end do
1000 if (keep378.eq.1) then
1001 i=1
1002 do while (i .lt. n)
1003 porder(i) = abs(porder(i))
1004 j = i+1
1005 if (parent(porder(i)).ne.0) then
1006 do while (porder(j) .gt. 0)
1007 j = j+1
1008 if (parent(abs(porder(j-1))).eq.0) exit
1009 if (j.gt.n) exit
1010 end do
1011 endif
1012 parent(porder(i)) = parent(porder(j-1))
1013 do k=i+1, j-1
1014 parent(porder(k)) = -porder(i)
1015 cc(porder(k)) = 0
1016 end do
1017 i = j
1018 end do
1019 porder(n) = abs(porder(n))
1020 do i=1,n-1
1021 f = abs(parent(i))
1022 if (f.eq.0) cycle
1023 if (cc(f).eq.0) then
1024 parent(i) = parent(f)
1025 endif
1026 end do
1027 endif
1028 do i=1,n
1029 f = parent(i)
1030 if (parent(i).gt.0) then
1031 parent(i) = -parent(i)
1032 endif
1033 end do
1034 return
1035 contains
1036 function setfind(setpath, p_leaf)
1037 implicit none
1038 integer :: setpath(:), p_leaf, setfind
1039 integer :: q, c, tmp
1040 q=p_leaf
1041 do while (setpath(q) .ne.q)
1042 q = setpath(q)
1043 end do
1044 c = p_leaf
1045 do while (c .ne.q)
1046 tmp = setpath(c)
1047 setpath(c) = q
1048 c = tmp
1049 end do
1050 setfind = q
1051 return
1052 end function setfind
1053 END SUBROUTINE mumps_ginp94_colcounts
1054 SUBROUTINE mumps_ginp94_postprocess_schur (
1055 & N, PARENT, COLCOUNT, PERM,
1056 & LISTVAR_SCHUR, SIZE_SCHUR
1057 & )
1058 IMPLICIT NONE
1059 INTEGER, INTENT(IN) :: N, SIZE_SCHUR
1060 INTEGER, INTENT(IN) :: PERM(N), LISTVAR_SCHUR(SIZE_SCHUR)
1061 INTEGER, INTENT(INOUT) :: PARENT(N), COLCOUNT(N)
1062 INTEGER I, FIRSTinSchur, PrincipalVarSchur
1063 firstinschur = n-size_schur+1
1064 principalvarschur = listvar_schur(1)
1065 DO i=1, n
1066 IF (i.EQ.principalvarschur) THEN
1067 parent(i) = 0
1068 colcount(i) = size_schur
1069 ELSE IF (parent(i).EQ.0.AND.perm(i).LT.firstinschur) THEN
1070 parent(i) = -principalvarschur
1071 ELSE IF (perm(i).GE.firstinschur) THEN
1072 parent(i) = -principalvarschur
1073 colcount(i) = 0
1074 ELSE IF (perm(-parent(i)).GE.firstinschur) THEN
1075 parent(i) = -principalvarschur
1076 ENDIF
1077 ENDDO
1078 RETURN
1079 END SUBROUTINE mumps_ginp94_postprocess_schur
mumps_abort
#define mumps_abort
Definition VE_Metis.h:25
add_node
subroutine add_node(n, ancestor, parent, j, i)
Definition ana_AMDMF.F:827
mumps_ginp94_postorder
subroutine mumps_ginp94_postorder(parent, n, porder, son, brother, stack, info)
Definition ana_AMDMF.F:853
mumps_ginp94_elim_tree
subroutine mumps_ginp94_elim_tree(n, iptr, jcn, ljcn, invperm, perm, parent, ancestor, info)
Definition ana_AMDMF.F:803
mumps_ginp94_postprocess_schur
subroutine mumps_ginp94_postprocess_schur(n, parent, colcount, perm, listvar_schur, size_schur)
Definition ana_AMDMF.F:1058
mumps_ginp94_colcounts
subroutine mumps_ginp94_colcounts(n, ljcn, iptr, jcn, parent, porder, cc, sizeofblocks_provided, sizeofblocks, keep378, fst_desc, iporder, prev_p, prev_nbr, setpath, info)
Definition ana_AMDMF.F:903
mumps_wrap_ginp94
subroutine mumps_wrap_ginp94(n, ipe, iw, liw8, perm, sizeofblocks, keep60, listvar_schur, size_schur, keep378, colcount, parent, porder, iwtmp1, iwtmp2, iwtmp3, iwtmp4, iwtmp5, info)
Definition ana_AMDMF.F:757
mumps_symqamd_new
subroutine mumps_symqamd_new(job, thresh, ndense, n, iwlen, pe, pfree, len, iw, nv, elen, last, ncmpa, degree, head, next, w, perm, complem_list, size_complem_list, agg6)
Definition ana_AMDMF.F:20
setfind
integer function setfind(setpath, p_leaf)
Definition ana_AMDMF.F:1037
min
#define min(a, b)
Definition macros.h:20