← Back to branch summary

~ubuntu-branches/ubuntu/lucid/mumps/lucid

« back to all changes in this revision

Viewing changes to src/dmumps_part7.F

  • Committer: Bazaar Package Importer
  • Author(s): Adam C. Powell, IV
  • Date: 2009-12-07 17:56:51 UTC
  • mfrom: (1.1.1 upstream)
  • Revision ID: james.westby@ubuntu.com-20091207175651-ftogh061hebcqzty
Tags: 4.9.2.dfsg-1
http://bugs.debian.org/554159
http://bugs.debian.org/557699
http://bugs.debian.org/555759
* New upstream release (closes: #554159).
* Changed -lblas to -lblas-3gf in Makefile.*.inc (closes: #557699).
* Linking tests to shared instead of static libs (closes: #555759).

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/dmumps_part7.F
1
1
C
2
 
C  This file is part of MUMPS 4.8.4, built on Mon Dec 15 15:31:38 UTC 2008
 
2
C  This file is part of MUMPS 4.9.2, built on Thu Nov  5 07:05:08 UTC 2009
3
3
C
4
4
C
5
5
C  This version of MUMPS is provided to you free of charge. It is public
6
6
C  domain, based on public domain software developed during the Esprit IV
7
7
C  European project PARASOL (1996-1999) by CERFACS, ENSEEIHT-IRIT and RAL.
8
8
C  Since this first public domain version in 1999, the developments are
9
 
C  supported by the following institutions: CERFACS, ENSEEIHT-IRIT, and
10
 
C  INRIA.
 
9
C  supported by the following institutions: CERFACS, CNRS, INPT(ENSEEIHT)-
 
10
C  IRIT, and INRIA.
11
11
C
12
 
C  Main contributors are Patrick Amestoy, Iain Duff, Abdou Guermouche,
13
 
C  Jacko Koster, Jean-Yves L'Excellent, and Stephane Pralet.
 
12
C  Current development team includes Patrick Amestoy, Alfredo Buttari,
 
13
C  Abdou Guermouche, Jean-Yves L'Excellent, Bora Ucar.
14
14
C
15
15
C  Up-to-date copies of the MUMPS package can be obtained
16
16
C  from the Web pages:
23
23
C
24
24
C  User documentation of any code that uses this software can
25
25
C  include this complete notice. You can acknowledge (using
26
 
C  references [1], [2], and [3]) the contribution of this package
 
26
C  references [1] and [2]) the contribution of this package
27
27
C  in any scientific publication dependent upon the use of the
28
28
C  package. You shall use reasonable endeavours to notify
29
29
C  the authors of the package of this publication.
30
30
C
31
 
C   [1] P. R. Amestoy, I. S. Duff and  J.-Y. L'Excellent,
32
 
C   Multifrontal parallel distributed symmetric and unsymmetric solvers,
33
 
C   in Comput. Methods in Appl. Mech. Eng., 184,  501-520 (2000).
34
 
C
35
 
C   [2] P. R. Amestoy, I. S. Duff, J. Koster and  J.-Y. L'Excellent,
 
31
C   [1] P. R. Amestoy, I. S. Duff, J. Koster and  J.-Y. L'Excellent,
36
32
C   A fully asynchronous multifrontal solver using distributed dynamic
37
33
C   scheduling, SIAM Journal of Matrix Analysis and Applications,
38
34
C   Vol 23, No 1, pp 15-41 (2001).
39
35
C
40
 
C   [3] P. R. Amestoy and A. Guermouche and J.-Y. L'Excellent and
 
36
C   [2] P. R. Amestoy and A. Guermouche and J.-Y. L'Excellent and
41
37
C   S. Pralet, Hybrid scheduling for the parallel solution of linear
42
38
C   systems. Parallel Computing Vol 32 (2), pp 136-156 (2006).
43
39
C
52
48
      INTEGER ICNTL(40), KEEP(500), MPG
53
49
      KEEP(111)=ICNTL(25)
54
50
      IF (KEEP(111) < -1 .OR.
55
 
     *      KEEP(111).GT.KEEP(112)+KEEP(17)) THEN
 
51
     &      KEEP(111).GT.KEEP(112)+KEEP(17)) THEN
56
52
          KEEP(111)=0
57
53
      ENDIF
58
54
      IF (KEEP(19).EQ.0.AND.KEEP(110).EQ.0) THEN
59
55
        IF(KEEP(111).NE.0.AND.MPG.GT.0) THEN
60
56
          WRITE( MPG,'(A)')
61
 
     *'** Warning: ICNTL(25) option disabled because'
 
57
     &'** Warning: ICNTL(25) option disabled because'
62
58
          WRITE( MPG,'(A)')
63
 
     *'** null space was not required during factorization'
 
59
     &'** null space was not required during factorization'
64
60
        ENDIF
65
61
        KEEP(111)=0
66
62
      ENDIF
67
63
      IF (ICNTL(9).NE.1) THEN
68
64
        IF (KEEP(111).NE.0.AND. MPG.GT.0) THEN
69
65
          WRITE(MPG,'(A)')
70
 
     *'** Warning: ICNTL(25) option disabled because'
 
66
     &'** Warning: ICNTL(25) option disabled because'
71
67
          WRITE( MPG,'(A)')
72
 
     *'** it is not available for the transposed system'
 
68
     &'** it is not available for the transposed system'
73
69
        ENDIF
74
70
      ENDIF
75
71
      RETURN
92
88
      RETURN
93
89
      END SUBROUTINE DMUMPS_636
94
90
      SUBROUTINE DMUMPS_279( PHASE, MBLOCK, NBLOCK, 
95
 
     *           SIZE_ROOT_ARG,
96
 
     *           LOCAL_M, LOCAL_N, ROOT_OWNER, KEEP,KEEP8,
97
 
     *           LIWK_RR, LWK_RR )
 
91
     &           SIZE_ROOT_ARG,
 
92
     &           LOCAL_M, LOCAL_N, ROOT_OWNER, KEEP,KEEP8,
 
93
     &           LIWK_RR, LWK_RR )
98
94
      IMPLICIT NONE
99
95
      INTEGER, INTENT(IN) :: PHASE, SIZE_ROOT_ARG
100
96
      INTEGER, INTENT(IN) :: MBLOCK, NBLOCK, LOCAL_M, LOCAL_N
123
119
        END IF
124
120
      ENDIF
125
121
      RETURN
126
 
      END
 
122
      END SUBROUTINE DMUMPS_279
127
123
      SUBROUTINE DMUMPS_333(N,PERM,X,RN01)
128
124
      INTEGER N,PERM(N),I
129
125
      DOUBLE PRECISION RN01(N),X(N)
134
130
      X(I)=RN01(I)
135
131
200   CONTINUE
136
132
      RETURN
137
 
      END
 
133
      END SUBROUTINE DMUMPS_333
138
134
      SUBROUTINE DMUMPS_334(N,PERM,X,RN01)
139
135
      INTEGER N,PERM(N),I
140
136
      DOUBLE PRECISION RN01(N),X(N)
145
141
      X(I)=RN01(I)
146
142
200   CONTINUE
147
143
      RETURN
148
 
      END
149
 
      DOUBLE PRECISION FUNCTION DMUMPS_104(N,AA,WR01)
150
 
      INTEGER N
151
 
      DOUBLE PRECISION AA(N*N),WR01(N)
152
 
      INTEGER I,J,JSTRT
153
 
      DOUBLE PRECISION ZERO
154
 
      DOUBLE PRECISION VALUE
155
 
      PARAMETER(ZERO=0.0D0)
156
 
      VALUE=ZERO
157
 
      CALL DMUMPS_117(N,ZERO,WR01,1)
158
 
      JSTRT=1
159
 
      DO I=1,N
160
 
        DO J=JSTRT,JSTRT+N-1
161
 
          WR01(I)=WR01(I)+abs(AA(J))
162
 
        END DO
163
 
        JSTRT=JSTRT+N
164
 
      END DO
165
 
      DO I=1,N
166
 
        VALUE=max(abs(VALUE),abs(WR01(I)))
167
 
      END DO
168
 
      DMUMPS_104=VALUE
169
 
      RETURN
170
 
      END
 
144
      END SUBROUTINE DMUMPS_334
171
145
      SUBROUTINE DMUMPS_117(N,ALPHA,DX,INCX)
172
146
      INTEGER N,INCX
173
147
      DOUBLE PRECISION ALPHA,DX(*)
202
176
        END DO
203
177
        RETURN
204
178
      END IF
205
 
      END
 
179
      END SUBROUTINE DMUMPS_117
206
180
      SUBROUTINE DMUMPS_146( MYID, root, N, IROOT,
207
 
     *           COMM, IW, LIW, IFREE,
208
 
     *           A, LA, PTRAST, PTLUST_S, PTRFAC,
209
 
     *           STEP, INFO, LDLT, QR,
210
 
     *           WK, LWK, KEEP,KEEP8)
 
181
     &           COMM, IW, LIW, IFREE,
 
182
     &           A, LA, PTRAST, PTLUST_S, PTRFAC,
 
183
     &           STEP, INFO, LDLT, QR,
 
184
     &           WK, LWK, KEEP,KEEP8)
211
185
      IMPLICIT NONE
212
186
      INCLUDE 'dmumps_root.h'
213
187
      INCLUDE 'mpif.h'
214
188
      TYPE ( DMUMPS_ROOT_STRUC ) :: root
215
 
      INTEGER N, IROOT, COMM, LIW, LA, MYID, LWK, LIWK, IFREE
 
189
      INTEGER N, IROOT, COMM, LIW, MYID, LIWK, IFREE
 
190
      INTEGER(8) :: LA
 
191
      INTEGER(8) :: LWK
216
192
      DOUBLE PRECISION WK( LWK )
217
193
      INTEGER KEEP(500)
218
194
      INTEGER*8 KEEP8(150)
219
 
      INTEGER PTRAST(KEEP(28)), PTLUST_S(KEEP(28)), PTRFAC(KEEP(28)),
220
 
     *STEP(N), IW( LIW )
 
195
      INTEGER(8) :: PTRAST(KEEP(28))
 
196
      INTEGER(8) :: PTRFAC(KEEP(28))
 
197
      INTEGER PTLUST_S(KEEP(28)), STEP(N), IW( LIW )
221
198
      INTEGER INFO( 2 ), LDLT, QR
222
199
      DOUBLE PRECISION A( LA )
223
 
      INTEGER IOLDPS, IAPOS
 
200
      INTEGER IOLDPS
 
201
      INTEGER(8) :: IAPOS
224
202
      INTEGER LOCAL_M, LOCAL_N, LPIV, IERR, allocok,i
225
203
      INCLUDE 'mumps_headers.h'
226
204
        IF ( .NOT. root%yes ) RETURN
227
205
        IF ( KEEP(60) .NE. 0 ) THEN
228
206
          IF ((LDLT == 1 .OR. LDLT == 2) .AND. KEEP(60) == 3 ) THEN
229
207
            CALL DMUMPS_320( WK, root%MBLOCK,
230
 
     *      root%MYROW, root%MYCOL, root%NPROW, root%NPCOL,
231
 
     *      root%SCHUR_POINTER(1),
232
 
     *      root%SCHUR_LLD, root%SCHUR_NLOC,
233
 
     *      root%TOT_ROOT_SIZE, MYID, COMM )
 
208
     &      root%MYROW, root%MYCOL, root%NPROW, root%NPCOL,
 
209
     &      root%SCHUR_POINTER(1),
 
210
     &      root%SCHUR_LLD, root%SCHUR_NLOC,
 
211
     &      root%TOT_ROOT_SIZE, MYID, COMM )
234
212
          ENDIF
235
213
        RETURN
236
214
        ENDIF
254
232
          CALL MUMPS_ABORT()
255
233
        END IF
256
234
        CALL DESCINIT( root%DESCRIPTOR, root%TOT_ROOT_SIZE,
257
 
     *      root%TOT_ROOT_SIZE, root%MBLOCK, root%NBLOCK,
258
 
     *      0, 0, root%CNTXT_BLACS, LOCAL_M, IERR )
 
235
     &      root%TOT_ROOT_SIZE, root%MBLOCK, root%NBLOCK,
 
236
     &      0, 0, root%CNTXT_BLACS, LOCAL_M, IERR )
259
237
        IF ( LDLT.EQ.2 ) THEN
260
238
            IF(root%MBLOCK.NE.root%NBLOCK) THEN
261
239
              WRITE(*,*) ' Error: symmetrization only works for'
262
240
              WRITE(*,*) ' square block sizes, MBLOCK/NBLOCK=',
263
 
     *        root%MBLOCK, root%NBLOCK
 
241
     &        root%MBLOCK, root%NBLOCK
264
242
              CALL MUMPS_ABORT()
265
243
            END IF
266
 
            IF ( LWK .LT.
267
 
     *        min(root%MBLOCK * root%NBLOCK,root%TOT_ROOT_SIZE*
268
 
     *        root%TOT_ROOT_SIZE ) ) THEN
269
 
              WRITE(*,*) 'Not enough workspace for symmetrization.'
270
 
              CALL MUMPS_ABORT()
 
244
            IF ( LWK .LT. min(
 
245
     &           int(root%MBLOCK,8) * int(root%NBLOCK,8),
 
246
     &           int(root%TOT_ROOT_SIZE,8)* int(root%TOT_ROOT_SIZE,8 )
 
247
     &         )) THEN
 
248
               WRITE(*,*) 'Not enough workspace for symmetrization.'
 
249
               CALL MUMPS_ABORT()
271
250
            END IF
272
251
            CALL DMUMPS_320( WK, root%MBLOCK,
273
 
     *      root%MYROW, root%MYCOL, root%NPROW, root%NPCOL,
274
 
     *      A( IAPOS ), LOCAL_M, LOCAL_N,
275
 
     *      root%TOT_ROOT_SIZE, MYID, COMM )
 
252
     &      root%MYROW, root%MYCOL, root%NPROW, root%NPCOL,
 
253
     &      A( IAPOS ), LOCAL_M, LOCAL_N,
 
254
     &      root%TOT_ROOT_SIZE, MYID, COMM )
276
255
        END IF
277
256
        IF (LDLT.EQ.0.OR.LDLT.EQ.2) THEN
278
257
          CALL PDGETRF( root%TOT_ROOT_SIZE, root%TOT_ROOT_SIZE,
279
 
     *      A( IAPOS ),
280
 
     *      1, 1, root%DESCRIPTOR, root%IPIV(1), IERR )
 
258
     &      A( IAPOS ),
 
259
     &      1, 1, root%DESCRIPTOR, root%IPIV(1), IERR )
281
260
          IF ( IERR .GT. 0 ) THEN
282
261
              INFO(1)=-10
283
262
              INFO(2)=IERR-1
284
263
          END IF
285
264
        ELSE
286
265
          CALL PDPOTRF('L',root%TOT_ROOT_SIZE,A(IAPOS),
287
 
     *      1,1,root%DESCRIPTOR,IERR)
 
266
     &      1,1,root%DESCRIPTOR,IERR)
288
267
            IF ( IERR .GT. 0 ) THEN
289
 
              INFO(1)=-10
 
268
              INFO(1)=-40
290
269
              INFO(2)=IERR-1
291
270
            END IF
292
271
        END IF
293
272
        RETURN
294
273
      END SUBROUTINE DMUMPS_146
295
274
      SUBROUTINE DMUMPS_556(
296
 
     *     N,PIV,FRERE,FILS,NFSIZ,IKEEP,
297
 
     *     NCST,KEEP,KEEP8,id)
 
275
     &     N,PIV,FRERE,FILS,NFSIZ,IKEEP,
 
276
     &     NCST,KEEP,KEEP8,id)
298
277
      USE DMUMPS_STRUC_DEF
299
278
      IMPLICIT NONE
300
279
      TYPE (DMUMPS_STRUC) :: id
312
291
         P2 = PIV(I+1)
313
292
         K1 = IKEEP(P1,1)
314
293
         IF(K1 .GT. 0) THEN
315
 
            V1 = (abs(id%A(K1)*(id%ROWSCA(P1)**2)).GE.1.0D-1)
 
294
            V1 = (abs(id%A(K1))*(id%ROWSCA(P1)**2).GE.1.0D-1)
316
295
         ELSE
317
296
            V1 = .FALSE.
318
297
         ENDIF
319
298
         K2 = IKEEP(P2,1)
320
299
         IF(K2 .GT. 0) THEN
321
 
            V2 = (abs(id%A(K2)*(id%ROWSCA(P2)**2)).GE.1.0D-1)
 
300
            V2 = (abs(id%A(K2))*(id%ROWSCA(P2)**2).GE.1.0D-1)
322
301
         ELSE
323
302
            V2 = .FALSE.
324
303
         ENDIF
365
344
      ENDDO
366
345
      END SUBROUTINE DMUMPS_556
367
346
      SUBROUTINE DMUMPS_550(N,NCMP,N11,N22,PIV,
368
 
     *     INVPERM,PERM)
 
347
     &     INVPERM,PERM)
369
348
      IMPLICIT NONE
370
349
      INTEGER N11,N22,N,NCMP
371
350
      INTEGER, intent(in) :: PIV(N),PERM(N)
399
378
      RETURN
400
379
      END SUBROUTINE DMUMPS_550
401
380
      SUBROUTINE DMUMPS_547(
402
 
     *     N,NZ, IRN, ICN, PIV,
403
 
     *     NCMP, IW, LW, IPE, LEN, IQ, 
404
 
     *     FLAG, ICMP, IWFR,
405
 
     *     IERROR, KEEP,KEEP8, ICNTL)
 
381
     &     N,NZ, IRN, ICN, PIV,
 
382
     &     NCMP, IW, LW, IPE, LEN, IQ, 
 
383
     &     FLAG, ICMP, IWFR,
 
384
     &     IERROR, KEEP,KEEP8, ICNTL)
406
385
      IMPLICIT NONE
407
386
      INTEGER N,NZ,NCMP,LW,IWFR,IERROR
408
387
      INTEGER ICNTL(40),KEEP(500)
444
423
         I = ICMP(I)
445
424
         J = ICMP(J)
446
425
         IF ((I.GT.N).OR.(J.GT.N).OR.(I.LT.1)
447
 
     *        .OR.(J.LT.1)) THEN
 
426
     &        .OR.(J.LT.1)) THEN
448
427
            IERROR = IERROR + 1
449
428
         ELSE
450
429
            IF (I.NE.J) THEN
541
520
      RETURN
542
521
      END SUBROUTINE DMUMPS_547
543
522
      SUBROUTINE DMUMPS_551(
544
 
     *     N, NE, IP, IRN, SCALING,LSC,CPERM, DIAG,
545
 
     *     ICNTL, WEIGHT,MARKED,FLAG,
546
 
     *     PIV_OUT, INFO)
 
523
     &     N, NE, IP, IRN, SCALING,LSC,CPERM, DIAG,
 
524
     &     ICNTL, WEIGHT,MARKED,FLAG,
 
525
     &     PIV_OUT, INFO)
547
526
      IMPLICIT NONE
548
527
      INTEGER N, NE, ICNTL(10), INFO(10),LSC,LWEIGHT
549
528
      INTEGER CPERM(N),PIV_OUT(N),IP(N+1),IRN(NE), DIAG(N)
553
532
      INTEGER I,J,K,L,CUR_EL,CUR_EL_PATH,CUR_EL_PATH_NEXT,BEST_BEG
554
533
      INTEGER L1,L2,PTR_SET1,PTR_SET2,TUP,T22,INTER,MERGE
555
534
      DOUBLE PRECISION BEST_SCORE,CUR_VAL,TMP,VAL
556
 
      DOUBLE PRECISION INITSCORE, DMUMPS_UPDATESCORE, 
557
 
     *     DMUMPS_UPDATE_INVERSE, DMUMPS_METRIC2x2
 
535
      DOUBLE PRECISION INITSCORE, DMUMPS_739, 
 
536
     &     DMUMPS_740, DMUMPS_741
558
537
      LOGICAL VRAI,FAUX,MAX_CARD_DIAG,USE_SCALING
559
538
      INTEGER JOB_DEF,DONE
560
539
      INTEGER SUM,PROD,STRUCT,MA47,MAGNITUDE
561
 
      INTEGER ZERO,ONE
 
540
      DOUBLE PRECISION ZERO,ONE
562
541
      PARAMETER (JOB_DEF = 5, DONE = -1,
563
 
     *     SUM = 1, PROD = 2, STRUCT=1, MA47=2, MAGNITUDE=3,
564
 
     *     VRAI = .TRUE., FAUX = .FALSE.)
 
542
     &     SUM = 1, PROD = 2, STRUCT=1, MA47=2, MAGNITUDE=3,
 
543
     &     VRAI = .TRUE., FAUX = .FALSE.)
565
544
      PARAMETER(ZERO = 0.0D0, ONE = 1.0D0)
566
545
      MAX_CARD_DIAG = .TRUE.
567
546
      NUM1 = 0
585
564
      ENDIF
586
565
      IF(ICNTL(2) .GT. 2 .OR. ICNTL(2) .LE. 0) THEN
587
566
         WRITE(*,*)
588
 
     *        'ERROR: WRONG VALUE FOR ICNTL(2) = ',ICNTL(2)
 
567
     &        'ERROR: WRONG VALUE FOR ICNTL(2) = ',ICNTL(2)
589
568
         INFO(1) = -1
590
569
         RETURN
591
570
      ENDIF
592
571
      T22 = ICNTL(1)
593
572
      IF(ICNTL(1) .LT. 0 .OR. ICNTL(1) .GT. 2) THEN
594
573
         WRITE(*,*)
595
 
     *        'ERROR: WRONG VALUE FOR ICNTL(1) = ',ICNTL(1)
 
574
     &        'ERROR: WRONG VALUE FOR ICNTL(1) = ',ICNTL(1)
596
575
         INFO(1) = -1
597
576
         RETURN
598
577
      ENDIF
620
599
         IF(USE_SCALING) THEN
621
600
            VAL = -SCALING(CUR_EL_PATH) - SCALING(CUR_EL+N)
622
601
         ENDIF
623
 
         CUR_VAL = DMUMPS_METRIC2x2(
624
 
     *        CUR_EL,CUR_EL_PATH,
625
 
     *        IRN(PTR_SET1),IRN(PTR_SET2),
626
 
     *        L1,L2,
627
 
     *        VAL,DIAG,N,FLAG,FAUX,T22)
 
602
         CUR_VAL = DMUMPS_741(
 
603
     &        CUR_EL,CUR_EL_PATH,
 
604
     &        IRN(PTR_SET1),IRN(PTR_SET2),
 
605
     &        L1,L2,
 
606
     &        VAL,DIAG,N,FLAG,FAUX,T22)
628
607
         WEIGHT(PATH_LENGTH+1) = 
629
 
     *        DMUMPS_UPDATESCORE(WEIGHT(1),CUR_VAL,TUP)
 
608
     &        DMUMPS_739(WEIGHT(1),CUR_VAL,TUP)
630
609
         DO
631
610
            IF(CUR_EL_PATH .EQ. CUR_EL) EXIT
632
611
            PATH_LENGTH = PATH_LENGTH+1
638
617
            PTR_SET2 = IP(CUR_EL_PATH_NEXT)
639
618
            IF(USE_SCALING) THEN
640
619
               VAL = -SCALING(CUR_EL_PATH_NEXT) 
641
 
     *              - SCALING(CUR_EL_PATH+N)
 
620
     &              - SCALING(CUR_EL_PATH+N)
642
621
            ENDIF
643
 
            CUR_VAL = DMUMPS_METRIC2x2(
644
 
     *           CUR_EL_PATH,CUR_EL_PATH_NEXT,
645
 
     *           IRN(PTR_SET1),IRN(PTR_SET2),
646
 
     *           L1,L2,
647
 
     *           VAL,DIAG,N,FLAG,VRAI,T22)
 
622
            CUR_VAL = DMUMPS_741(
 
623
     &           CUR_EL_PATH,CUR_EL_PATH_NEXT,
 
624
     &           IRN(PTR_SET1),IRN(PTR_SET2),
 
625
     &           L1,L2,
 
626
     &           VAL,DIAG,N,FLAG,VRAI,T22)
648
627
            WEIGHT(PATH_LENGTH+1) = 
649
 
     *           DMUMPS_UPDATESCORE(WEIGHT(PATH_LENGTH-1),CUR_VAL,TUP)
 
628
     &           DMUMPS_739(WEIGHT(PATH_LENGTH-1),CUR_VAL,TUP)
650
629
            CUR_EL_PATH = CUR_EL_PATH_NEXT
651
630
         ENDDO
652
631
         IF(mod(PATH_LENGTH,2) .EQ. 1) THEN
683
662
            BEST_SCORE = WEIGHT(PATH_LENGTH-1)
684
663
            CUR_EL_PATH = CPERM(CUR_EL)
685
664
            DO I=1,(PATH_LENGTH/2)-1
686
 
               TMP = DMUMPS_UPDATESCORE(WEIGHT(PATH_LENGTH),
687
 
     *              WEIGHT(2*I-1),TUP)
688
 
               TMP = DMUMPS_UPDATE_INVERSE(TMP,WEIGHT(2*I),TUP)
 
665
               TMP = DMUMPS_739(WEIGHT(PATH_LENGTH),
 
666
     &              WEIGHT(2*I-1),TUP)
 
667
               TMP = DMUMPS_740(TMP,WEIGHT(2*I),TUP)
689
668
               IF(TMP .GT. BEST_SCORE) THEN
690
669
                  BEST_SCORE = TMP
691
670
                  BEST_BEG = CUR_EL_PATH
692
671
               ENDIF
693
672
               CUR_EL_PATH = CPERM(CUR_EL_PATH)
694
 
               TMP = DMUMPS_UPDATESCORE(WEIGHT(PATH_LENGTH+1),
695
 
     *              WEIGHT(2*I),TUP)
696
 
               TMP = DMUMPS_UPDATE_INVERSE(TMP,WEIGHT(2*I+1),TUP)
 
673
               TMP = DMUMPS_739(WEIGHT(PATH_LENGTH+1),
 
674
     &              WEIGHT(2*I),TUP)
 
675
               TMP = DMUMPS_740(TMP,WEIGHT(2*I+1),TUP)
697
676
               IF(TMP .GT. BEST_SCORE) THEN
698
677
                  BEST_SCORE = TMP
699
678
                  BEST_BEG = CUR_EL_PATH
730
709
      INFO(4) = NUM2
731
710
      RETURN
732
711
      END SUBROUTINE DMUMPS_551
733
 
      FUNCTION DMUMPS_UPDATESCORE(A,B,T)
734
 
      IMPLICIT NONE
735
 
      DOUBLE PRECISION DMUMPS_UPDATESCORE
736
 
      DOUBLE PRECISION A,B
737
 
      INTEGER T
738
 
      INTEGER SUM,PROD
739
 
      PARAMETER(SUM = 1,PROD = 2)
740
 
      IF(T .EQ. SUM) THEN
741
 
         DMUMPS_UPDATESCORE = A+B
742
 
      ELSE
743
 
         DMUMPS_UPDATESCORE = A*B
744
 
      ENDIF
745
 
      END FUNCTION DMUMPS_UPDATESCORE
746
 
      FUNCTION DMUMPS_UPDATE_INVERSE(A,B,T)
747
 
      IMPLICIT NONE
748
 
      DOUBLE PRECISION DMUMPS_UPDATE_INVERSE
749
 
      DOUBLE PRECISION A,B
750
 
      INTEGER T
751
 
      INTEGER SUM,PROD
752
 
      PARAMETER(SUM = 1,PROD = 2)
753
 
      IF(T .EQ. SUM) THEN
754
 
         DMUMPS_UPDATE_INVERSE = A-B
755
 
      ELSE
756
 
         DMUMPS_UPDATE_INVERSE = A/B
757
 
      ENDIF
758
 
      END FUNCTION DMUMPS_UPDATE_INVERSE
759
 
      FUNCTION DMUMPS_METRIC2x2(CUR_EL,CUR_EL_PATH,
760
 
     *     SET1,SET2,L1,L2,VAL,DIAG,N,FLAG,FLAGON,T)
761
 
      IMPLICIT NONE
762
 
      DOUBLE PRECISION DMUMPS_METRIC2x2
 
712
      FUNCTION DMUMPS_739(A,B,T)
 
713
      IMPLICIT NONE
 
714
      DOUBLE PRECISION DMUMPS_739
 
715
      DOUBLE PRECISION A,B
 
716
      INTEGER T
 
717
      INTEGER SUM,PROD
 
718
      PARAMETER(SUM = 1,PROD = 2)
 
719
      IF(T .EQ. SUM) THEN
 
720
         DMUMPS_739 = A+B
 
721
      ELSE
 
722
         DMUMPS_739 = A*B
 
723
      ENDIF
 
724
      END FUNCTION DMUMPS_739
 
725
      FUNCTION DMUMPS_740(A,B,T)
 
726
      IMPLICIT NONE
 
727
      DOUBLE PRECISION DMUMPS_740
 
728
      DOUBLE PRECISION A,B
 
729
      INTEGER T
 
730
      INTEGER SUM,PROD
 
731
      PARAMETER(SUM = 1,PROD = 2)
 
732
      IF(T .EQ. SUM) THEN
 
733
         DMUMPS_740 = A-B
 
734
      ELSE
 
735
         DMUMPS_740 = A/B
 
736
      ENDIF
 
737
      END FUNCTION DMUMPS_740
 
738
      FUNCTION DMUMPS_741(CUR_EL,CUR_EL_PATH,
 
739
     &     SET1,SET2,L1,L2,VAL,DIAG,N,FLAG,FLAGON,T)
 
740
      IMPLICIT NONE
 
741
      DOUBLE PRECISION DMUMPS_741
763
742
      INTEGER CUR_EL,CUR_EL_PATH,L1,L2,N
764
743
      INTEGER SET1(L1),SET2(L2),DIAG(N),FLAG(N)
765
744
      DOUBLE PRECISION VAL
782
761
            ENDIF
783
762
         ENDDO
784
763
         MERGE = L1 + L2 - INTER
785
 
         DMUMPS_METRIC2x2 = dble(INTER) / dble(MERGE)
 
764
         DMUMPS_741 = dble(INTER) / dble(MERGE)
786
765
      ELSE IF (T .EQ. MA47) THEN
787
766
         MERGE = 3
788
767
         IF(DIAG(CUR_EL) .NE. 0) MERGE = 2
789
768
         IF(DIAG(CUR_EL_PATH) .NE. 0) MERGE = MERGE - 2
790
769
         IF(MERGE .EQ. 0) THEN
791
 
            DMUMPS_METRIC2x2 = L1+L2-2
792
 
            DMUMPS_METRIC2x2 = -(DMUMPS_METRIC2x2**2)/2
 
770
            DMUMPS_741 = dble(L1+L2-2)
 
771
            DMUMPS_741 = -(DMUMPS_741**2)/2.0D0
793
772
         ELSE IF(MERGE .EQ. 1) THEN
794
 
            DMUMPS_METRIC2x2 = - dble(L1+L2-4) * dble(L1-2)
 
773
            DMUMPS_741 = - dble(L1+L2-4) * dble(L1-2)
795
774
         ELSE IF(MERGE .EQ. 2) THEN
796
 
            DMUMPS_METRIC2x2 = - dble(L1+L2-4) * dble(L2-2)
 
775
            DMUMPS_741 = - dble(L1+L2-4) * dble(L2-2)
797
776
         ELSE
798
 
            DMUMPS_METRIC2x2 = - dble(L1-2) * dble(L2-2)
 
777
            DMUMPS_741 = - dble(L1-2) * dble(L2-2)
799
778
         ENDIF
800
779
      ELSE
801
 
         DMUMPS_METRIC2x2 = VAL
 
780
         DMUMPS_741 = VAL
802
781
      ENDIF
803
782
      RETURN
804
783
      END FUNCTION 
805
784
      SUBROUTINE DMUMPS_622(NA, NCMP,
806
 
     *      INVPERM,PERM, 
807
 
     *      LISTVAR_SCHUR, SIZE_SCHUR, AOTOA)
 
785
     &      INVPERM,PERM, 
 
786
     &      LISTVAR_SCHUR, SIZE_SCHUR, AOTOA)
808
787
      IMPLICIT NONE
809
788
      INTEGER, INTENT(IN):: SIZE_SCHUR, LISTVAR_SCHUR(SIZE_SCHUR)
810
789
      INTEGER, INTENT(IN):: NA, NCMP
824
803
      RETURN
825
804
      END SUBROUTINE DMUMPS_622
826
805
      SUBROUTINE DMUMPS_623
827
 
     * (NA,N,NZ, IRN, ICN, IW, LW, IPE, LEN,
828
 
     * IQ, FLAG, IWFR,
829
 
     * NRORM, NIORM, IFLAG,IERROR, ICNTL, 
830
 
     * symmetry, SYM, MedDens, NBQD, AvgDens, 
831
 
     * LISTVAR_SCHUR, SIZE_SCHUR, ATOAO, AOTOA)
 
806
     & (NA,N,NZ, IRN, ICN, IW, LW, IPE, LEN,
 
807
     & IQ, FLAG, IWFR,
 
808
     & NRORM, NIORM, IFLAG,IERROR, ICNTL, 
 
809
     & symmetry, SYM, MedDens, NBQD, AvgDens, 
 
810
     & LISTVAR_SCHUR, SIZE_SCHUR, ATOAO, AOTOA)
832
811
      IMPLICIT NONE
833
812
      INTEGER, INTENT(IN)  :: NA,N,NZ,LW
834
813
      INTEGER, INTENT(IN)  :: SIZE_SCHUR, LISTVAR_SCHUR(SIZE_SCHUR)
835
814
      INTEGER, INTENT(IN)  :: IRN(NZ), ICN(NZ) 
836
815
      INTEGER, INTENT(IN)  :: ICNTL(40), SYM
837
 
      INTEGER, INTENT(OUT) :: IFLAG,IERROR,NRORM,NIORM,IWFR
 
816
      INTEGER, INTENT(INOUT) :: IFLAG
 
817
      INTEGER, INTENT(OUT) :: IERROR,NRORM,NIORM,IWFR
838
818
      INTEGER, INTENT(OUT) :: AOTOA(N)
839
819
      INTEGER, INTENT(OUT) :: ATOAO(NA)  
840
820
      INTEGER, INTENT(OUT) :: LEN(N), IPE(N+1)
868
848
        I = IRN(K)
869
849
        J = ICN(K)
870
850
        IF ((I.GT.NA).OR.(J.GT.NA).OR.(I.LT.1)
871
 
     *                          .OR.(J.LT.1)) THEN
 
851
     &                          .OR.(J.LT.1)) THEN
872
852
           IERROR = IERROR + 1
873
853
        ELSE
874
854
          I = ATOAO(I)
893
873
           I = IRN(K)
894
874
           J = ICN(K)
895
875
           IF ((I.GT.NA).OR.(J.GT.NA).OR.(I.LT.1)
896
 
     *                            .OR.(J.LT.1)) THEN
 
876
     &                            .OR.(J.LT.1)) THEN
897
877
            NBERR = NBERR + 1
898
878
            IF (NBERR.LE.10)  THEN
899
879
               IF (mod(K,10).GT.3 .OR. mod(K,10).EQ.0 .OR.
900
 
     *             (10.LE.K .AND. K.LE.20)) THEN
 
880
     &             (10.LE.K .AND. K.LE.20)) THEN
901
881
                 WRITE (MP,'(I8,A,I8,A,I8,A)')
902
 
     *             K,'th entry (in row',I,' and column',J,') ignored'
 
882
     &             K,'th entry (in row',I,' and column',J,') ignored'
903
883
               ELSE
904
884
                 IF (mod(K,10).EQ.1) WRITE(MP,'(I8,A,I8,A,I8,A)')
905
 
     *             K,'st entry (in row',I,' and column',J,') ignored'
 
885
     &             K,'st entry (in row',I,' and column',J,') ignored'
906
886
                 IF (mod(K,10).EQ.2) WRITE(MP,'(I8,A,I8,A,I8,A)')
907
 
     *             K,'nd entry (in row',I,' and column',J,') ignored'
 
887
     &             K,'nd entry (in row',I,' and column',J,') ignored'
908
888
                 IF (mod(K,10).EQ.3) WRITE(MP,'(I8,A,I8,A,I8,A)')
909
 
     *             K,'rd entry (in row',I,' and column',J,') ignored'
 
889
     &             K,'rd entry (in row',I,' and column',J,') ignored'
910
890
               ENDIF
911
891
            ELSE
912
892
               GO TO 100
932
912
         I = IRN(K)
933
913
         J = ICN(K)
934
914
         IF ((I.GT.NA).OR.(J.GT.NA).OR.(I.LT.1)
935
 
     *                          .OR.(J.LT.1)) CYCLE
 
915
     &                          .OR.(J.LT.1)) CYCLE
936
916
         I = ATOAO(I)
937
917
         J = ATOAO(J)
938
918
         IF ((I.LT.0).OR.(J.LT.0)) CYCLE  
998
978
      IWFR = IPE(N+1)
999
979
      IF (SYM.EQ.0) THEN
1000
980
      RSYM =  dble(NDIAGA+2*NZOFFA - (IWFR-1))/
1001
 
     *            dble(NZOFFA+NDIAGA) 
1002
 
      symmetry = nint (100.0*RSYM)
 
981
     &            dble(NZOFFA+NDIAGA) 
 
982
      symmetry = nint (100.0D0*RSYM)
1003
983
         IF (MPG .GT. 0)
1004
984
     &  write(MPG,'(A,I5)') 
1005
985
     &  ' ... Structural symmetry (in percent)=', symmetry
1076
1056
      ENDDO
1077
1057
      RETURN
1078
1058
      END SUBROUTINE DMUMPS_549
 
1059
      SUBROUTINE DMUMPS_548(N,PE,NV,WORK)
 
1060
      IMPLICIT NONE
 
1061
      INTEGER N
 
1062
      INTEGER PE(N),NV(N),WORK(N)
 
1063
      INTEGER I,FATHER,LEN,K,NEWSON,NEWFATHER
 
1064
      DO I=1,N
 
1065
         IF(NV(I) .GT. 0) CYCLE
 
1066
         LEN = 1
 
1067
         WORK(LEN) = I
 
1068
         FATHER = -PE(I)
 
1069
         DO
 
1070
            IF(NV(FATHER) .GT. 0) THEN
 
1071
               NEWSON = FATHER
 
1072
               EXIT
 
1073
            ENDIF
 
1074
            LEN = LEN + 1
 
1075
            WORK(LEN) = FATHER
 
1076
            NV(FATHER) = 1
 
1077
            FATHER = -PE(FATHER)
 
1078
         ENDDO
 
1079
         NEWFATHER = -PE(FATHER)
 
1080
         PE(WORK(LEN)) = -NEWFATHER
 
1081
         PE(NEWSON) = -WORK(1)
 
1082
      ENDDO      
 
1083
      END SUBROUTINE DMUMPS_548