1C
2C  This file is part of MUMPS 5.1.2, released
3C  on Mon Oct  2 07:37:01 UTC 2017
4C
5C
6C  Copyright 1991-2017 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
7C  University of Bordeaux.
8C
9C  This version of MUMPS is provided to you free of charge. It is
10C  released under the CeCILL-C license:
11C  http://www.cecill.info/licences/Licence_CeCILL-C_V1-en.html
12C
13      SUBROUTINE DMUMPS_UPDATEDETER(PIV, DETER, NEXP)
14      IMPLICIT NONE
15      DOUBLE PRECISION, intent(in) :: PIV
16      DOUBLE PRECISION, intent(inout) :: DETER
17      INTEGER, intent(inout) :: NEXP
18      DETER=DETER*fraction(PIV)
19      NEXP=NEXP+exponent(PIV)+exponent(DETER)
20      DETER=fraction(DETER)
21      RETURN
22      END SUBROUTINE DMUMPS_UPDATEDETER
23      SUBROUTINE DMUMPS_UPDATEDETER_SCALING(PIV, DETER, NEXP)
24      IMPLICIT NONE
25      DOUBLE PRECISION, intent(in) :: PIV
26      DOUBLE PRECISION, intent(inout) :: DETER
27      INTEGER, intent(inout) :: NEXP
28      DETER=DETER*fraction(PIV)
29      NEXP=NEXP+exponent(PIV)+exponent(DETER)
30      DETER=fraction(DETER)
31      RETURN
32      END SUBROUTINE DMUMPS_UPDATEDETER_SCALING
33      SUBROUTINE DMUMPS_GETDETER2D(BLOCK_SIZE,IPIV,
34     &                      MYROW, MYCOL, NPROW, NPCOL,
35     &                      A, LOCAL_M, LOCAL_N, N, MYID,
36     &                      DETER,NEXP,SYM)
37      IMPLICIT NONE
38      INTEGER, intent (in)    :: SYM
39      INTEGER, intent (inout) :: NEXP
40      DOUBLE PRECISION, intent (inout) :: DETER
41      INTEGER, intent (in)    :: BLOCK_SIZE, NPROW, NPCOL,
42     &                           LOCAL_M, LOCAL_N, N
43      INTEGER, intent (in)    :: MYROW, MYCOL, MYID, IPIV(LOCAL_M)
44      DOUBLE PRECISION, intent(in) :: A(*)
45      INTEGER  I,IMX,DI,NBLOCK,IBLOCK,ILOC,JLOC,
46     &         ROW_PROC,COL_PROC, K
47      DI = LOCAL_M + 1
48      NBLOCK = ( N - 1 ) / BLOCK_SIZE
49      DO IBLOCK = 0, NBLOCK
50        ROW_PROC = mod( IBLOCK, NPROW )
51        IF ( MYROW.EQ.ROW_PROC ) THEN
52          COL_PROC = mod( IBLOCK, NPCOL )
53          IF ( MYCOL.EQ.COL_PROC ) THEN
54            ILOC = ( IBLOCK / NPROW ) * BLOCK_SIZE
55            JLOC = ( IBLOCK / NPCOL ) * BLOCK_SIZE
56            I   =   ILOC + JLOC *  LOCAL_M + 1
57            IMX = min(ILOC+BLOCK_SIZE,LOCAL_M)
58     &            + (min(JLOC+BLOCK_SIZE,LOCAL_N)-1)*LOCAL_M
59     &            + 1
60            K=1
61            DO WHILE ( I .LT. IMX )
62              CALL DMUMPS_UPDATEDETER(A(I),DETER,NEXP)
63              IF (SYM.NE.1) THEN
64                IF (IPIV(ILOC+K) .NE. IBLOCK*BLOCK_SIZE+K) THEN
65                  DETER = -DETER
66                ENDIF
67              ENDIF
68              K = K + 1
69              I = I + DI
70            END DO
71          END IF
72        END IF
73      END DO
74      RETURN
75      END SUBROUTINE DMUMPS_GETDETER2D
76      SUBROUTINE DMUMPS_DETER_REDUCTION(
77     &           COMM, DETER_IN, NEXP_IN,
78     &           DETER_OUT, NEXP_OUT, NPROCS)
79      IMPLICIT NONE
80      INTEGER, intent(in) :: COMM, NPROCS
81      DOUBLE PRECISION, intent(in) :: DETER_IN
82      INTEGER,intent(in) :: NEXP_IN
83      DOUBLE PRECISION,intent(out):: DETER_OUT
84      INTEGER,intent(out):: NEXP_OUT
85      INTEGER            :: IERR_MPI
86      EXTERNAL DMUMPS_DETERREDUCE_FUNC
87      INTEGER TWO_SCALARS_TYPE, DETERREDUCE_OP
88      DOUBLE PRECISION :: INV(2)
89      DOUBLE PRECISION :: OUTV(2)
90      INCLUDE 'mpif.h'
91      IF (NPROCS .EQ. 1) THEN
92        DETER_OUT = DETER_IN
93        NEXP_OUT  = NEXP_IN
94        RETURN
95      ENDIF
96      CALL MPI_TYPE_CONTIGUOUS(2, MPI_DOUBLE_PRECISION,
97     &                         TWO_SCALARS_TYPE,
98     &                         IERR_MPI)
99      CALL MPI_TYPE_COMMIT(TWO_SCALARS_TYPE, IERR_MPI)
100      CALL MPI_OP_CREATE(DMUMPS_DETERREDUCE_FUNC,
101     &                   .TRUE.,
102     &                   DETERREDUCE_OP,
103     &                   IERR_MPI)
104      INV(1)=DETER_IN
105      INV(2)=dble(NEXP_IN)
106      CALL MPI_ALLREDUCE( INV, OUTV, 1, TWO_SCALARS_TYPE,
107     &                    DETERREDUCE_OP, COMM, IERR_MPI)
108      CALL MPI_OP_FREE(DETERREDUCE_OP, IERR_MPI)
109      CALL MPI_TYPE_FREE(TWO_SCALARS_TYPE, IERR_MPI)
110      DETER_OUT = OUTV(1)
111      NEXP_OUT  = int(OUTV(2))
112      RETURN
113      END SUBROUTINE DMUMPS_DETER_REDUCTION
114      SUBROUTINE DMUMPS_DETERREDUCE_FUNC(INV, INOUTV, NEL, DATATYPE)
115      IMPLICIT NONE
116#if defined(WORKAROUNDINTELILP64MPI2INTEGER)
117      INTEGER(4), INTENT(IN)    :: NEL, DATATYPE
118#else
119      INTEGER, INTENT(IN)    :: NEL, DATATYPE
120#endif
121      DOUBLE PRECISION, INTENT(IN)    :: INV    ( 2 * NEL )
122      DOUBLE PRECISION, INTENT(INOUT) :: INOUTV ( 2 * NEL )
123      INTEGER I, TMPEXPIN, TMPEXPINOUT
124      DO I = 1, NEL
125        TMPEXPIN    = int(INV   (I*2))
126        TMPEXPINOUT = int(INOUTV(I*2))
127        CALL DMUMPS_UPDATEDETER(INV(I*2-1),
128     &                          INOUTV(I*2-1),
129     &                          TMPEXPINOUT)
130        TMPEXPINOUT = TMPEXPINOUT + TMPEXPIN
131        INOUTV(I*2) = dble(TMPEXPINOUT)
132      ENDDO
133      RETURN
134      END SUBROUTINE DMUMPS_DETERREDUCE_FUNC
135      SUBROUTINE DMUMPS_DETER_SQUARE(DETER, NEXP)
136      IMPLICIT NONE
137      INTEGER, intent (inout) :: NEXP
138      DOUBLE PRECISION, intent (inout) :: DETER
139      DETER=DETER*DETER
140      NEXP=NEXP+NEXP
141      RETURN
142      END SUBROUTINE DMUMPS_DETER_SQUARE
143      SUBROUTINE DMUMPS_DETER_SCALING_INVERSE(DETER, NEXP)
144      IMPLICIT NONE
145      INTEGER, intent (inout) :: NEXP
146      DOUBLE PRECISION, intent (inout) :: DETER
147      DETER=1.0D0/DETER
148      NEXP=-NEXP
149      RETURN
150      END SUBROUTINE DMUMPS_DETER_SCALING_INVERSE
151      SUBROUTINE DMUMPS_DETER_SIGN_PERM(DETER, N, VISITED, PERM)
152      IMPLICIT NONE
153      DOUBLE PRECISION, intent(inout) :: DETER
154      INTEGER, intent(in)    :: N
155      INTEGER, intent(inout) :: VISITED(N)
156      INTEGER, intent(in)    :: PERM(N)
157      INTEGER I, J, K
158      K = 0
159      DO I = 1, N
160        IF (VISITED(I) .GT. N) THEN
161          VISITED(I)=VISITED(I)-N-N-1
162          CYCLE
163        ENDIF
164        J = PERM(I)
165        DO WHILE (J.NE.I)
166          VISITED(J) = VISITED(J) + N + N + 1
167          K = K + 1
168          J = PERM(J)
169        ENDDO
170      ENDDO
171      IF (mod(K,2).EQ.1) THEN
172        DETER = -DETER
173      ENDIF
174      RETURN
175      END SUBROUTINE DMUMPS_DETER_SIGN_PERM
176