cfac_determinant.F

Go to the documentation of this file.

C
C  This file is part of MUMPS 5.5.1, released
C  on Tue Jul 12 13:17:24 UTC 2022
C
C
C  Copyright 1991-2022 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
C  Mumps Technologies, 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  (see doc/CeCILL-C_V1-en.txt, doc/CeCILL-C_V1-fr.txt, and
C  https://cecill.info/licences/Licence_CeCILL-C_V1-en.html)
C
      SUBROUTINE cmumps_updatedeter(PIV, DETER, NEXP)
      IMPLICIT NONE
      COMPLEX, intent(in) :: PIV
      COMPLEX, intent(inout) :: DETER
      INTEGER, intent(inout) :: NEXP
      REAL R_PART, C_PART
      INTEGER NEXP_LOC
      deter=deter*piv
      r_part=real(deter)
      c_part=aimag(deter)
      nexp_loc = exponent(abs(r_part)+abs(c_part))
      nexp = nexp + nexp_loc
      r_part=scale(r_part, -nexp_loc)
      c_part=scale(c_part, -nexp_loc)
      deter=cmplx(r_part,c_part,kind=kind(deter))
      RETURN
      END SUBROUTINE cmumps_updatedeter
      SUBROUTINE cmumps_updatedeter_scaling(PIV, DETER, NEXP)
      IMPLICIT NONE
      REAL, intent(in) :: PIV
      REAL, intent(inout) :: DETER
      INTEGER, intent(inout) :: NEXP
      deter=deter*fraction(piv)
      nexp=nexp+exponent(piv)+exponent(deter)
      deter=fraction(deter)
      RETURN
      END SUBROUTINE cmumps_updatedeter_scaling
      SUBROUTINE cmumps_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
      COMPLEX, intent (inout) :: DETER
      INTEGER, intent (in)    :: BLOCK_SIZE, NPROW, NPCOL,
     &                           local_m, local_n, n
      INTEGER, intent (in)    :: MYROW, MYCOL, MYID, IPIV(LOCAL_M)
      COMPLEX, 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 cmumps_updatedeter(a(i),deter,nexp)
              IF (sym.EQ.1) THEN
               CALL cmumps_updatedeter(a(i),deter,nexp)
              ENDIF
              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 cmumps_getdeter2d
      SUBROUTINE cmumps_deter_reduction(
     &           COMM, DETER_IN, NEXP_IN,
     &           DETER_OUT, NEXP_OUT, NPROCS)
      IMPLICIT NONE
      INTEGER, intent(in) :: COMM, NPROCS
      COMPLEX, intent(in) :: DETER_IN
      INTEGER,intent(in) :: NEXP_IN
      COMPLEX,intent(out):: DETER_OUT
      INTEGER,intent(out):: NEXP_OUT
      INTEGER            :: IERR_MPI
      EXTERNAL cmumps_deterreduce_func
      INTEGER TWO_SCALARS_TYPE, DETERREDUCE_OP
      COMPLEX :: INV(2)
      COMPLEX :: 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_complex,
     &                         two_scalars_type,
     &                         ierr_mpi)
      CALL mpi_type_commit(two_scalars_type, ierr_mpi)
      CALL mpi_op_create(cmumps_deterreduce_func,  
     &                   .true.,             
     &                   deterreduce_op,     
     &                   ierr_mpi)
      inv(1)=deter_in
      inv(2)=cmplx(nexp_in,kind=kind(inv))
      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 cmumps_deter_reduction
      SUBROUTINE cmumps_deterreduce_func(INV, INOUTV, NEL, DATATYPE)
      IMPLICIT NONE
#if defined(WORKAROUNDINTELILP64MPI2INTEGER) || defined(WORKAROUNDILP64MPICUSTOMREDUCE)
      INTEGER(4), INTENT(IN)    :: NEL, DATATYPE
#else
      INTEGER, INTENT(IN)    :: NEL, DATATYPE
#endif
      COMPLEX, INTENT(IN)    :: INV    ( 2 * NEL )
      COMPLEX, INTENT(INOUT) :: INOUTV ( 2 * NEL )
      INTEGER I, TMPEXPIN, TMPEXPINOUT
      DO i = 1, nel
        tmpexpin    = int(inv(i*2))
        tmpexpinout = int(inoutv(i*2))
        CALL cmumps_updatedeter(inv(i*2-1), 
     &                          inoutv(i*2-1), 
     &                          tmpexpinout)   
        tmpexpinout = tmpexpinout + tmpexpin
        inoutv(i*2) = cmplx(tmpexpinout,kind=kind(inoutv))
      ENDDO
      RETURN
      END SUBROUTINE cmumps_deterreduce_func
      SUBROUTINE cmumps_deter_square(DETER, NEXP)
      IMPLICIT NONE
      INTEGER, intent (inout) :: NEXP
      COMPLEX, intent (inout) :: DETER
      deter=deter*deter
      nexp=nexp+nexp
      RETURN
      END SUBROUTINE cmumps_deter_square
      SUBROUTINE cmumps_deter_scaling_inverse(DETER, NEXP)
      IMPLICIT NONE
      INTEGER, intent (inout) :: NEXP
      REAL, intent (inout) :: DETER
      deter=1.0e0/deter
      nexp=-nexp
      RETURN
      END SUBROUTINE cmumps_deter_scaling_inverse
      SUBROUTINE cmumps_deter_sign_perm(DETER, N, VISITED, PERM)
      IMPLICIT NONE
      COMPLEX, 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 cmumps_deter_sign_perm
      SUBROUTINE cmumps_par_root_minmax_piv_upd (
     &                      BLOCK_SIZE,IPIV,
     &                      MYROW, MYCOL, NPROW, NPCOL,
     &                      A, LOCAL_M, LOCAL_N, N, MYID,
     &                      DKEEP, KEEP, SYM)
       USE cmumps_fac_front_aux_m, 
     &      only : cmumps_update_minmax_pivot 
      IMPLICIT NONE
      INTEGER, intent (in)    :: BLOCK_SIZE, NPROW, NPCOL,
     &                           local_m, local_n, n, sym
      INTEGER, intent (in)    :: MYROW, MYCOL, MYID, IPIV(LOCAL_M)
      COMPLEX, intent(in) :: A(*)
      REAL, INTENT(INOUT) :: DKEEP(230)
      INTEGER, INTENT(IN) :: KEEP(500)
      INTEGER  I,IMX,DI,NBLOCK,IBLOCK,ILOC,JLOC,
     &         ROW_PROC,COL_PROC, K
      REAL :: ABSPIVOT
      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 )
             IF (sym.NE.1) THEN
              abspivot =  abs(a(i))
             ELSE
              abspivot =  abs(a(i)*a(i))
             ENDIF
             CALL cmumps_update_minmax_pivot 
     &           ( abspivot,
     &             dkeep, keep, .false.)
              k = k + 1
              i = i + di
            END DO
          END IF
        END IF
      END DO
      RETURN
      END SUBROUTINE cmumps_par_root_minmax_piv_upd