File: dfac_determinant.F

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