sfac_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 smumps_updatedeter(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 smumps_updatedeter
      SUBROUTINE smumps_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 smumps_updatedeter_scaling
      SUBROUTINE smumps_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
      REAL, intent (inout) :: DETER
      INTEGER, intent (in)    :: BLOCK_SIZE, NPROW, NPCOL,
     &                           local_m, local_n, n
      INTEGER, intent (in)    :: MYROW, MYCOL, MYID, IPIV(LOCAL_M)
      REAL, 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 smumps_updatedeter(a(i),deter,nexp)
              IF (sym.EQ.1) THEN
               CALL smumps_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 smumps_getdeter2d
      SUBROUTINE smumps_deter_reduction(
     &           COMM, DETER_IN, NEXP_IN,
     &           DETER_OUT, NEXP_OUT, NPROCS)
      IMPLICIT NONE
      INTEGER, intent(in) :: COMM, NPROCS
      REAL, intent(in) :: DETER_IN
      INTEGER,intent(in) :: NEXP_IN
      REAL,intent(out):: DETER_OUT
      INTEGER,intent(out):: NEXP_OUT
      INTEGER            :: IERR_MPI
      EXTERNAL smumps_deterreduce_func
      INTEGER TWO_SCALARS_TYPE, DETERREDUCE_OP
      REAL :: INV(2)
      REAL :: 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_real,
     &                         two_scalars_type,
     &                         ierr_mpi)
      CALL mpi_type_commit(two_scalars_type, ierr_mpi)
      CALL mpi_op_create(smumps_deterreduce_func,  
     &                   .true.,             
     &                   deterreduce_op,     
     &                   ierr_mpi)
      inv(1)=deter_in
      inv(2)=real(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 smumps_deter_reduction
      SUBROUTINE smumps_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
      REAL, INTENT(IN)    :: INV    ( 2 * NEL )
      REAL, INTENT(INOUT) :: INOUTV ( 2 * NEL )
      INTEGER I, TMPEXPIN, TMPEXPINOUT
      DO i = 1, nel
        tmpexpin    = int(inv   (i*2))
        tmpexpinout = int(inoutv(i*2))
        CALL smumps_updatedeter(inv(i*2-1), 
     &                          inoutv(i*2-1), 
     &                          tmpexpinout)   
        tmpexpinout = tmpexpinout + tmpexpin
        inoutv(i*2) = real(tmpexpinout)
      ENDDO
      RETURN
      END SUBROUTINE smumps_deterreduce_func
      SUBROUTINE smumps_deter_square(DETER, NEXP)
      IMPLICIT NONE
      INTEGER, intent (inout) :: NEXP
      REAL, intent (inout) :: DETER
      deter=deter*deter
      nexp=nexp+nexp
      RETURN
      END SUBROUTINE smumps_deter_square
      SUBROUTINE smumps_deter_scaling_inverse(DETER, NEXP)
      IMPLICIT NONE
      INTEGER, intent (inout) :: NEXP
      REAL, intent (inout) :: DETER
      deter=1.0e0/deter
      nexp=-nexp
      RETURN
      END SUBROUTINE smumps_deter_scaling_inverse
      SUBROUTINE smumps_deter_sign_perm(DETER, N, VISITED, PERM)
      IMPLICIT NONE
      REAL, 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 smumps_deter_sign_perm
      SUBROUTINE smumps_par_root_minmax_piv_upd (
     &                      BLOCK_SIZE,IPIV,
     &                      MYROW, MYCOL, NPROW, NPCOL,
     &                      A, LOCAL_M, LOCAL_N, N, MYID,
     &                      DKEEP, KEEP, SYM)
       USE smumps_fac_front_aux_m, 
     &      only : smumps_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)
      REAL, 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 smumps_update_minmax_pivot 
     &           ( abspivot,
     &             dkeep, keep, .false.)
              k = k + 1
              i = i + di
            END DO
          END IF
        END IF
      END DO
      RETURN
      END SUBROUTINE smumps_par_root_minmax_piv_upd