OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cana_LDLT_preprocess.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 cmumps_set_constraints(
15 & N,PIV,FRERE,FILS,NFSIZ,IKEEP,
16 & NCST,KEEP,KEEP8, ROWSCA
17 & )
18 USE cmumps_struc_def
19 IMPLICIT NONE
20 INTEGER, INTENT(IN) :: N
21 INTEGER, INTENT(OUT) :: NCST
22 INTEGER :: PIV(N),FRERE(N),FILS(N),NFSIZ(N),IKEEP(N)
23 INTEGER :: KEEP(500)
24 INTEGER(8) :: KEEP8(150)
25 REAL :: ROWSCA(N)
26 INTEGER I,P11,P1,P2,K1,K2,NLOCKED
27 LOGICAL V1,V2
28 ncst = 0
29 nlocked = 0
30 p11 = keep(93)
31 DO i=keep(93)-1,1,-2
32 p1 = piv(i)
33 p2 = piv(i+1)
34 k1 = ikeep(p1)
35 IF (k1 .NE. 0) THEN
36 v1 = (k1+2*exponent(rowsca(p1)) .GE. -3)
37 ELSE
38 v1 = .false.
39 ENDIF
40 k2 = ikeep(p2)
41 IF (k2 .NE. 0) THEN
42 v2 = (k2+exponent(rowsca(p2)**2) .GE. -3)
43 ELSE
44 v2 = .false.
45 ENDIF
46 IF(v1 .AND. v2) THEN
47 piv(p11) = p1
48 p11 = p11 - 1
49 piv(p11) = p2
50 p11 = p11 - 1
51 ELSE IF(v1) THEN
52 ncst = ncst+1
53 frere(ncst) = p1
54 ncst = ncst+1
55 frere(ncst) = p2
56 ELSE IF(v2) THEN
57 ncst = ncst+1
58 frere(ncst) = p2
59 ncst = ncst+1
60 frere(ncst) = p1
61 ELSE
62 nlocked = nlocked + 1
63 fils(nlocked) = p1
64 nlocked = nlocked + 1
65 fils(nlocked) = p2
66 ENDIF
67 ENDDO
68 DO i=1,nlocked
69 piv(i) = fils(i)
70 ENDDO
71 keep(94) = keep(94) + keep(93) - nlocked
72 keep(93) = nlocked
73 DO i=1,ncst
74 nlocked = nlocked + 1
75 piv(nlocked) = frere(i)
76 ENDDO
77 DO i=1,keep(93)/2
78 nfsiz(i) = 0
79 ENDDO
80 DO i=(keep(93)/2)+1,(keep(93)/2)+ncst,2
81 nfsiz(i) = i+1
82 nfsiz(i+1) = -1
83 ENDDO
84 DO i=(keep(93)/2)+ncst+1,(keep(93)/2)+keep(94)
85 nfsiz(i) = 0
86 ENDDO
87 END SUBROUTINE cmumps_set_constraints
88 SUBROUTINE cmumps_expand_permutation(N,NCMP,N11,N22,PIV,
89 & INVPERM,PERM)
90 IMPLICIT NONE
91 INTEGER N11,N22,N,NCMP
92 INTEGER, intent(in) :: PIV(N),PERM(N)
93 INTEGER, intent (out):: INVPERM(N)
94 INTEGER CMP_POS,EXP_POS,I,J,N2,K
95 n2 = n22/2
96 exp_pos = 1
97 DO cmp_pos=1,ncmp
98 j = perm(cmp_pos)
99 IF(j .LE. n2) THEN
100 k = 2*j-1
101 i = piv(k)
102 invperm(i) = exp_pos
103 exp_pos = exp_pos+1
104 k = k+1
105 i = piv(k)
106 invperm(i) = exp_pos
107 exp_pos = exp_pos+1
108 ELSE
109 k = n2 + j
110 i = piv(k)
111 invperm(i) = exp_pos
112 exp_pos = exp_pos+1
113 ENDIF
114 ENDDO
115 DO k=n22+n11+1,n
116 i = piv(k)
117 invperm(i) = exp_pos
118 exp_pos = exp_pos+1
119 ENDDO
120 RETURN
121 END SUBROUTINE cmumps_expand_permutation
122 SUBROUTINE cmumps_ldlt_compress(
123 & N,NZ, IRN, ICN, PIV,
124 & NCMP, IW, LW, IPE, LEN, IQ,
125 & FLAG, ICMP, IWFR,
126 & IERROR, KEEP,KEEP8, ICNTL,INPLACE64_GRAPH_COPY)
127 IMPLICIT NONE
128 INTEGER, intent(in) :: N
129 INTEGER(8), intent(in) :: NZ, LW
130 INTEGER, intent(in) :: IRN(NZ), ICN(NZ), PIV(N)
131 INTEGER, intent(in) :: ICNTL(60)
132 INTEGER, intent(in) :: KEEP(500)
133 INTEGER(8), intent(in) :: KEEP8(150)
134 INTEGER, intent(out) :: NCMP, IERROR
135 INTEGER(8), intent(out) :: IWFR, IPE(N+1)
136 INTEGER, intent(out) :: IW(LW)
137 INTEGER, intent(out) :: LEN(N)
138 INTEGER(8), intent(out) :: IQ(N)
139 INTEGER, intent(out) :: FLAG(N), ICMP(N)
140 LOGICAL, intent(inout) :: INPLACE64_GRAPH_COPY
141 INTEGER :: MP, N11, N22
142 INTEGER :: I, J, N1, K
143 INTEGER(8) :: NDUP, L, K8, K1, K2, LAST
144 mp = icntl(2)
145 ierror = 0
146 n22 = keep(93)
147 n11 = keep(94)
148 ncmp = n22/2 + n11
149 DO i=1,ncmp
150 ipe(i) = 0
151 ENDDO
152 k = 1
153 DO i=1,n22/2
154 j = piv(k)
155 icmp(j) = i
156 k = k + 1
157 j = piv(k)
158 icmp(j) = i
159 k = k + 1
160 ENDDO
161 k = n22/2 + 1
162 DO i=n22+1,n22+n11
163 j = piv(i)
164 icmp(j) = k
165 k = k + 1
166 ENDDO
167 DO i=n11+n22+1,n
168 j = piv(i)
169 icmp(j) = 0
170 ENDDO
171 DO k8=1,nz
172 i = irn(k8)
173 j = icn(k8)
174 IF ((i.GT.n).OR.(j.GT.n).OR.(i.LT.1)
175 & .OR.(j.LT.1)) THEN
176 ierror = ierror + 1
177 ELSE
178 i = icmp(i)
179 j = icmp(j)
180 IF ((i.NE.0).AND.(j.NE.0).AND.(i.NE.j)) THEN
181 ipe(i) = ipe(i) + 1_8
182 ipe(j) = ipe(j) + 1_8
183 ENDIF
184 ENDIF
185 ENDDO
186 iq(1) = 1_8
187 n1 = ncmp - 1
188 IF (n1.GT.0) THEN
189 DO i=1,n1
190 iq(i+1) = ipe(i) + iq(i)
191 ENDDO
192 ENDIF
193 last = max(ipe(ncmp)+iq(ncmp)-1_8,iq(ncmp))
194 DO i = 1,ncmp
195 flag(i) = 0
196 ipe(i) = iq(i)
197 ENDDO
198 iw(1:last) = 0
199 iwfr = last + 1_8
200 DO k8=1,nz
201 i = irn(k8)
202 j = icn(k8)
203 IF ((i.GT.n).OR.(j.GT.n).OR.(i.LT.1)
204 & .OR.(j.LT.1)) cycle
205 i = icmp(i)
206 j = icmp(j)
207 IF (i.NE.j) THEN
208 IF (i.LT.j) THEN
209 IF ((i.GE.1).AND.(j.LE.n)) THEN
210 iw(iq(i)) = -j
211 iq(i) = iq(i) + 1_8
212 ENDIF
213 ELSE
214 IF ((j.GE.1).AND.(i.LE.n)) THEN
215 iw(iq(j)) = -i
216 iq(j) = iq(j) + 1_8
217 ENDIF
218 ENDIF
219 ENDIF
220 ENDDO
221 ndup = 0_8
222 DO i=1,ncmp
223 k1 = ipe(i)
224 k2 = iq(i) -1_8
225 IF (k1.GT.k2) THEN
226 len(i) = 0
227 ELSE
228 DO k8=k1,k2
229 j = -iw(k8)
230 IF (j.LE.0) GO TO 250
231 l = iq(j)
232 iq(j) = l + 1_8
233 IF (flag(j).EQ.i) THEN
234 ndup = ndup + 1_8
235 iw(l) = 0
236 iw(k8) = 0
237 ELSE
238 iw(l) = i
239 iw(k8) = j
240 flag(j) = i
241 ENDIF
242 ENDDO
243 250 len(i) = int(iq(i) - ipe(i))
244 ENDIF
245 ENDDO
246 IF (ndup.NE.0_8) THEN
247 iwfr = 1_8
248 DO i=1,ncmp
249 k1 = ipe(i)
250 IF (len(i).EQ.0) THEN
251 ipe(i) = iwfr
252 cycle
253 ENDIF
254 k2 = k1 + len(i) - 1
255 l = iwfr
256 ipe(i) = iwfr
257 DO k8=k1,k2
258 IF (iw(k8).NE.0) THEN
259 iw(iwfr) = iw(k8)
260 iwfr = iwfr + 1_8
261 ENDIF
262 ENDDO
263 len(i) = int(iwfr - l)
264 ENDDO
265 ENDIF
266 ipe(ncmp+1) = ipe(ncmp) + int(len(ncmp),8)
267 iwfr = ipe(ncmp+1)
268 inplace64_graph_copy = (lw.GE.2*(iwfr-1_8))
269 RETURN
270 END SUBROUTINE cmumps_ldlt_compress
271 SUBROUTINE cmumps_sym_mwm(
272 & N, NE, IP, IRN, SCALING,LSC,CPERM, DIAG,
273 & ICNTL, WEIGHT,MARKED,FLAG,
274 & PIV_OUT, INFO)
275 IMPLICIT NONE
276 INTEGER, INTENT(IN) :: N
277 INTEGER(8), INTENT(IN) :: NE
278 INTEGER :: ICNTL(10), INFO(10),LSC
279 INTEGER :: CPERM(N),PIV_OUT(N), IRN(NE), DIAG(N)
280 INTEGER(8), INTENT(IN) :: IP(N+1)
281 REAL :: SCALING(LSC),WEIGHT(N+2)
282 INTEGER :: MARKED(N),FLAG(N)
283 INTEGER :: NUM1,NUM2,NUMTOT,PATH_LENGTH,NLAST
284 INTEGER :: I,BEST_BEG, CUR_EL,CUR_EL_PATH,CUR_EL_PATH_NEXT
285 INTEGER :: L1,L2,TUP,T22
286 INTEGER(8) :: PTR_SET1,PTR_SET2
287 REAL :: BEST_SCORE,CUR_VAL,TMP,VAL
288 REAL INITSCORE, CMUMPS_UPDATESCORE,
289 & cmumps_update_inverse, cmumps_metric2x2
290 LOGICAL VRAI,FAUX,MAX_CARD_DIAG,USE_SCALING
291 INTEGER SUM
292 REAL ZERO,ONE
293 parameter(sum = 1, vrai = .true., faux = .false.)
294 parameter(zero = 0.0e0, one = 1.0e0)
295 max_card_diag = .true.
296 num1 = 0
297 num2 = 0
298 numtot = 0
299 nlast = n
300 info = 0
301 marked = 1
302 flag = 0
303 val = one
304 IF(lsc .GT. 1) THEN
305 use_scaling = .true.
306 ELSE
307 use_scaling = .false.
308 ENDIF
309 tup = icntl(2)
310 IF(tup .EQ. sum) THEN
311 initscore = zero
312 ELSE
313 initscore = one
314 ENDIF
315 IF(icntl(2) .GT. 2 .OR. icntl(2) .LE. 0) THEN
316 WRITE(*,*)
317 & 'ERROR: WRONG VALUE FOR ICNTL(2) = ',icntl(2)
318 info(1) = -1
319 RETURN
320 ENDIF
321 t22 = icntl(1)
322 IF(icntl(1) .LT. 0 .OR. icntl(1) .GT. 2) THEN
323 WRITE(*,*)
324 & 'ERROR: WRONG VALUE FOR ICNTL(1) = ',icntl(1)
325 info(1) = -1
326 RETURN
327 ENDIF
328 DO cur_el=1,n
329 IF(marked(cur_el) .LE. 0) THEN
330 cycle
331 ENDIF
332 IF(cperm(cur_el) .LT. 0) THEN
333 marked(cur_el) = -1
334 cycle
335 ENDIF
336 path_length = 2
337 cur_el_path = cperm(cur_el)
338 IF(cur_el_path .EQ. cur_el) THEN
339 marked(cur_el) = -1
340 cycle
341 ENDIF
342 marked(cur_el) = 0
343 weight(1) = initscore
344 weight(2) = initscore
345 l1 = int(ip(cur_el+1)-ip(cur_el))
346 l2 = int(ip(cur_el_path+1)-ip(cur_el_path))
347 ptr_set1 = ip(cur_el)
348 ptr_set2 = ip(cur_el_path)
349 IF(use_scaling) THEN
350 val = -scaling(cur_el_path) - scaling(cur_el+n)
351 ENDIF
352 cur_val = cmumps_metric2x2(
353 & cur_el,cur_el_path,
354 & irn(ptr_set1),irn(ptr_set2),
355 & l1,l2,
356 & val,diag,n,flag,faux,t22)
357 weight(path_length+1) =
358 & cmumps_updatescore(weight(1),cur_val,tup)
359 DO
360 IF(cur_el_path .EQ. cur_el) EXIT
361 path_length = path_length+1
362 marked(cur_el_path) = 0
363 cur_el_path_next = cperm(cur_el_path)
364 l1 = int(ip(cur_el_path+1)-ip(cur_el_path))
365 l2 = int(ip(cur_el_path_next+1)-ip(cur_el_path_next))
366 ptr_set1 = ip(cur_el_path)
367 ptr_set2 = ip(cur_el_path_next)
368 IF(use_scaling) THEN
369 val = -scaling(cur_el_path_next)
370 & - scaling(cur_el_path+n)
371 ENDIF
372 cur_val = cmumps_metric2x2(
373 & cur_el_path,cur_el_path_next,
374 & irn(ptr_set1),irn(ptr_set2),
375 & l1,l2,
376 & val,diag,n,flag,vrai,t22)
377 weight(path_length+1) =
378 & cmumps_updatescore(weight(path_length-1),cur_val,tup)
379 cur_el_path = cur_el_path_next
380 ENDDO
381 IF(mod(path_length,2) .EQ. 1) THEN
382 IF(weight(path_length+1) .GE. weight(path_length)) THEN
383 cur_el_path = cperm(cur_el)
384 ELSE
385 cur_el_path = cur_el
386 ENDIF
387 DO i=1,(path_length-1)/2
388 num2 = num2+1
389 piv_out(num2) = cur_el_path
390 cur_el_path = cperm(cur_el_path)
391 num2 = num2+1
392 piv_out(num2) = cur_el_path
393 cur_el_path = cperm(cur_el_path)
394 ENDDO
395 numtot = numtot + path_length - 1
396 ELSE
397 IF(max_card_diag) THEN
398 cur_el_path = cperm(cur_el)
399 IF(diag(cur_el) .NE. 0) THEN
400 best_beg = cur_el_path
401 GOTO 1000
402 ENDIF
403 DO i=1,(path_length/2)
404 cur_el_path_next = cperm(cur_el_path)
405 IF(diag(cur_el_path) .NE. 0) THEN
406 best_beg = cur_el_path_next
407 GOTO 1000
408 ENDIF
409 ENDDO
410 ENDIF
411 best_beg = cur_el
412 best_score = weight(path_length-1)
413 cur_el_path = cperm(cur_el)
414 DO i=1,(path_length/2)-1
415 tmp = cmumps_updatescore(weight(path_length),
416 & weight(2*i-1),tup)
417 tmp = cmumps_update_inverse(tmp,weight(2*i),tup)
418 IF(tmp .GT. best_score) THEN
419 best_score = tmp
420 best_beg = cur_el_path
421 ENDIF
422 cur_el_path = cperm(cur_el_path)
423 tmp = cmumps_updatescore(weight(path_length+1),
424 & weight(2*i),tup)
425 tmp = cmumps_update_inverse(tmp,weight(2*i+1),tup)
426 IF(tmp .GT. best_score) THEN
427 best_score = tmp
428 best_beg = cur_el_path
429 ENDIF
430 cur_el_path = cperm(cur_el_path)
431 ENDDO
432 1000 cur_el_path = best_beg
433 DO i=1,(path_length/2)-1
434 num2 = num2+1
435 piv_out(num2) = cur_el_path
436 cur_el_path = cperm(cur_el_path)
437 num2 = num2+1
438 piv_out(num2) = cur_el_path
439 cur_el_path = cperm(cur_el_path)
440 ENDDO
441 numtot = numtot + path_length - 2
442 marked(cur_el_path) = -1
443 ENDIF
444 ENDDO
445 DO i=1,n
446 IF(marked(i) .LT. 0) THEN
447 IF(diag(i) .EQ. 0) THEN
448 piv_out(nlast) = i
449 nlast = nlast - 1
450 ELSE
451 num1 = num1 + 1
452 piv_out(num2+num1) = i
453 numtot = numtot + 1
454 ENDIF
455 ENDIF
456 ENDDO
457 info(2) = numtot
458 info(3) = num1
459 info(4) = num2
460 RETURN
461 END SUBROUTINE cmumps_sym_mwm
462 FUNCTION cmumps_updatescore(A,B,T)
463 IMPLICIT NONE
464 REAL cmumps_updatescore
465 REAL a,b
466 INTEGER t
467 INTEGER sum
468 parameter(sum = 1)
469 IF(t .EQ. sum) THEN
470 cmumps_updatescore = a+b
471 ELSE
472 cmumps_updatescore = a*b
473 ENDIF
474 END FUNCTION cmumps_updatescore
475 FUNCTION cmumps_update_inverse(A,B,T)
476 IMPLICIT NONE
477 REAL cmumps_update_inverse
478 REAL a,b
479 INTEGER t
480 INTEGER sum
481 parameter(sum = 1)
482 IF(t .EQ. sum) THEN
483 cmumps_update_inverse = a-b
484 ELSE
485 cmumps_update_inverse = a/b
486 ENDIF
487 END FUNCTION cmumps_update_inverse
488 FUNCTION cmumps_metric2x2(CUR_EL,CUR_EL_PATH,
489 & SET1,SET2,L1,L2,VAL,DIAG,N,FLAG,FLAGON,T)
490 IMPLICIT NONE
491 REAL cmumps_metric2x2
492 INTEGER cur_el,cur_el_path,l1,l2,n
493 INTEGER set1(l1),set2(l2),diag(n),flag(n)
494 REAL val
495 LOGICAL flagon
496 INTEGER t
497 INTEGER i,inter,merge
498 INTEGER struct,ma47
499 parameter(struct=0,ma47=1)
500 IF(t .EQ. struct) THEN
501 IF(.NOT. flagon) THEN
502 DO i=1,l1
503 flag(set1(i)) = cur_el
504 ENDDO
505 ENDIF
506 inter = 0
507 DO i=1,l2
508 IF(flag(set2(i)) .EQ. cur_el) THEN
509 inter = inter + 1
510 flag(set2(i)) = cur_el_path
511 ENDIF
512 ENDDO
513 merge = l1 + l2 - inter
514 cmumps_metric2x2 = real(inter) / real(merge)
515 ELSE IF (t .EQ. ma47) THEN
516 merge = 3
517 IF(diag(cur_el) .NE. 0) merge = 2
518 IF(diag(cur_el_path) .NE. 0) merge = merge - 2
519 IF(merge .EQ. 0) THEN
520 cmumps_metric2x2 = real(l1+l2-2)
521 cmumps_metric2x2 = -(cmumps_metric2x2**2)/2.0e0
522 ELSE IF(merge .EQ. 1) THEN
523 cmumps_metric2x2 = - real(l1+l2-4) * real(l1-2)
524 ELSE IF(merge .EQ. 2) THEN
525 cmumps_metric2x2 = - real(l1+l2-4) * real(l2-2)
526 ELSE
527 cmumps_metric2x2 = - real(l1-2) * real(l2-2)
528 ENDIF
529 ELSE
530 cmumps_metric2x2 = val
531 ENDIF
532 RETURN
533 END FUNCTION
534 SUBROUTINE cmumps_expand_perm_schur(NA, NCMP,
535 & INVPERM,PERM,
536 & LISTVAR_SCHUR, SIZE_SCHUR, AOTOA)
537 IMPLICIT NONE
538 INTEGER, INTENT(IN):: SIZE_SCHUR, LISTVAR_SCHUR(SIZE_SCHUR)
539 INTEGER, INTENT(IN):: NA, NCMP
540 INTEGER, INTENT(IN):: AOTOA(NCMP), PERM(NCMP)
541 INTEGER, INTENT(OUT):: INVPERM(NA)
542 INTEGER CMP_POS, IO, I, K, IPOS
543 DO cmp_pos=1, ncmp
544 io = perm(cmp_pos)
545 invperm(aotoa(io)) = cmp_pos
546 ENDDO
547 ipos = ncmp
548 DO k =1, size_schur
549 i = listvar_schur(k)
550 ipos = ipos+1
551 invperm(i) = ipos
552 ENDDO
553 RETURN
554 END SUBROUTINE cmumps_expand_perm_schur
555 SUBROUTINE cmumps_gnew_schur
556 & (na, n, nz, irn, icn, iw, lw, ipe, len,
557 & iq, flag, iwfr,
558 & nrorm, niorm, iflag,ierror, icntl,
559 & symmetry, sym, nbqd, avgdens,
560 & keep264, keep265,
561 & listvar_schur, size_schur, atoao, aotoa,
562 & inplace64_graph_copy
563 & )
564 IMPLICIT NONE
565 INTEGER, intent(in) :: NA
566 INTEGER, intent(in) :: N, SYM
567 INTEGER(8), intent(in) :: NZ, LW
568 INTEGER, intent(in) :: ICNTL(60)
569 INTEGER, intent(in) :: IRN(NZ), ICN(NZ)
570 INTEGER, INTENT(IN) :: SIZE_SCHUR, LISTVAR_SCHUR(SIZE_SCHUR)
571 INTEGER, intent(out) :: IERROR, symmetry
572 INTEGER, intent(out) :: NBQD, AvgDens
573 INTEGER, intent(out) :: LEN(N), IW(LW)
574 INTEGER(8), intent(out):: IWFR
575 INTEGER(8), intent(out):: NRORM, NIORM
576 INTEGER(8), intent(out):: IPE(N+1)
577 INTEGER, INTENT(OUT) :: AOTOA(N)
578 INTEGER, INTENT(OUT) :: ATOAO(NA)
579 INTEGER, intent(inout) :: IFLAG, KEEP264
580 INTEGER, intent(in) :: KEEP265
581 INTEGER(8), intent(out):: IQ(N)
582 INTEGER, intent(out) :: FLAG(N)
583 LOGICAL, intent(inout) :: INPLACE64_GRAPH_COPY
584 INTEGER :: MP, MPG, I, J, N1
585 INTEGER :: NBERR, THRESH, IAO
586 INTEGER(8) :: K8, K1, K2, LAST, NDUP
587 INTEGER(8) :: NZOFFA, NDIAGA, L, N8
588 REAL :: RSYM
589 INTRINSIC nint
590 mp = icntl(2)
591 mpg= icntl(3)
592 atoao(1:na) = 0
593 DO i = 1, size_schur
594 atoao(listvar_schur(i)) = -1
595 ENDDO
596 iao = 0
597 DO i= 1, na
598 IF (atoao(i).LT.0) cycle
599 iao = iao +1
600 atoao(i) = iao
601 aotoa(iao) = i
602 ENDDO
603 nzoffa = 0_8
604 ndiaga = 0
605 ierror = 0
606 n8 = int(n,8)
607 DO i=1,n+1
608 ipe(i) = 0_8
609 ENDDO
610 IF (keep264.EQ.0) THEN
611 IF ((sym.EQ.0).AND.(keep265.EQ.-1)) THEN
612 DO k8=1_8,nz
613 i = irn(k8)
614 j = icn(k8)
615 IF ((i.GT.na).OR.(j.GT.na).OR.(i.LT.1)
616 & .OR.(j.LT.1)) THEN
617 ierror = ierror + 1
618 ELSE
619 i = atoao(i)
620 j = atoao(j)
621 IF ((i.LT.0).OR.(j.LT.0)) cycle
622 IF (i.NE.j) THEN
623 ipe(i) = ipe(i) + 1_8
624 nzoffa = nzoffa + 1_8
625 ELSE
626 ndiaga = ndiaga + 1_8
627 ENDIF
628 ENDIF
629 ENDDO
630 ELSE
631 DO k8=1_8,nz
632 i = irn(k8)
633 j = icn(k8)
634 IF ((i.GT.na).OR.(j.GT.na).OR.(i.LT.1)
635 & .OR.(j.LT.1)) THEN
636 ierror = ierror + 1
637 ELSE
638 i = atoao(i)
639 j = atoao(j)
640 IF ((i.LT.0).OR.(j.LT.0)) cycle
641 IF (i.NE.j) THEN
642 ipe(i) = ipe(i) + 1_8
643 ipe(j) = ipe(j) + 1_8
644 nzoffa = nzoffa + 1_8
645 ELSE
646 ndiaga = ndiaga + 1_8
647 ENDIF
648 ENDIF
649 ENDDO
650 ENDIF
651 IF (ierror.GE.1) THEN
652 keep264 = 0
653 ELSE
654 keep264 = 1
655 ENDIF
656 ELSE
657 IF ((sym.EQ.0).AND.(keep265.EQ.-1)) THEN
658 DO k8=1_8,nz
659 i = irn(k8)
660 j = icn(k8)
661 i = atoao(i)
662 j = atoao(j)
663 IF ((i.LT.0).OR.(j.LT.0)) cycle
664 IF (i.EQ.j) THEN
665 ndiaga = ndiaga + 1_8
666 ELSE
667 ipe(i) = ipe(i) + 1_8
668 nzoffa = nzoffa + 1_8
669 ENDIF
670 ENDDO
671 ELSE
672 DO k8=1_8,nz
673 i = irn(k8)
674 j = icn(k8)
675 i = atoao(i)
676 j = atoao(j)
677 IF ((i.LT.0).OR.(j.LT.0)) cycle
678 IF (i.NE.j) THEN
679 ipe(i) = ipe(i) + 1_8
680 ipe(j) = ipe(j) + 1_8
681 nzoffa = nzoffa + 1_8
682 ELSE
683 ndiaga = ndiaga + 1_8
684 ENDIF
685 ENDDO
686 ENDIF
687 ENDIF
688 niorm = nzoffa + 3_8*n8
689 IF (ierror.GE.1) THEN
690 nberr = 0
691 IF (mod(iflag,2) .EQ. 0) iflag = iflag+1
692 IF ((mp.GT.0).AND.(icntl(4).GE.2)) THEN
693 WRITE (mp,99999)
694 DO 70 k8=1_8,nz
695 i = irn(k8)
696 j = icn(k8)
697 IF ((i.GT.na).OR.(j.GT.na).OR.(i.LT.1)
698 & .OR.(j.LT.1)) THEN
699 nberr = nberr + 1
700 IF (nberr.LE.10) THEN
701 IF (mod(k8,10_8).GT.3_8 .OR. mod(k8,10_8).EQ.0_8 .OR.
702 & (10_8.LE.k8 .AND. k8.LE.20_8)) THEN
703 WRITE (mp,'(I16,A,I10,A,I10,A)')
704 & k8,'th entry (in row',i,' and column',j,') ignored'
705 ELSE
706 IF (mod(k8,10_8).EQ.1_8)
707 & WRITE(mp,'(I16,A,I10,A,I10,A)')
708 & k8,'st entry (in row',i,' and column',j,') ignored'
709 IF (mod(k8,10_8).EQ.2_8)
710 & WRITE(mp,'(I16,A,I10,A,I10,A)')
711 & k8,'nd entry (in row',i,' and column',j,') ignored'
712 IF (mod(k8,10_8).EQ.3_8)
713 & WRITE(mp,'(I16,A,I10,A,I10,A)')
714 & k8,'rd entry (in row',i,' and column',j,') ignored'
715 ENDIF
716 ELSE
717 GO TO 100
718 ENDIF
719 ENDIF
720 70 CONTINUE
721 ENDIF
722 ENDIF
723 100 nrorm = niorm - 2_8*n8
724 iq(1) = 1_8
725 n1 = n - 1
726 IF (n1.GT.0) THEN
727 DO i=1,n1
728 iq(i+1) = ipe(i) + iq(i)
729 ENDDO
730 ENDIF
731 last = max(ipe(n)+iq(n)-1,iq(n))
732 flag(1:n) = 0
733 ipe(1:n) = iq(1:n)
734 iw(1:last) = 0
735 iwfr = last + 1_8
736 IF (keep264 .EQ. 0) THEN
737 IF ((sym.EQ.0).AND.(keep265.EQ.-1)) THEN
738 DO k8=1_8,nz
739 i = irn(k8)
740 j = icn(k8)
741 IF ((i.GT.na).OR.(j.GT.na).OR.(i.LT.1)
742 & .OR.(j.LT.1)) cycle
743 i = atoao(i)
744 j = atoao(j)
745 IF ((i.LT.0).OR.(j.LT.0)) cycle
746 IF (i.NE.j) THEN
747 IF ((j.GE.1).AND.(i.LE.n)) THEN
748 iw(iq(i)) = j
749 iq(i) = iq(i) + 1
750 ENDIF
751 ENDIF
752 ENDDO
753 ELSE IF (keep265.EQ.1) THEN
754 DO k8=1_8,nz
755 i = irn(k8)
756 j = icn(k8)
757 IF ((i.GT.na).OR.(j.GT.na).OR.(i.LT.1)
758 & .OR.(j.LT.1)) cycle
759 i = atoao(i)
760 j = atoao(j)
761 IF ((i.LT.0).OR.(j.LT.0)) cycle
762 IF (i.NE.j) THEN
763 IF ((j.GE.1).AND.(i.LE.n)) THEN
764 iw(iq(j)) = i
765 iq(j) = iq(j) + 1
766 iw(iq(i)) = j
767 iq(i) = iq(i) + 1
768 ENDIF
769 ENDIF
770 ENDDO
771 ELSE
772 DO k8=1_8,nz
773 i = irn(k8)
774 j = icn(k8)
775 IF ((i.GT.na).OR.(j.GT.na).OR.(i.LT.1)
776 & .OR.(j.LT.1)) cycle
777 i = atoao(i)
778 j = atoao(j)
779 IF ((i.LT.0).OR.(j.LT.0)) cycle
780 IF (i.NE.j) THEN
781 IF (i.LT.j) THEN
782 IF ((i.GE.1).AND.(j.LE.n)) THEN
783 iw(iq(i)) = -j
784 iq(i) = iq(i) + 1
785 ENDIF
786 ELSE
787 IF ((j.GE.1).AND.(i.LE.n)) THEN
788 iw(iq(j)) = -i
789 iq(j) = iq(j) + 1
790 ENDIF
791 ENDIF
792 ENDIF
793 ENDDO
794 ENDIF
795 ELSE
796 IF ((sym.EQ.0).AND.(keep265.EQ.-1)) THEN
797 DO k8=1_8,nz
798 i = irn(k8)
799 j = icn(k8)
800 i = atoao(i)
801 j = atoao(j)
802 IF ((i.LT.0).OR.(j.LT.0)) cycle
803 IF (i.NE.j) THEN
804 iw(iq(i)) = j
805 iq(i) = iq(i) + 1
806 ENDIF
807 ENDDO
808 ELSE IF (keep265.EQ.1) THEN
809 DO k8=1_8,nz
810 i = irn(k8)
811 j = icn(k8)
812 i = atoao(i)
813 j = atoao(j)
814 IF ((i.LT.0).OR.(j.LT.0)) cycle
815 IF (i.NE.j) THEN
816 iw(iq(j)) = i
817 iq(j) = iq(j) + 1
818 iw(iq(i)) = j
819 iq(i) = iq(i) + 1
820 ENDIF
821 ENDDO
822 ELSE
823 DO k8=1_8,nz
824 i = irn(k8)
825 j = icn(k8)
826 i = atoao(i)
827 j = atoao(j)
828 IF ((i.LT.0).OR.(j.LT.0)) cycle
829 IF (i.NE.j) THEN
830 IF (i.LT.j) THEN
831 iw(iq(i)) = -j
832 iq(i) = iq(i) + 1
833 ELSE
834 iw(iq(j)) = -i
835 iq(j) = iq(j) + 1
836 ENDIF
837 ENDIF
838 ENDDO
839 ENDIF
840 ENDIF
841 IF (keep265.EQ.0) THEN
842 ndup = 0_8
843 DO i=1,n
844 k1 = ipe(i)
845 k2 = iq(i) - 1_8
846 IF (k1.GT.k2) THEN
847 len(i) = 0
848 ELSE
849 DO k8=k1,k2
850 j = -iw(k8)
851 IF (j.LE.0) EXIT
852 IF (flag(j).EQ.i) THEN
853 ndup = ndup + 1_8
854 iw(k8) = 0
855 ELSE
856 l = iq(j)
857 iq(j) = l + 1
858 iw(l) = i
859 iw(k8) = j
860 flag(j) = i
861 ENDIF
862 END DO
863 len(i) = int(iq(i) - ipe(i))
864 ENDIF
865 ENDDO
866 IF (ndup.NE.0_8) THEN
867 iwfr = 1_8
868 DO i=1,n
869 IF (len(i).EQ.0) THEN
870 ipe(i) = iwfr
871 cycle
872 ENDIF
873 k1 = ipe(i)
874 k2 = k1 + len(i) - 1
875 l = iwfr
876 ipe(i) = iwfr
877 DO 270 k8=k1,k2
878 IF (iw(k8).NE.0) THEN
879 iw(iwfr) = iw(k8)
880 iwfr = iwfr + 1_8
881 ENDIF
882 270 CONTINUE
883 len(i) = int(iwfr - l)
884 ENDDO
885 ENDIF
886 ipe(n+1) = ipe(n) + int(len(n),8)
887 iwfr = ipe(n+1)
888 ELSE
889 ipe(1) = 1_8
890 DO i = 1, n
891 len(i) = int(iq(i) - ipe(i))
892 ENDDO
893 DO i = 1, n
894 ipe(i+1) = ipe(i) + int(len(i),8)
895 ENDDO
896 iwfr = ipe(n+1)
897 ENDIF
898 symmetry = 100
899 IF (sym.EQ.0) THEN
900 rsym = real(ndiaga+2_8*nzoffa - (iwfr-1_8))/
901 & real(nzoffa+ndiaga)
902 IF ((keep265.EQ.0) .AND. (nzoffa - (iwfr-1_8)).EQ.0_8) THEN
903 ENDIF
904 symmetry = nint(100.0e0*rsym)
905 IF ((mpg .GT. 0).AND.(icntl(4).GE.2) )
906 & write(mpg,'(A,A,I5)')
907 & ' Case of Schur:',
908 & ' structural symmetry (in percent) of interior block=',
909 & symmetry
910 IF (mp.GT.0 .AND. mpg.NE.mp.AND. (icntl(4).GE.2) )
911 & write(mp,'(A,A,I5)')
912 & ' Case of Schur:',
913 & ' structural symmetry (in percent) of interior block=',
914 & symmetry
915 ELSE
916 symmetry = 100
917 ENDIF
918 inplace64_graph_copy = (lw.GE.2*(iwfr-1))
919 avgdens = nint(real(iwfr-1_8)/real(n))
920 thresh = avgdens*50 - avgdens/10 + 1
921 nbqd = 0
922 IF (n.GT.2) THEN
923 DO i= 1, n
924 j = max(len(i),1)
925 IF (j.GT.thresh) nbqd = nbqd+1
926 ENDDO
927 ENDIF
928 IF (mpg .GT. 0.AND.(icntl(4).GE.2))
929 & write(mpg,'(A,1I5)')
930 & ' Average density of rows/columns =', avgdens
931 IF (mp.GT.0 .AND. mpg.NE.mp.AND.(icntl(4).GE.2))
932 & write(mpg,'(A,1I5)')
933 & ' Average density of rows/columns =', avgdens
934 RETURN
93599999 FORMAT (/'*** Warning message from analysis routine ***')
936 END SUBROUTINE cmumps_gnew_schur
937 SUBROUTINE cmumps_get_perm_from_pe(N,PE,INVPERM,NFILS,WORK)
938 IMPLICIT NONE
939 INTEGER N
940 INTEGER PE(N),INVPERM(N),NFILS(N),WORK(N)
941 INTEGER I,FATHER,STKLEN,STKPOS,PERMPOS,CURVAR
942 nfils = 0
943 DO i=1,n
944 father = -pe(i)
945 IF(father .NE. 0) nfils(father) = nfils(father) + 1
946 ENDDO
947 stklen = 0
948 permpos = 1
949 DO i=1,n
950 IF(nfils(i) .EQ. 0) THEN
951 stklen = stklen + 1
952 work(stklen) = i
953 invperm(i) = permpos
954 permpos = permpos + 1
955 ENDIF
956 ENDDO
957 DO stkpos = 1,stklen
958 curvar = work(stkpos)
959 father = -pe(curvar)
960 DO
961 IF(father .EQ. 0) EXIT
962 IF(nfils(father) .EQ. 1) THEN
963 invperm(father) = permpos
964 father = -pe(father)
965 permpos = permpos + 1
966 ELSE
967 nfils(father) = nfils(father) - 1
968 EXIT
969 ENDIF
970 ENDDO
971 ENDDO
972 RETURN
973 END SUBROUTINE cmumps_get_perm_from_pe
974 SUBROUTINE cmumps_get_elim_tree(N,PE,NV,WORK)
975 IMPLICIT NONE
976 INTEGER N
977 INTEGER PE(N),NV(N),WORK(N)
978 INTEGER I,FATHER,LEN,NEWSON,NEWFATHER
979 DO i=1,n
980 IF(nv(i) .GT. 0) cycle
981 len = 1
982 work(len) = i
983 father = -pe(i)
984 DO
985 IF(nv(father) .GT. 0) THEN
986 newson = father
987 EXIT
988 ENDIF
989 len = len + 1
990 work(len) = father
991 nv(father) = 1
992 father = -pe(father)
993 ENDDO
994 newfather = -pe(father)
995 pe(work(len)) = -newfather
996 pe(newson) = -work(1)
997 ENDDO
998 END SUBROUTINE cmumps_get_elim_tree
cmumps_ldlt_compress
subroutine cmumps_ldlt_compress(n, nz, irn, icn, piv, ncmp, iw, lw, ipe, len, iq, flag, icmp, iwfr, ierror, keep, keep8, icntl, inplace64_graph_copy)
Definition cana_LDLT_preprocess.F:127
cmumps_set_constraints
subroutine cmumps_set_constraints(n, piv, frere, fils, nfsiz, ikeep, ncst, keep, keep8, rowsca)
Definition cana_LDLT_preprocess.F:18
cmumps_get_perm_from_pe
subroutine cmumps_get_perm_from_pe(n, pe, invperm, nfils, work)
Definition cana_LDLT_preprocess.F:938
cmumps_get_elim_tree
subroutine cmumps_get_elim_tree(n, pe, nv, work)
Definition cana_LDLT_preprocess.F:975
cmumps_updatescore
real function cmumps_updatescore(a, b, t)
Definition cana_LDLT_preprocess.F:463
cmumps_gnew_schur
subroutine cmumps_gnew_schur(na, n, nz, irn, icn, iw, lw, ipe, len, iq, flag, iwfr, nrorm, niorm, iflag, ierror, icntl, symmetry, sym, nbqd, avgdens, keep264, keep265, listvar_schur, size_schur, atoao, aotoa, inplace64_graph_copy)
Definition cana_LDLT_preprocess.F:564
cmumps_expand_permutation
subroutine cmumps_expand_permutation(n, ncmp, n11, n22, piv, invperm, perm)
Definition cana_LDLT_preprocess.F:90
cmumps_sym_mwm
subroutine cmumps_sym_mwm(n, ne, ip, irn, scaling, lsc, cperm, diag, icntl, weight, marked, flag, piv_out, info)
Definition cana_LDLT_preprocess.F:275
cmumps_update_inverse
real function cmumps_update_inverse(a, b, t)
Definition cana_LDLT_preprocess.F:476
cmumps_expand_perm_schur
subroutine cmumps_expand_perm_schur(na, ncmp, invperm, perm, listvar_schur, size_schur, aotoa)
Definition cana_LDLT_preprocess.F:537
cmumps_metric2x2
real function cmumps_metric2x2(cur_el, cur_el_path, set1, set2, l1, l2, val, diag, n, flag, flagon, t)
Definition cana_LDLT_preprocess.F:490
max
#define max(a, b)
Definition macros.h:21
merge
subroutine merge(x, itab, itabm1, cmerge, imerge, imerge2, iadmerge2, nmerge_tot)
Definition merge.F:36
cmumps_struc_def
Definition cmumps_struc_def.F:14