zfac_driver.F File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	zmumps_fac_driver (id)
subroutine	zmumps_print_allocated_mem (prok, prokg, print_maxavg, mp, mpg, info16, infog18, infog19, nslaves, irank, keep)
subroutine	zmumps_avgmax_stat8 (prokg, mpg, val, nslaves, print_maxavg, comm, msg)
subroutine	zmumps_extract_schur_redrhs (id)

Function/Subroutine Documentation

zmumps_avgmax_stat8()

subroutine zmumps_avgmax_stat8	(	logical, intent(in)	prokg,
		integer, intent(in)	mpg,
		integer(8), intent(in)	val,
		integer, intent(in)	nslaves,
		logical, intent(in)	print_maxavg,
		integer, intent(in)	comm,
		character*48	msg )

Definition at line 3877 of file zfac_driver.F.

      IMPLICIT NONE
      include 'mpif.h'
      LOGICAL, intent(in) :: PROKG
      INTEGER, intent(in) :: MPG
      INTEGER(8), intent(in) :: VAL
      INTEGER, intent(in) :: NSLAVES
      LOGICAL, intent(in) :: PRINT_MAXAVG
      INTEGER, intent(in) :: COMM
      CHARACTER*48 MSG 
C  Local
      INTEGER(8) MAX_VAL
      INTEGER IERR, MASTER
      DOUBLE PRECISION LOC_VAL, AVG_VAL
      parameter(master=0)
C
      CALL mumps_reducei8( val, max_val, mpi_max, master, comm)
      loc_val = dble(val)/dble(nslaves)
      CALL mpi_reduce( loc_val, avg_val, 1, mpi_double_precision,
     &                 mpi_sum, master, comm, ierr )
      IF (prokg) THEN
        IF (print_maxavg) THEN
          WRITE(mpg,100) " Average", msg, int(avg_val,8)
        ELSE
          WRITE(mpg,110)  msg, max_val
        ENDIF
      ENDIF
      RETURN
 100  FORMAT(a8,a48,i18)
 110  FORMAT(a48,i18)

zmumps_extract_schur_redrhs()

subroutine zmumps_extract_schur_redrhs

(

type(zmumps_struc)

id

)

Definition at line 3910 of file zfac_driver.F.

      USE zmumps_struc_def
      IMPLICIT NONE
C
C  Purpose
C  =======
C
C     Extract the Schur and possibly also the reduced right-hand side
C     (if Fwd in facto) from the processor working on Schur and copy
C     it into the user datastructures id%SCHUR and id%REDRHS on the host.
C     This routine assumes that the integer list of the Schur has not
C     been permuted and still corresponds to LISTVAR_SCHUR.
C
C     If the Schur is centralized, the master of the Schur holds the
C     Schur and possibly also the reduced right-hand side.
C     If the Schur is distribued (already built in user's datastructure),
C     then the master of the Schur may hold the reduced right-hand side,
C     in which case it is available in root%RHS_CNTR_MASTER_ROOT.
C     
      TYPE(ZMUMPS_STRUC) :: id
C
C  Local variables
C  ===============
C
      include 'mpif.h'
      include 'mumps_tags.h'
      include 'mumps_headers.h'
      INTEGER :: STATUS(MPI_STATUS_SIZE)
      INTEGER :: IERR
      INTEGER, PARAMETER :: MASTER = 0
      INTEGER :: ID_SCHUR, SIZE_SCHUR, LD_SCHUR, IB, BL4
      INTEGER(4) :: I4 ! 32-bit even in 64-bit version
      INTEGER :: ROW_LENGTH, I
      INTEGER(8) :: SURFSCHUR8, BL8, SHIFT8
      INTEGER(8) :: ISCHUR_SRC, ISCHUR_DEST, ISCHUR_SYM, ISCHUR_UNS
C
C  External functions
C  ==================
C
      INTEGER MUMPS_PROCNODE
      EXTERNAL mumps_procnode
C     Quick return in case factorization did not terminate correctly
      IF (id%INFO(1) .LT. 0) RETURN
C     Quick return if Schur option off
      IF (id%KEEP(60) .EQ. 0) RETURN
C     Get Schur id
      id_schur =mumps_procnode(
     &    id%PROCNODE_STEPS(id%STEP(max(id%KEEP(20),id%KEEP(38)))),
     &    id%KEEP(199))
      IF ( id%KEEP( 46 )  .NE. 1 ) THEN
        id_schur = id_schur + 1
      END IF
C     Get size of Schur
      IF (id%MYID.EQ.id_schur) THEN
        IF (id%KEEP(60).EQ.1) THEN
C         Sequential Schur
          ld_schur =
     &    id%IS(id%PTLUST_S(id%STEP(id%KEEP(20)))+2+id%KEEP(ixsz))
          size_schur = ld_schur - id%KEEP(253)
        ELSE
C         Parallel Schur
          ld_schur   = -999999 ! not used
          size_schur = id%root%TOT_ROOT_SIZE
        ENDIF
      ELSE IF (id%MYID .EQ. master) THEN
        size_schur = id%KEEP(116)
        ld_schur = -44444 ! Not used
      ELSE
C       Proc is not concerned with Schur, return
        RETURN
      ENDIF
      surfschur8 = int(size_schur,8)*int(size_schur,8)
C     =================================
C     Case of parallel Schur: if REDRHS
C     was requested, obtain it directly
C     from id%root%RHS_CNTR_MASTER_ROOT
C     =================================
      IF (id%KEEP(60) .GT. 1) THEN
        IF (id%KEEP(221).EQ.1 .AND. id%KEEP(252).GT.0) THEN
          DO i = 1, id%KEEP(253)
            IF (id_schur.EQ.master) THEN ! Necessarily = id%MYID
              CALL zcopy(size_schur,
     &             id%root%RHS_CNTR_MASTER_ROOT((i-1)*size_schur+1), 1,
     &             id%REDRHS((i-1)*id%LREDRHS+1), 1)
            ELSE
              IF (id%MYID.EQ.id_schur) THEN
C               Send
                CALL mpi_send(
     &             id%root%RHS_CNTR_MASTER_ROOT((i-1)*size_schur+1),
     &             size_schur,
     &             mpi_double_complex,
     &             master, tag_schur,
     &             id%COMM, ierr )
              ELSE ! MYID.EQ.MASTER
C               Receive
                CALL mpi_recv( id%REDRHS((i-1)*id%LREDRHS+1),
     &             size_schur,
     &             mpi_double_complex, id_schur, tag_schur,
     &             id%COMM, status, ierr )
              ENDIF
            ENDIF
          ENDDO
C         ------------------------------
C         In case of parallel Schur, we
C         free root%RHS_CNTR_MASTER_ROOT
C         ------------------------------
          IF (id%MYID.EQ.id_schur) THEN
            DEALLOCATE(id%root%RHS_CNTR_MASTER_ROOT)
            NULLIFY   (id%root%RHS_CNTR_MASTER_ROOT)
          ENDIF
        ENDIF
C       return because this is all we need to do
C       in case of parallel Schur complement
        RETURN
      ENDIF
C     ============================
C     Centralized Schur complement
C     ============================
C     PTRAST has been freed at the moment of calling this
C     routine. Schur is available through
C     PTRFAC(IW( PTLUST_S( STEP(KEEP(20)) ) + 4 +KEEP(IXSZ) ))
      IF (id%KEEP(252).EQ.0) THEN
C       CASE 1 (ORIGINAL CODE):
C       Schur is contiguous on ID_SCHUR
        IF ( id_schur .EQ. master ) THEN ! Necessarily equals id%MYID
C         ---------------------
C         Copy Schur complement
C         ---------------------
          CALL zmumps_copyi8size( surfschur8,
     &      id%S(id%PTRFAC(id%STEP(id%KEEP(20)))),
     &      id%SCHUR(1) )
        ELSE
C         -----------------------------------------
C         The processor responsible of the Schur
C         complement sends it to the host processor
C         Use blocks to avoid too large messages.
C         -----------------------------------------
          bl8=int(huge(i4)/id%KEEP(35)/10,8)
          DO ib=1, int((surfschur8+bl8-1_8) / bl8)
            shift8 = int(ib-1,8) * bl8                ! Where to send
            bl4    = int(min(bl8,surfschur8-shift8)) ! Size of block
            IF ( id%MYID .eq. id_schur ) THEN
C             Send Schur complement
              CALL mpi_send( id%S( shift8 +
     &          id%PTRFAC(id%IS(id%PTLUST_S(id%STEP(id%KEEP(20)))
     &                    +4+id%KEEP(ixsz)))),
     &          bl4,
     &          mpi_double_complex,
     &          master, tag_schur,
     &          id%COMM, ierr )
            ELSE IF ( id%MYID .eq. master ) THEN
C             Receive Schur complement
              CALL mpi_recv( id%SCHUR(1_8 + shift8),
     &                     bl4,
     &                     mpi_double_complex, id_schur, tag_schur,
     &                     id%COMM, status, ierr )
            END IF
          ENDDO
        END IF
      ELSE
C       CASE 2 (Fwd in facto): Schur is not contiguous on ID_SCHUR,
C       process it row by row.
C
C       2.1: We first centralize Schur complement into id%SCHUR
        ischur_src = id%PTRFAC(id%IS(id%PTLUST_S(id%STEP(id%KEEP(20)))
     &               +4+id%KEEP(ixsz)))
        ischur_dest= 1_8
        DO i=1, size_schur
          row_length = size_schur
          IF (id_schur.EQ.master) THEN ! Necessarily = id%MYID
            CALL zcopy(row_length, id%S(ischur_src), 1,
     &                 id%SCHUR(ischur_dest),1)
          ELSE
            IF (id%MYID.EQ.id_schur) THEN
C             Send
              CALL mpi_send( id%S(ischur_src), row_length,
     &        mpi_double_complex,
     &        master, tag_schur,
     &        id%COMM, ierr )
            ELSE
C             Recv
              CALL mpi_recv( id%SCHUR(ischur_dest),
     &                   row_length,
     &                   mpi_double_complex, id_schur, tag_schur,
     &                   id%COMM, status, ierr )
            ENDIF
          ENDIF
          ischur_src = ischur_src+int(ld_schur,8)
          ischur_dest= ischur_dest+int(size_schur,8)
        ENDDO
C       2.2: Get REDRHS on host
C       2.2.1: Symmetric => REDRHS is available in last KEEP(253)
C              rows of Schur structure on ID_SCHUR
C       2.2.2: Unsymmetric => REDRHS corresponds to last KEEP(253)
C              columns. However it must be transposed.
        IF (id%KEEP(221).EQ.1) THEN ! Implies Fwd in facto
          ischur_sym = id%PTRFAC(id%IS(id%PTLUST_S(id%STEP(id%KEEP(20)))
     &                    +4+id%KEEP(ixsz))) + int(size_schur,8) *
     &                    int(ld_schur,8)
          ischur_uns =
     &                 id%PTRFAC(id%IS(id%PTLUST_S(id%STEP(id%KEEP(20)))
     &                    +4+id%KEEP(ixsz))) + int(size_schur,8)
          ischur_dest = 1_8
          DO i = 1, id%KEEP(253)
            IF (id_schur .EQ. master) THEN ! necessarily = id%MYID
              IF (id%KEEP(50) .EQ. 0) THEN
                CALL zcopy(size_schur, id%S(ischur_uns), ld_schur,
     &                     id%REDRHS(ischur_dest), 1)
              ELSE
                CALL zcopy(size_schur, id%S(ischur_sym), 1,
     &                     id%REDRHS(ischur_dest), 1)
              ENDIF
            ELSE
              IF (id%MYID .NE. master) THEN
                IF (id%KEEP(50) .EQ. 0) THEN
C                 Use id%S(ISCHUR_SYM) as temporary contig. workspace
C                 of size SIZE_SCHUR. 
                  CALL zcopy(size_schur, id%S(ischur_uns), ld_schur,
     &            id%S(ischur_sym), 1)
                ENDIF
                CALL mpi_send(id%S(ischur_sym), size_schur,
     &          mpi_double_complex, master, tag_schur,
     &          id%COMM, ierr )
              ELSE
                CALL mpi_recv(id%REDRHS(ischur_dest),
     &          size_schur, mpi_double_complex, id_schur, tag_schur,
     &          id%COMM, status, ierr )
              ENDIF
            ENDIF
            IF (id%KEEP(50).EQ.0) THEN
              ischur_uns = ischur_uns + int(ld_schur,8)
            ELSE
              ischur_sym = ischur_sym + int(ld_schur,8)
            ENDIF
            ischur_dest = ischur_dest + int(id%LREDRHS,8)
          ENDDO
        ENDIF
      ENDIF
      RETURN

zmumps_fac_driver()

subroutine zmumps_fac_driver

(

type(zmumps_struc), target

id

)

Definition at line 14 of file zfac_driver.F.

      USE zmumps_buf
      USE zmumps_load
      USE zmumps_ooc
      USE zmumps_struc_def
      USE zmumps_facsol_l0omp_m, ONLY: zmumps_free_l0_omp_factors,
     &  zmumps_init_l0_omp_factors
      USE zmumps_lr_stats
      USE zmumps_lr_data_m, only: zmumps_blr_init_module, 
     &                            zmumps_blr_end_module
     &                          , zmumps_blr_mod_to_struc
      USE mumps_front_data_mgt_m
#if ! defined(NO_FDM_DESCBAND)
      USE mumps_fac_descband_data_m
#endif
#if ! defined(NO_FDM_MAPROW)
      USE mumps_fac_maprow_data_m
#endif
!$    USE OMP_LIB
C     Derived datatype to pass pointers with implicit interfaces
      USE zmumps_fac_s_is_pointers_m, ONLY : s_is_pointers_t
      IMPLICIT NONE
C
C  Purpose
C  =======
C
C  Performs scaling, sorting in arrowhead, then
C  distributes the matrix, and perform
C  factorization.
C
C
      INTERFACE
      SUBROUTINE zmumps_anorminf(id, ANORMINF, LSCAL, EFF_SIZE_SCHUR)
      USE zmumps_struc_def
      TYPE (ZMUMPS_STRUC), TARGET :: id
      DOUBLE PRECISION, INTENT(OUT) :: ANORMINF
      LOGICAL, INTENT(IN) :: LSCAL
      INTEGER, INTENT(IN) :: EFF_SIZE_SCHUR
      END SUBROUTINE zmumps_anorminf
      SUBROUTINE zmumps_free_id_data_modules(id_FDM_F_ENCODING,
     &  id_BLRARRAY_ENCODING, KEEP8, K34)
#     if defined(MUMPS_F2003)
      CHARACTER, DIMENSION(:), POINTER, intent(inout) ::
     &                                            id_blrarray_encoding
      CHARACTER, DIMENSION(:), POINTER, intent(inout) ::
     &                                            id_fdm_f_encoding
#     else
      CHARACTER, DIMENSION(:), POINTER :: id_BLRARRAY_ENCODING
      CHARACTER, DIMENSION(:), POINTER :: id_FDM_F_ENCODING
#     endif
      INTEGER(8), intent(inout) :: KEEP8(150)
      INTEGER, intent(in) :: K34
      END SUBROUTINE zmumps_free_id_data_modules
      END INTERFACE
C
C  Parameters
C  ==========
C
      TYPE(ZMUMPS_STRUC), TARGET :: id
C
C  MPI
C  ===
C
      include 'mpif.h'
      include 'mumps_tags.h'
      INTEGER :: STATUS(MPI_STATUS_SIZE)
      INTEGER :: IERR
      INTEGER, PARAMETER :: MASTER = 0
C
C  Local variables
C  ===============
C
      include 'mumps_headers.h'
      INTEGER(8) :: NSEND8, NSEND_TOT8
      INTEGER(8) :: NLOCAL8, NLOCAL_TOT8
      INTEGER(4) :: I4
      INTEGER :: LDPTRAR, NELT_arg, NBRECORDS
      INTEGER :: ITMP, JTMP
      INTEGER :: KEEP464COPY, KEEP465COPY
      DOUBLE PRECISION    :: RATIOK465
      INTEGER(8) :: KEEP826_SAVE
      INTEGER(8) :: K67, K68, K70, K74, K75
      INTEGER(8) ITMP8
      INTEGER  MUMPS_PROCNODE
      EXTERNAL mumps_procnode
      INTEGER MP, LP, MPG, allocok
      LOGICAL PROK, PROKG, LSCAL, LPOK, COMPUTE_ANORMINF
C     Reception buffer
      INTEGER    :: ZMUMPS_LBUFR, ZMUMPS_LBUFR_BYTES
      INTEGER(8) :: ZMUMPS_LBUFR_BYTES8 ! for intermediate computation
      INTEGER, ALLOCATABLE, DIMENSION(:) :: BUFR
C     Size of send buffers (in bytes)
      INTEGER    :: ZMUMPS_LBUF, ZMUMPS_LBUF_INT
      INTEGER(8) :: ZMUMPS_LBUF8 ! for intermediate computation
C
      INTEGER PTRIST, PTRWB, MAXELT_SIZE,
     &     itloc, ipool, k28, lpool
      INTEGER IRANK, ID_ROOT
      INTEGER KKKK
      INTEGER(8) :: NZ_locMAX8
      INTEGER(8) MEMORY_MD_ARG
      INTEGER(8) MAXS_BASE8, MAXS_BASE_RELAXED8
      DOUBLE PRECISION CNTL4, AVG_FLOPS
      INTEGER MIN_PERLU, MAXIS_ESTIM
      INTEGER SUM_INFO22_THIS_NODE, MAX_SUM_INFO22_THIS_NODE
C
      TYPE (S_IS_POINTERS_T) :: S_IS_POINTERS
      INTEGER   MAXIS
      INTEGER(8) :: MAXS
C     For S argument to arrowhead routines:
      INTEGER(8) :: MAXS_ARG
      COMPLEX(kind=8), TARGET :: S_DUMMY_ARG(1)
      COMPLEX(kind=8), POINTER, DIMENSION(:) :: S_PTR_ARG
      INTEGER NB_THREADS, NOMP
      DOUBLE PRECISION TIMEAVG, TIMEMAX,
     &                 flopavg, flopmax
      DOUBLE PRECISION TMPTIME, TMPFLOP
      INTEGER NPIV_CRITICAL_PATH, EFF_SIZE_SCHUR
      DOUBLE PRECISION TIME, TIMEET
      DOUBLE PRECISION ZERO, ONE, MONE
      parameter( zero = 0.0d0, one = 1.0d0, mone = -1.0d0)
      COMPLEX(kind=8) CZERO
      parameter( czero = (0.0d0, 0.0d0) )
      INTEGER PERLU, TOTAL_MBYTES, K231, K232, K233, BLR_STRAT
      INTEGER, PARAMETER :: IDUMMY = -9999
      LOGICAL, PARAMETER :: BDUMMY =.false.
      INTEGER, PARAMETER :: PANEL_TABSIZE = 20
      INTEGER COLOUR, COMM_FOR_SCALING ! For Simultaneous scaling
      INTEGER LIWK, LWK_REAL
      INTEGER(8) :: LWK
C     I_AM_SLAVE: used to determine if proc has the role of a slave
C     WK_USER_PROVIDED is set to true when WK_USER is provided by user
      LOGICAL I_AM_SLAVE, PERLU_ON, WK_USER_PROVIDED, EARLYT3ROOTINS
      LOGICAL PRINT_MAXAVG, PRINT_NODEINFO
      DOUBLE PRECISION :: ANORMINF, SEUIL, SEUIL_LDLT_NIV2, Thresh_Seuil
      DOUBLE PRECISION :: CNTL1, CNTL3, CNTL5, CNTL6, EPS
      INTEGER N, LPN_LIST,POSBUF
      INTEGER, DIMENSION (:), ALLOCATABLE :: ITMP2
      INTEGER I,K
      INTEGER(8) :: ITEMP8
      INTEGER    :: PARPIV_T1
      INTEGER FRONTWISE
C temporary variables for collecting stats from all processors
      DOUBLE PRECISION :: TMP_MRY_LU_FR
      DOUBLE PRECISION :: TMP_MRY_LU_LRGAIN
      DOUBLE PRECISION :: TMP_MRY_CB_FR
      DOUBLE PRECISION :: TMP_MRY_CB_LRGAIN
      DOUBLE PRECISION :: TMP_FLOP_LRGAIN
      DOUBLE PRECISION :: TMP_FLOP_TRSM
      DOUBLE PRECISION :: TMP_FLOP_PANEL
      DOUBLE PRECISION :: TMP_FLOP_FRFRONTS
      DOUBLE PRECISION :: TMP_FLOP_TRSM_FR
      DOUBLE PRECISION :: TMP_FLOP_TRSM_LR
      DOUBLE PRECISION :: TMP_FLOP_UPDATE_FR
      DOUBLE PRECISION :: TMP_FLOP_UPDATE_LR
      DOUBLE PRECISION :: TMP_FLOP_UPDATE_LRLR3
      DOUBLE PRECISION :: TMP_FLOP_COMPRESS
      DOUBLE PRECISION :: TMP_FLOP_DECOMPRESS
      DOUBLE PRECISION :: TMP_FLOP_MIDBLK_COMPRESS
      DOUBLE PRECISION :: TMP_FLOP_FRSWAP_COMPRESS
      DOUBLE PRECISION :: TMP_FLOP_ACCUM_COMPRESS
      DOUBLE PRECISION :: TMP_FLOP_CB_COMPRESS
      DOUBLE PRECISION :: TMP_FLOP_CB_DECOMPRESS
      DOUBLE PRECISION :: TMP_FLOP_FACTO_FR
      INTEGER :: TMP_CNT_NODES
      DOUBLE PRECISION :: TMP_TIME_UPDATE
      DOUBLE PRECISION :: TMP_TIME_UPDATE_LRLR1
      DOUBLE PRECISION :: TMP_TIME_UPDATE_LRLR2
      DOUBLE PRECISION :: TMP_TIME_UPDATE_LRLR3
      DOUBLE PRECISION :: TMP_TIME_UPDATE_FRLR
      DOUBLE PRECISION :: TMP_TIME_UPDATE_FRFR
      DOUBLE PRECISION :: TMP_TIME_COMPRESS
      DOUBLE PRECISION :: TMP_TIME_MIDBLK_COMPRESS
      DOUBLE PRECISION :: TMP_TIME_FRSWAP_COMPRESS
      DOUBLE PRECISION :: TMP_TIME_CB_COMPRESS
      DOUBLE PRECISION :: TMP_TIME_PANEL
      DOUBLE PRECISION :: TMP_TIME_FAC_I
      DOUBLE PRECISION :: TMP_TIME_FAC_MQ
      DOUBLE PRECISION :: TMP_TIME_FAC_SQ
      DOUBLE PRECISION :: TMP_TIME_LRTRSM
      DOUBLE PRECISION :: TMP_TIME_FRTRSM
      DOUBLE PRECISION :: TMP_TIME_FRFRONTS
      DOUBLE PRECISION :: TMP_TIME_LR_MODULE
      DOUBLE PRECISION :: TMP_TIME_DIAGCOPY
      DOUBLE PRECISION :: TMP_TIME_DECOMP
      DOUBLE PRECISION :: TMP_TIME_DECOMP_UCFS
      DOUBLE PRECISION :: TMP_TIME_DECOMP_ASM1
      DOUBLE PRECISION :: TMP_TIME_DECOMP_LOCASM2
      DOUBLE PRECISION :: TMP_TIME_DECOMP_MAPLIG1
      DOUBLE PRECISION :: TMP_TIME_DECOMP_ASMS2S
      DOUBLE PRECISION :: TMP_TIME_DECOMP_ASMS2M
C
C  Workspace.
C
      INTEGER, DIMENSION(:), ALLOCATABLE :: IWK
      COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE :: WK
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: WK_REAL
      INTEGER(8), DIMENSION(:), ALLOCATABLE :: IWK8
      INTEGER, DIMENSION(:), ALLOCATABLE :: BURP
      INTEGER, DIMENSION(:), ALLOCATABLE :: BUCP
      INTEGER, DIMENSION(:), ALLOCATABLE :: BURS
      INTEGER, DIMENSION(:), ALLOCATABLE :: BUCS
      INTEGER BUREGISTRE(12)
      INTEGER BUINTSZ, BURESZ, BUJOB
      INTEGER BUMAXMN, M, SCMYID, SCNPROCS
      DOUBLE PRECISION    SCONEERR, SCINFERR
C
C  Parameters arising from the structure
C  =====================================
C
      INTEGER, POINTER ::  JOB
*     Control parameters: see description in ZMUMPSID
      DOUBLE PRECISION,DIMENSION(:),POINTER::RINFO, RINFOG
      DOUBLE PRECISION,DIMENSION(:),POINTER::    CNTL
      INTEGER,DIMENSION(:),POINTER:: INFOG, KEEP
      INTEGER, DIMENSION(:), POINTER :: MYIRN_loc, MYJCN_loc
      COMPLEX(kind=8), DIMENSION(:), POINTER :: MYA_loc
      INTEGER, TARGET :: DUMMYIRN_loc(1), DUMMYJCN_loc(1)
      COMPLEX(kind=8), TARGET :: DUMMYA_loc(1)
      INTEGER,DIMENSION(:),POINTER::ICNTL
      EXTERNAL mumps_get_pool_length
      INTEGER MUMPS_GET_POOL_LENGTH
      INTEGER(8) :: TOTAL_BYTES
      INTEGER(8) :: I8TMP, LWK_USER_SUM8
C
C  External references
C  ===================
      INTEGER numroc
      EXTERNAL numroc
      INTEGER:: NWORKING
      LOGICAL:: MEM_EFF_ALLOCATED
      INTEGER   :: TOTAL_MBYTES_UNDER_L0
      INTEGER(8):: TOTAL_BYTES_UNDER_L0
C  Fwd in facto:
      COMPLEX(kind=8), DIMENSION(:), POINTER :: RHS_MUMPS
      LOGICAL :: RHS_MUMPS_ALLOCATED
      INTEGER :: NB_ACTIVE_FRONTS_ESTIM
      INTEGER :: NB_FRONTS_F_ESTIM
C
C 
      job=>id%JOB
      rinfo=>id%RINFO
      rinfog=>id%RINFOG
      cntl=>id%CNTL
      infog=>id%INFOG
      keep=>id%KEEP
      icntl=>id%ICNTL
      IF (id%KEEP8(29) .NE. 0) THEN
        myirn_loc=>id%IRN_loc
        myjcn_loc=>id%JCN_loc
        mya_loc=>id%A_loc
      ELSE
        myirn_loc=>dummyirn_loc
        myjcn_loc=>dummyjcn_loc
        mya_loc=>dummya_loc
      ENDIF
      n = id%N
      eps = epsilon( zero )
C     TIMINGS: reset to 0
      id%DKEEP(92)=0.0d0
      id%DKEEP(93)=0.0d0
      id%DKEEP(94)=0.0d0
      id%DKEEP(97)=0.0d0
      id%DKEEP(98)=0.0d0
      id%DKEEP(56)=0.0d0
C     Count of MPI messages: reset to 0
      id%KEEP(266)=0
      id%KEEP(267)=0
C     MIN/MAX pivots reset to 0
      id%DKEEP(19)=huge(0.0d0)
      id%DKEEP(20)=huge(0.0d0)
      id%DKEEP(21)=0.0d0
C     Number of symmetric swaps
      id%KEEP8(80)=0_8
C     Largest increase of internal panel size
      id%KEEP(425) =0
C
      print_maxavg = .NOT.(id%NSLAVES.EQ.1 .AND. keep(46).EQ.1)
C     Print per node informtation only in case ther are several
C     compute nodes (id%KEEP(412): #MPI procs on comupte node)
      print_nodeinfo = print_maxavg .AND. id%NPROCS .NE. id%KEEP(412)
C
C     BEGIN CASE OF ALLOCATED DATA FROM PREVIOUS CALLS
C     Data from factorization is now freed asap
C     id%S, id%IS
      IF (id%KEEP8(24).EQ.0_8) THEN
C       -- deallocate only when not provided/allocated by the user
         IF (associated(id%S)) THEN
            DEALLOCATE(id%S)
            id%KEEP8(23)=0_8
            NULLIFY(id%S)
         ENDIF
      ENDIF
      IF (associated(id%IS)) THEN
        DEALLOCATE(id%IS)
        NULLIFY(id%IS)
      ENDIF
C     Free BLR factors, if any
      CALL zmumps_free_id_data_modules(id%FDM_F_ENCODING,
     &            id%BLRARRAY_ENCODING, id%KEEP8(1), id%KEEP(34))
      IF (associated(id%root%RG2L_ROW))THEN
        DEALLOCATE(id%root%RG2L_ROW)
        NULLIFY(id%root%RG2L_ROW)
      ENDIF
      IF (associated(id%root%RG2L_COL))THEN
        DEALLOCATE(id%root%RG2L_COL)
        NULLIFY(id%root%RG2L_COL)
      ENDIF
      IF (associated( id%PTLUST_S )) THEN
        DEALLOCATE(id%PTLUST_S)
        NULLIFY(id%PTLUST_S)
      ENDIF
      IF (associated(id%PTRFAC)) THEN
        DEALLOCATE(id%PTRFAC)
        NULLIFY(id%PTRFAC)
      END IF
      IF (associated(id%RHSCOMP)) THEN
        DEALLOCATE(id%RHSCOMP)
        NULLIFY(id%RHSCOMP)
        id%KEEP8(25)=0_8
      ENDIF
      IF (associated(id%POSINRHSCOMP_ROW)) THEN
        DEALLOCATE(id%POSINRHSCOMP_ROW)
        NULLIFY(id%POSINRHSCOMP_ROW)
      ENDIF
      IF (id%POSINRHSCOMP_COL_ALLOC) THEN
        DEALLOCATE(id%POSINRHSCOMP_COL)
        NULLIFY(id%POSINRHSCOMP_COL)
        id%POSINRHSCOMP_COL_ALLOC = .false.
      ENDIF
C
C     END CASE OF ALLOCATED DATA FROM PREVIOUS CALLS
C
C     Related to forward in facto functionality (referred to as "Fwd in facto")
      NULLIFY(rhs_mumps)
      rhs_mumps_allocated = .false.
C     -----------------------------------------------------------------------
C     Set WK_USER_PROVIDED to true when workspace WK_USER is provided by user
C     We can accept WK_USER to be provided on only one proc and 
C     different values of WK_USER per processor
C     
      IF (id%KEEP8(24).GT.0_8) THEN 
C                We nullify S so that later when we test
C                if (associated(S) we can free space and reallocate it).
           NULLIFY(id%S)
      ENDIF
C
C     --  KEEP8(24) can now then be reset safely
      wk_user_provided = (id%LWK_USER.NE.0)
      IF (wk_user_provided) THEN
          IF (id%LWK_USER.GT.0) THEN
            id%KEEP8(24) = int(id%LWK_USER,8)
          ELSE
            id%KEEP8(24) = -int(id%LWK_USER,8)* 1000000_8 
          ENDIF
      ELSE
          id%KEEP8(24) = 0_8
      ENDIF
C     Compute sum of LWK_USER provided by user 
      lwk_user_sum8 = 0_8
      CALL mpi_reduce ( id%KEEP8(24), lwk_user_sum8, 1, mpi_integer8,
     &                  mpi_sum, master, id%COMM, ierr )
C
C     KEEP8(26) might be modified
C       (element entry format) 
C       but need be restore for 
C       future factorisation
C       with different scaling option
C 
      keep826_save = id%KEEP8(26)
C     In case of loop on factorization with
C     different scaling options, initialize
C     DKEEP(4:5) to 0.
      id%DKEEP(4)=-1.0d0
      id%DKEEP(5)=-1.0d0
C  Mapping information used during solve. In case of several facto+solve
C  it has to be recomputed. In case of several solves with the same
C  facto, it is not recomputed.
      IF (associated(id%IPTR_WORKING)) THEN
        DEALLOCATE(id%IPTR_WORKING)
        NULLIFY(id%IPTR_WORKING)
      END IF
      IF (associated(id%WORKING)) THEN 
        DEALLOCATE(id%WORKING)
        NULLIFY(id%WORKING)
      END IF
C
C  Units for printing
C  MP: diagnostics
C  LP: errors
C
      lp  = icntl( 1 )
      mp  = icntl( 2 )
      mpg = icntl( 3 )
      lpok    = ((lp.GT.0).AND.(id%ICNTL(4).GE.1))
      prok    = ((mp.GT.0).AND.(id%ICNTL(4).GE.2))
      prokg   = ( mpg .GT. 0 .and. id%MYID .eq. master )
      prokg   = (prokg.AND.(id%ICNTL(4).GE.2))
      IF ( prok ) WRITE( mp, 130 )
      IF ( prokg ) WRITE( mpg, 130 )
C     -------------------------------------
C     Depending on the type of parallelism,
C     the master can now (soon) potentially
C     have the role of a slave
C     -------------------------------------
      i_am_slave = ( id%MYID .ne. master  .OR.
     &             ( id%MYID .eq. master .AND.
     &               keep(46) .eq. 1 ) )
C
C  Prepare work for out-of-core
C
      IF (id%MYID .EQ. master .AND. keep(201) .NE. -1) THEN
C       Note that if KEEP(201)=-1, then we have decided
C       at analysis phase that factors will not be stored
C       (neither in memory nor on disk). In that case,
C       ICNTL(22) is ignored.
C       -- ICNTL(22) must be set before facto phase 
C          (=1 OOC on; =0 OOC off)
C          and cannot be changed for subsequent solve phases.
        keep(201)=id%ICNTL(22)
        IF (keep(201) .NE. 0) THEN
#         if defined(OLD_OOC_NOPANEL)
            keep(201)=2
#         else
            keep(201)=1
#         endif
        ENDIF
      ENDIF
C     ----------------------
C     Broadcast KEEP options
C     defined for facto:
C     ----------------------
      CALL mpi_bcast( keep(12), 1, mpi_integer,
     &                master, id%COMM, ierr )
      CALL mpi_bcast( keep(19), 1, mpi_integer,
     &                master, id%COMM, ierr )
      CALL mpi_bcast( keep(21), 1, mpi_integer,
     &                master, id%COMM, ierr )
      CALL mpi_bcast( keep(201), 1, mpi_integer,
     &                master, id%COMM, ierr )
      CALL mpi_bcast( keep(459), 1, mpi_integer,
     &                master, id%COMM, ierr )
      CALL mpi_bcast( keep(460), 1, mpi_integer,
     &                master, id%COMM, ierr )
      IF ( keep(459) .GE. panel_tabsize ) THEN
        IF ( lpok ) THEN
        WRITE(lp,'(A,I4,A,I3)') " ** WARNING ** KEEP(459)=",keep(459),
     &  " too large, resetting to",panel_tabsize-1
        ENDIF
        keep(459) = panel_tabsize - 1
      ENDIF
      perlu = keep(12)
      IF (id%MYID.EQ.master) THEN
C       KEEP(50)  case
C       ==============
C
C       KEEP(50)  = 0 : matrix is unsymmetric
C       KEEP(50) /= 0 : matrix is symmetric
C       KEEP(50) = 1 : Ask L L^T on the root. Matrix is PSD.
C       KEEP(50) = 2 : Ask for L U on the root
C       KEEP(50) = 3 ... L D L^T ??
C     
        cntl1 = id%CNTL(1)
C       ---------------------------------------
C       For symmetric (non general) matrices 
C       set (directly) CNTL1 = 0.0
C       ---------------------------------------
        keep(17)=0
        IF ( keep(50) .eq. 1 ) THEN
          IF (cntl1 .ne. zero ) THEN
            IF ( prokg ) THEN
              WRITE(mpg,'(A)')
     & '** Warning : SPD solver called, resetting CNTL(1) to 0.0D0'
            END IF
          END IF
          cntl1 = zero
        END IF
C       CNTL1 threshold value must be between 
C       0.0 and 1.0 (for SYM=0) and 0.5 (for SYM=1,2)
        IF (cntl1.GT.one)   cntl1=one
        IF (cntl1.LT.zero)  cntl1=zero
        IF (keep(50).NE.0.AND.cntl1.GT.0.5d0) THEN
          cntl1 = 0.5d0
        ENDIF
        parpiv_t1 = id%KEEP(268)
        IF (parpiv_t1.EQ.77) THEN
         parpiv_t1 = 0 
#if defined(__ve__)
         parpiv_t1 = -2
#endif
        ENDIF
        IF (parpiv_t1.EQ.-3) THEN
          parpiv_t1 = 0
        ENDIF
        IF ((parpiv_t1.LT.-3).OR.(parpiv_t1.GT.1)) THEN
C        out of range values
         parpiv_t1 =0
        ENDIF
C       note that KEEP(50).EQ.1 => CNTL1=0.0
        IF (cntl1.EQ.0.0.OR.(keep(50).eq.1)) parpiv_t1 = 0 
C
        IF (parpiv_t1.EQ.-2) THEN
         IF (keep(19).NE.0) THEN
C         switch off PARPIV_T1 if RR activated
C         but do NOT switch off PARPIV_1 with null pivot detection
          parpiv_t1 = 0
         ENDIF
        ENDIF
        id%KEEP(269) = parpiv_t1
      ENDIF
      CALL mpi_bcast(cntl1, 1, mpi_double_precision,
     &             master, id%COMM, ierr)
        CALL mpi_bcast( keep(269), 1, mpi_integer,
     &                  master, id%COMM, ierr )
        IF (id%MYID.EQ.master) THEN
C         -----------------------------------------------------
C         Decoding of ICNTL(35) for factorization: same as
C         at analysis except that we store a copy of ICNTL(35)
C         in KEEP(486) instead of KEEP(494) and need to check
C         compatibility of KEEP(486) and KEEP(494): If LR was
C         not activated during analysis, it cannot be activated 
C         at factorization.
C         ------------------------------------------------------
          id%KEEP(486) = id%ICNTL(35)
          IF (id%KEEP(486).EQ.1) THEN
C         -- Automatic BLR option setting
           id%KEEP(486)= 2
          ENDIF
          IF ( id%KEEP(486).EQ.4) id%KEEP(486)=0
          IF ((id%KEEP(486).LT.0).OR.(id%KEEP(486).GT.4)) THEN
C           Out of range values treated as 0
            id%KEEP(486) = 0
          ENDIF
          IF ((keep(486).NE.0).AND.(keep(494).EQ.0)) THEN
C           To activate BLR during factorization,
C           ICNTL(35) must have been set at analysis.
            IF (lpok)  THEN
              WRITE(lp,'(A)') 
     &      " *** Error with BLR setting "
              WRITE(lp,'(A)') " *** BLR was not activated during ",
     &      " analysis but is requested during factorization."
            ENDIF
            id%INFO(1)=-54
            id%INFO(2)=0
            GOTO 105
          ENDIF
          keep464copy  = id%ICNTL(38)
          IF (keep464copy.LT.0.OR.keep464copy.GT.1000) THEN
C          Out of range values treated as 1000
           keep464copy = 1000
          ENDIF
          IF (id%KEEP(461).LT.1) THEN
            id%KEEP(461) = 10
          ENDIF
          keep465copy=0
          IF (id%ICNTL(36).EQ.1.OR.id%ICNTL(36).EQ.3) THEN
            IF (cntl1.EQ.zero .OR. keep(468).LE.1) THEN
              keep(475) = 3
            ELSE IF ( (keep(269).GT.0).OR. (keep(269).EQ.-2)) THEN
              keep(475) = 2
            ELSE IF (keep(468).EQ.2) THEN
              keep(475) = 2
            ELSE
              keep(475) = 1
            ENDIF
          ELSE
            keep(475) = 0
          ENDIF
          keep(481)=0
          IF (id%ICNTL(36).LT.0 .OR. id%ICNTL(36).GE.2) THEN
C           Only options 1 and 2 are allowed
            keep(475) = 0
          ENDIF
C         K489 is set according to ICNTL(37)
          IF (id%ICNTL(37).EQ.0.OR.id%ICNTL(37).EQ.1) THEN
            keep(489) = id%ICNTL(37)
          ELSE
C           Other values treated as zero
            keep(489) = 0
          ENDIF
          IF (keep(79).GE.1) THEN
C          CompressCB incompatible with type4,5,6 nodes
           keep(489)=0
          ENDIF
          keep(489)=0
C         id%KEEP(476) \in [1,100] 
          IF ((id%KEEP(476).GT.100).OR.(id%KEEP(476).LT.1)) THEN
              id%KEEP(476)=  50
          ENDIF
C         id%KEEP(477) \in [1,100] 
          IF ((id%KEEP(477).GT.100).OR.(id%KEEP(477).LT.1)) THEN
              id%KEEP(477)=  100
          ENDIF
C         id%KEEP(483) \in [1,100] 
          IF ((id%KEEP(483).GT.100).OR.(id%KEEP(483).LT.1)) THEN
              id%KEEP(483)= 50
          ENDIF
C         id%KEEP(484) \in [1,100] 
          IF ((id%KEEP(484).GT.100).OR.(id%KEEP(484).LT.1)) THEN
              id%KEEP(484)=  50
          ENDIF
C         id%KEEP(480)=0,2,3,4,5,6
          IF ((id%KEEP(480).GT.6).OR.(id%KEEP(480).LT.0)
     &                           .OR.(id%KEEP(480).EQ.1)) THEN
               id%KEEP(480)=0
          ENDIF
C         id%KEEP(473)=0 or 1
          IF ((id%KEEP(473).NE.0).AND.(id%KEEP(473).NE.1)) THEN
               id%KEEP(473)=0
          ENDIF
C         id%KEEP(474)=0,1,2,3
          IF ((id%KEEP(474).GT.3).OR.(id%KEEP(474).LT.0)) THEN
               id%KEEP(474)=0
          ENDIF
C         id%KEEP(479)>0
          IF (id%KEEP(479).LE.0) THEN
               id%KEEP(479)=1
          ENDIF
         IF (id%KEEP(474).NE.0.AND.id%KEEP(480).EQ.0) THEN
              id%KEEP(474) = 0
          ENDIF
          IF (id%KEEP(478).NE.0.AND.id%KEEP(480).LT.4) THEN
              id%KEEP(478) = 0
          ENDIF
          IF (id%KEEP(480).GE.5 .OR.
     &           (id%KEEP(480).NE.0.AND.id%KEEP(474).EQ.3)) THEN
            IF (id%KEEP(475).LT.2) THEN
C             Reset to 3 if 5 or to 4 if 6
              id%KEEP(480) = id%KEEP(480) - 2
              write(*,*) ' Resetting KEEP(480) to ', id%KEEP(480)
            ENDIF
          ENDIF
 105    CONTINUE
        ENDIF  ! id%MYID .EQ. MASTER
        CALL mumps_propinfo( id%ICNTL(1), id%INFO(1),
     &                        id%COMM, id%MYID )
C
        IF (id%INFO(1).LT.0) GOTO 530
        CALL mpi_bcast( keep(473), 14, mpi_integer,
     &                  master, id%COMM, ierr )
        IF (keep(486).NE.0) THEN
          CALL mpi_bcast( keep(489), 1, mpi_integer,
     &                  master, id%COMM, ierr )
          CALL mpi_bcast( keep464copy, 1, mpi_integer,
     &                  master, id%COMM, ierr )
          CALL mpi_bcast( keep465copy, 1, mpi_integer,
     &                  master, id%COMM, ierr )
         ENDIF
        IF (id%MYID.EQ.master) THEN
          IF (keep(217).GT.2.OR.keep(217).LT.0) THEN
            keep(217)=0
          ENDIF
          keep(214)=keep(217)
          IF (keep(214).EQ.0) THEN
            IF (keep(201).NE.0) THEN ! OOC or no factors
              keep(214)=1
            ELSE
              keep(214)=2
            ENDIF
            IF (keep(486).EQ.2) THEN
              keep(214)=1
            ENDIF
          ENDIF
        ENDIF
        CALL mpi_bcast( keep(214), 1, mpi_integer,
     &                  master, id%COMM, ierr )
        IF (keep(201).NE.0) THEN
C         -- Low Level I/O strategy
          CALL mpi_bcast( keep(99), 1, mpi_integer,
     &                  master, id%COMM, ierr )
          CALL mpi_bcast( keep(205), 1, mpi_integer,
     &                  master, id%COMM, ierr )
          CALL mpi_bcast( keep(211), 1, mpi_integer,
     &                  master, id%COMM, ierr )
        ENDIF
C     Fwd in facto: explicitly forbid
C     sparse RHS and A-1 computation
      IF (id%KEEP(252).EQ.1 .AND. id%MYID.EQ.master) THEN
        IF (id%ICNTL(20).EQ.1) THEN ! out-of-range => 0
C         NB: in doc ICNTL(20) only accessed during solve
C         In practice, will have failed earlier if RHS not allocated.
C         Still it looks safer to keep this test.
          id%INFO(1)=-43
          id%INFO(2)=20
          IF (lpok) WRITE(lp,'(A)')
     &       ' ERROR: Sparse RHS is incompatible with forward',
     &       ' performed during factorization (ICNTL(32)=1)'
        ELSE IF (id%ICNTL(30).NE.0) THEN ! out-of-range => 1
          id%INFO(1)=-43
          id%INFO(2)=30
          IF (lpok) WRITE(lp,'(A)')
     &       ' ERROR: A-1 functionality incompatible with forward',
     &       ' performed during factorization (ICNTL(32)=1)'
        ELSE IF (id%ICNTL(9) .NE. 1) THEN
          id%INFO(1)=-43
          id%INFO(2)=9
          IF (lpok) WRITE(lp,'(A)')
     &       .NE.' ERROR: Transpose system (ICNTL(9)0) not ',
     &       ' compatible with forward performed during',
     &       ' factorization (ICNTL(32)=1)'
        ENDIF
      ENDIF
      CALL mumps_propinfo( id%ICNTL(1), id%INFO(1),
     &                        id%COMM, id%MYID )
C
      IF (id%INFO(1).LT.0) GOTO 530
C
C     The memory allowed is given by ICNTL(23) in Mbytes
C     0 means that nothing is provided.
C     Save memory available, ICNTL(23) in KEEP8(4)
C
      IF ( icntl(23) .GT. 0 ) THEN
        itmp = 1
      ELSE
        itmp = 0
      ENDIF
      CALL mpi_allreduce( itmp, jtmp, 1, mpi_integer,
     &                 mpi_sum, id%COMM, ierr)
      IF ( id%MYID.EQ.master ) THEN
C       Negative values considered 0
        itmp = max(icntl(23),0)
      END IF
      CALL mpi_bcast( itmp, 1, mpi_integer,
     &                master, id%COMM, ierr )
C     JTMP: nb of procs with nonzero ICNTL(23)
C     ITMP: value of ICNTL(23) on master
      IF ( itmp .GT. 0 .AND. jtmp .EQ. 1 ) THEN
C       ICNTL(23)>0 only on master
      ELSE
C       Local values of ICNTL(23) are used, note that
C       they could all be zeros
        itmp = icntl(23)
      ENDIF
C 
      itmp8 = int(itmp, 8)
      id%KEEP8(4) = itmp8 * 1000000_8   ! convert to nb of bytes
C     Compute \sum of memories allowed
      CALL mpi_reduce( id%KEEP8(4), itmp8, 1, mpi_integer8,
     &                  mpi_sum, master, id%COMM, ierr )
      itmp8 = itmp8 / 1000000_8 ! Use to print \sum_{ICNTL(23)}
      IF ( prokg ) THEN
          nworking = id%NSLAVES
          WRITE( mpg, 172 ) nworking, id%ICNTL(22), keep(486),
     &    keep(12), 
     &    id%KEEP8(111), keep(126), keep(127), keep(28),
     &    id%KEEP8(4)/1000000_8, itmp8, lwk_user_sum8, cntl1
          IF (keep(252).GT.0) 
     &    WRITE(mpg,173) keep(253)
          IF (keep(269).NE.0) 
     &    WRITE(mpg,174) keep(269)
      ENDIF
      IF (keep(201).LE.0) THEN
C       In-core version or no factors
        keep(ixsz)=xsize_ic
      ELSE IF (keep(201).EQ.2) THEN
C       OOC version, no panels
        keep(ixsz)=xsize_ooc_nopanel
      ELSE IF (keep(201).EQ.1) THEN
C     Panel versions:
        IF (keep(50).EQ.0) THEN
          keep(ixsz)=xsize_ooc_unsym
        ELSE
          keep(ixsz)=xsize_ooc_sym
        ENDIF
      ENDIF
      IF ( keep(486) .NE. 0 ) THEN !LR is activated
C       Stats initialization for LR
        CALL init_stats_global(id)
       END IF
C
*     **********************************
*     Begin intializations regarding the
*     computation of the determinant
*     **********************************
      IF (id%MYID.EQ.master) keep(258)=icntl(33)
      CALL mpi_bcast(keep(258), 1, mpi_integer,
     &               master, id%COMM, ierr)
      IF (keep(258) .NE. 0) THEN
        keep(259) = 0      ! Initial exponent of the local determinant
        keep(260) = 1      ! Number of permutations
        id%DKEEP(6)  = 1.0d0  ! real part of the local determinant
        id%DKEEP(7)  = 0.0d0  ! imaginary part of the local determinant
      ENDIF
*     ********************************
*     End intializations regarding the
*     computation of the determinant
*     ********************************
C
*     **********************
*     Begin of Scaling phase
*     **********************
C
C     SCALING MANAGEMENT
C     * Options 1, 3, 4 centralized only
C  
C     * Options 7, 8  : also works for distributed matrix
C
C     At this point, we have the scaling arrays allocated
C     on the master. They have been allocated on the master
C     inside the main MUMPS driver.
C
      CALL mpi_bcast(keep(52), 1, mpi_integer,
     &               master, id%COMM, ierr)
      lscal = ((keep(52) .GT. 0) .AND. (keep(52) .LE. 8))
      IF (lscal) THEN
C
        IF ( id%MYID.EQ.master ) THEN
          CALL mumps_secdeb(timeet)
        ENDIF
C       -----------------------
C       Retrieve parameters for
C       simultaneous scaling
C       -----------------------
        IF (keep(52) .EQ. 7) THEN 
C       -- Cheap setting of SIMSCALING (it is the default in 4.8.4)
           k231= keep(231)
           k232= keep(232)
           k233= keep(233)
        ELSEIF (keep(52) .EQ. 8) THEN
C       -- More expensive setting of SIMSCALING (it was the default in 4.8.1,2,3)
           k231= keep(239)
           k232= keep(240)
           k233= keep(241)
        ENDIF
        CALL mpi_bcast(id%DKEEP(3),1,mpi_double_precision,master,
     &       id%COMM,ierr)
C
        IF ( ((keep(52).EQ.7).OR.(keep(52).EQ.8)) .AND. 
     &       keep(54).NE.0 ) THEN
C         ------------------------------
C         Scaling for distributed matrix
C         We need to allocate scaling
C         arrays on all processors, not
C         only the master.
C         ------------------------------
           IF ( id%MYID .NE. master ) THEN
              IF ( associated(id%COLSCA))
     &             DEALLOCATE( id%COLSCA )
              IF ( associated(id%ROWSCA))
     &             DEALLOCATE( id%ROWSCA )
            ALLOCATE( id%COLSCA(n), stat=ierr)
            IF (ierr .GT.0) THEN
               id%INFO(1)=-13
               id%INFO(2)=n
            ENDIF
            ALLOCATE( id%ROWSCA(n), stat=ierr)
            IF (ierr .GT.0) THEN
               id%INFO(1)=-13
               id%INFO(2)=n
            ENDIF
         ENDIF
         m = n
         bumaxmn=m
         IF(n > bumaxmn) bumaxmn = n
         liwk = 4*bumaxmn
         ALLOCATE (iwk(liwk),burp(m),bucp(n),
     &            burs(2* (id%NPROCS)),bucs(2* (id%NPROCS)),
     &            stat=allocok)
         IF (allocok > 0) THEN
            id%INFO(1)=-13
            id%INFO(2)=liwk+m+n+4* (id%NPROCS)
         ENDIF
C        --- Propagate enventual error
         CALL mumps_propinfo( icntl(1), id%INFO(1),
     &        id%COMM, id%MYID )
         IF (id%INFO(1).LT.0) GOTO 517
C        -- estimation of memory and construction of partvecs
         bujob = 1
C        -- LWK not used
         lwk_real   = 1
         ALLOCATE(wk_real(lwk_real),
     &            stat=allocok)
         IF (allocok > 0) THEN
            id%INFO(1)=-13
            id%INFO(2)=lwk_real
         ENDIF
C        --- Propagate enventual error
         CALL mumps_propinfo( icntl(1), id%INFO(1),
     &        id%COMM, id%MYID )
         IF (id%INFO(1).LT.0) GOTO 517
         CALL zmumps_simscaleabs(
     &        myirn_loc(1), myjcn_loc(1), mya_loc(1),
     &        id%KEEP8(29),
     &        m, n,  id%NPROCS, id%MYID, id%COMM,
     &        burp, bucp,
     &        burs, bucs, buregistre,
     &        iwk, liwk,
     &        buintsz, buresz, bujob,
     &        id%ROWSCA(1), id%COLSCA(1), wk_real, lwk_real,
     &        id%KEEP(50),
     &        k231, k232, k233, 
     &        id%DKEEP(3),
     &        sconeerr, scinferr)
         IF(liwk < buintsz) THEN
            DEALLOCATE(iwk)
            liwk = buintsz
            ALLOCATE(iwk(liwk), stat=allocok)
            IF (allocok > 0) THEN
               id%INFO(1)=-13
               id%INFO(2)=liwk
            ENDIF
         ENDIF
         lwk_real = buresz
         DEALLOCATE(wk_real)
         ALLOCATE (wk_real(lwk_real), stat=allocok)
         IF (allocok > 0) THEN
            id%INFO(1)=-13
            id%INFO(2)=lwk_real
         ENDIF
C        --- Propagate enventual error
         CALL mumps_propinfo( icntl(1), id%INFO(1),
     &        id%COMM, id%MYID )
         IF (id%INFO(1).LT.0) GOTO 517
C        -- estimation of memory and construction of partvecs
         bujob = 2
         CALL zmumps_simscaleabs(
     &        myirn_loc(1), myjcn_loc(1), mya_loc(1),
     &        id%KEEP8(29),
     &        m, n,  id%NPROCS, id%MYID, id%COMM,
     &        burp, bucp,
     &        burs, bucs, buregistre,
     &        iwk, liwk,
     &        buintsz, buresz, bujob,
     &        id%ROWSCA(1), id%COLSCA(1), wk_real, lwk_real,
     &        id%KEEP(50),
     &        k231, k232, k233, 
     &        id%DKEEP(3),
     &        sconeerr, scinferr)
         id%DKEEP(4) = sconeerr
         id%DKEEP(5) = scinferr
CXXXX 
         DEALLOCATE(iwk, wk_real,burp,bucp,burs, bucs)
        ELSE IF ( keep(54) .EQ. 0 ) THEN
C         ------------------
C         Centralized matrix
C         ------------------
          IF ((keep(52).EQ.7).OR.(keep(52).EQ.8))  THEN
C             -------------------------------
C             Create a communicator of size 1
C             -------------------------------
              IF (id%MYID.EQ.master) THEN
                colour = 0
              ELSE
                colour = mpi_undefined
              ENDIF
              CALL mpi_comm_split( id%COMM, colour, 0,
     &             comm_for_scaling, ierr )
              IF (id%MYID.EQ.master) THEN
                 m = n
                 bumaxmn=n
                 IF(n > bumaxmn) bumaxmn = n
                 liwk = 1
                 ALLOCATE(iwk(liwk),burp(1),bucp(1),
     &                burs(1),bucs(1),
     &                stat=allocok)
                 IF (allocok > 0) THEN
                    id%INFO(1)=-13
                    id%INFO(2)=liwk+1+1+1+1
                    GOTO 400
                 ENDIF
                 lwk_real = m + n  
                 ALLOCATE (wk_real(lwk_real), stat=allocok)
                 IF (allocok > 0) THEN
                    id%INFO(1)=-13
                    id%INFO(2)=lwk_real
                    GOTO 400
                 ENDIF
                 CALL mpi_comm_rank(comm_for_scaling, scmyid, ierr)
                 CALL mpi_comm_size(comm_for_scaling, scnprocs, ierr)
                 bujob = 1
                 CALL zmumps_simscaleabs(
     &                id%IRN(1), id%JCN(1), id%A(1),
     &                id%KEEP8(28),
     &                m, n,  scnprocs, scmyid, comm_for_scaling,
     &                burp, bucp,
     &                burs, bucs, buregistre,
     &                iwk, liwk,
     &                buintsz, buresz, bujob,
     &                id%ROWSCA(1), id%COLSCA(1), wk_real, lwk_real,
     &                id%KEEP(50),
     &                k231, k232, k233, 
     &                id%DKEEP(3),
     &                sconeerr, scinferr)
                 IF(lwk_real < buresz) THEN
                    id%INFO(1) = -136
                    GOTO 400
                 ENDIF
                 bujob = 2
                 CALL zmumps_simscaleabs(id%IRN(1),
     &                id%JCN(1), id%A(1),
     &                id%KEEP8(28),
     &                m, n,  scnprocs, scmyid, comm_for_scaling,
     &                burp, bucp,
     &                burs, bucs, buregistre,
     &                iwk, liwk,
     &                buintsz, buresz, bujob,
     &                id%ROWSCA(1), id%COLSCA(1), wk_real, lwk_real,
     &                id%KEEP(50),
     &                k231, k232, k233, 
     &                id%DKEEP(3),
     &                sconeerr, scinferr)
                 id%DKEEP(4) = sconeerr
                 id%DKEEP(5) = scinferr
 400             CONTINUE
                 IF (allocated(wk_real)) DEALLOCATE(wk_real)
                 IF (allocated(iwk))     DEALLOCATE(iwk)
                 IF (allocated(burp))    DEALLOCATE(burp)
                 IF (allocated(bucp))    DEALLOCATE(bucp)
                 IF (allocated(burs))    DEALLOCATE(burs)
                 IF (allocated(bucs))    DEALLOCATE(bucs)
              ENDIF
C             Centralized matrix: make DKEEP(4:5) available to all processors
              CALL mpi_bcast( id%DKEEP(4),2,mpi_double_precision,
     &                        master, id%COMM, ierr )
              IF (id%MYID.EQ.master) THEN
C               Communicator should only be
C               freed on the master process
                CALL mpi_comm_free(comm_for_scaling, ierr)
              ENDIF
              CALL mumps_propinfo(icntl(1), id%INFO(1),
     &             id%COMM, id%MYID)
              IF (id%INFO(1).LT.0) GOTO 517
          ELSE IF (id%MYID.EQ.master) THEN
C           ----------------------------------
C           Centralized scaling, options 1 to 6
C           ----------------------------------
            IF (keep(52).GT.0 .AND. keep(52).LE.6) THEN
C             ---------------------
C             Allocate temporary
C             workspace for scaling
C             ---------------------
              IF ( keep(52) .eq. 5 .or. 
     &          keep(52) .eq. 6 ) THEN
C               We have an explicit copy of the original
C               matrix in complex format which should probably
C               be avoided (but do we want to keep all
C               those old scaling options ?)
                lwk = id%KEEP8(28)
              ELSE
                lwk = 1_8
              END IF
              lwk_real = 5 * n
              ALLOCATE( wk_real( lwk_real ), stat = ierr )
              IF ( ierr .GT. 0 ) THEN
                id%INFO(1) = -13
                id%INFO(2) = lwk_real
                GOTO 137
              END IF
              ALLOCATE( wk( lwk ), stat = ierr )
              IF ( ierr .GT. 0 ) THEN
                id%INFO(1) = -13
                CALL mumps_set_ierror(lwk, id%INFO(2))
                GOTO 137
              END IF
              CALL zmumps_fac_a(n, id%KEEP8(28), keep(52), id%A(1),
     &             id%IRN(1), id%JCN(1),
     &             id%COLSCA(1), id%ROWSCA(1),
     &             wk, lwk, wk_real, lwk_real, icntl(1), id%INFO(1) )
              DEALLOCATE( wk_real )
              DEALLOCATE( wk )
            ENDIF
          ENDIF
        ENDIF ! Scaling distributed matrices or centralized
        IF (keep(125).NE.0) THEN
C         ------------------------
C           If we enable the scaling of the |A11 A12| block
C           we et to 1 the scaling corresponding to the Schur
C           complement matrix A22
C         ------------------------
            IF ((keep(60).GT.0) .and. (keep(116).GT.0)) THEN
C            Schur is active, reset Schur entries to ONE
             IF ( ((keep(52).EQ.7).OR.(keep(52).EQ.8)) .AND.
     &          keep(54).NE.0 ) THEN
C             Scaling available on all procs
              DO i=1, n
               IF (id%SYM_PERM(i).GT.id%N-keep(116)) THEN
                 id%COLSCA(i) = one
                 id%ROWSCA(i) = one
               ENDIF
              ENDDO
             ELSE IF ( id%MYID .EQ. master) THEN
C             Scaling available on master
              DO i=1, n
               IF (id%SYM_PERM(i).GT.id%N-keep(116)) THEN
                 id%COLSCA(i) = one
                 id%ROWSCA(i) = one
               ENDIF
              ENDDO
             ENDIF
            ENDIF
        ENDIF 
        IF (id%MYID.EQ.master) THEN
            CALL mumps_secfin(timeet)
            id%DKEEP(92)=timeet
C         Print inf-norm after last KEEP(233) iterations of
C         scaling option KEEP(52)=7 or 8 (SimScale)
C
          IF (prokg.AND.(keep(52).EQ.7.OR.keep(52).EQ.8) 
     &             .AND. (k233+k231+k232).GT.0) THEN
           IF (k232.GT.0) WRITE(mpg, 166) id%DKEEP(4)
          ENDIF
        ENDIF
      ENDIF ! LSCAL
C
C       scaling might also be provided by the user
        lscal = (lscal .OR. (keep(52) .EQ. -1) .OR. keep(52) .EQ. -2)
        IF (lscal .AND. keep(258).NE.0 .AND. id%MYID .EQ. master) THEN
          DO i = 1, id%N
            CALL zmumps_updatedeter_scaling(id%ROWSCA(i),
     &           id%DKEEP(6),    ! determinant
     &           keep(259))   ! exponent of the determinant
          ENDDO
          IF (keep(50) .EQ. 0) THEN ! unsymmetric
            DO i = 1, id%N
              CALL zmumps_updatedeter_scaling(id%COLSCA(i),
     &           id%DKEEP(6),    ! determinant
     &           keep(259))   ! exponent of the determinant
            ENDDO
          ELSE
C           -----------------------------------------
C           In this case COLSCA = ROWSCA
C           Since determinant was initialized to 1,
C           compute square of the current determinant
C           rather than going through COLSCA.
C           -----------------------------------------
            CALL zmumps_deter_square(id%DKEEP(6), keep(259))
          ENDIF
C         Now we should have taken the
C         inverse of the scaling vectors
          CALL zmumps_deter_scaling_inverse(id%DKEEP(6), keep(259))
        ENDIF
C
C       ********************
C       End of Scaling phase
C       At this point: either (matrix is distributed and KEEP(52)=7 or 8)
C       in which case scaling arrays are allocated on all processors,
C       or scaling arrays are only on the host processor.
C       In case of distributed matrix input, we will free the scaling
C       arrays on procs with MYID .NE. 0 after the all-to-all distribution
C       of the original matrix.
C       ********************
C
 137  CONTINUE
C     Fwd in facto: in case of repeated factorizations
C     with different Schur options we prefer to free
C     systematically this array now than waiting for
C     the root node. We rely on the fact that it is
C     allocated or not during the solve phase so if
C     it was allocated in a 1st call to facto and not
C     in a second, we don't want the solve to think
C     it was allocated in the second call.
      IF (associated(id%root%RHS_CNTR_MASTER_ROOT)) THEN
        DEALLOCATE (id%root%RHS_CNTR_MASTER_ROOT)
        NULLIFY (id%root%RHS_CNTR_MASTER_ROOT)
      ENDIF
C     Fwd in facto: check that id%NRHS has not changed
      IF ( id%MYID.EQ.master.AND. keep(252).EQ.1 .AND.
     &      id%NRHS .NE. id%KEEP(253) ) THEN
C         Error: NRHS should not have
C         changed since the analysis
          id%INFO(1)=-42
          id%INFO(2)=id%KEEP(253)
      ENDIF
C     Fwd in facto: allocate and broadcast RHS_MUMPS
C     to make it available on all processors.
      IF (id%KEEP(252) .EQ. 1) THEN
          IF ( id%MYID.NE.master ) THEN
            id%KEEP(254) = n              ! Leading dimension
            id%KEEP(255) = n*id%KEEP(253) ! Tot size
            ALLOCATE(rhs_mumps(id%KEEP(255)),stat=ierr)
            IF (ierr > 0) THEN
               id%INFO(1)=-13
               id%INFO(2)=id%KEEP(255)
               IF (lpok)
     &         WRITE(lp,*) 'ERROR while allocating RHS on a slave'
               NULLIFY(rhs_mumps)
            ENDIF
            rhs_mumps_allocated = .true.
          ELSE 
C           Case of non working master
            id%KEEP(254)=id%LRHS              ! Leading dimension
            id%KEEP(255)=id%LRHS*(id%KEEP(253)-1)+id%N ! Tot size
            rhs_mumps=>id%RHS
            rhs_mumps_allocated = .false.
            IF (lscal) THEN
C             Scale before broadcast: apply row
C             scaling (remark that we assume no
C             transpose).
              DO k=1, id%KEEP(253)
                DO i=1, n
                  rhs_mumps( id%KEEP(254) * (k-1) + i )
     &          = rhs_mumps( id%KEEP(254) * (k-1) + i )
     &          * id%ROWSCA(i)
                ENDDO
              ENDDO
            ENDIF
          ENDIF
      ELSE
          id%KEEP(255)=1
          ALLOCATE(rhs_mumps(1),stat=ierr)
          IF (ierr > 0) THEN
             id%INFO(1)=-13
             id%INFO(2)=1
             IF (lpok)
     &            WRITE(lp,*) 'ERREUR while allocating RHS on a slave'
             NULLIFY(rhs_mumps)
          ENDIF
          rhs_mumps_allocated = .true.
      ENDIF
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF ( id%INFO(1).lt.0 ) GOTO 517
      IF (keep(252) .EQ. 1) THEN
C
C         Broadcast the columns of the right-hand side
C         one by one. Leading dimension is keep(254)=N
C         on procs with MYID > 0 but may be larger on
C         the master processor.
          DO i= 1, id%KEEP(253)
            CALL mpi_bcast(rhs_mumps((i-1)*id%KEEP(254)+1), n,
     &           mpi_double_complex, master,id%COMM,ierr)
          END DO
      ENDIF
C     Keep a copy of ICNTL(24) and make it
C     available on all working processors.
      keep(110)=id%ICNTL(24)
      CALL mpi_bcast(keep(110), 1, mpi_integer,
     &               master, id%COMM, ierr)
C     KEEP(110) defaults to 0 for out of range values
      IF (keep(110).NE.1) keep(110)=0
      IF (keep(219).NE.0) THEN
       CALL zmumps_buf_max_array_minsize(max(keep(108),1),ierr)
       IF (ierr .NE. 0) THEN
C      ------------------------
C      Error allocating ZMUMPS_BUF
C      ------------------------
          id%INFO(1) = -13
          id%INFO(2) = max(keep(108),1)
       END IF
      ENDIF
C     -----------------------------------------------
C     Depending on the option used for 
C       -detecting null pivots (ICNTL(24)/KEEP(110))
C         CNTL(3) is used to set DKEEP(1)
C               ( A row is considered as null if ||row|| < DKEEP(1) )
C         CNTL(5) is then used to define if a large 
C                 value is set on the diagonal or if a 1 is set
C                 and other values in the row are reset to zeros.
C         SEUIL* corresponds to the minimum required 
C                absolute value of pivot.
C         SEUIL_LDLT_NIV2 is used only in the 
C                case of SYM=2 within a niv2 node for which 
C                we have only a partial view of the fully summed rows.
      IF (id%MYID .EQ. master) cntl3 = id%CNTL(3)
      CALL mpi_bcast(cntl3, 1, mpi_double_precision,
     &               master, id%COMM, ierr)
      IF (id%MYID .EQ. master) cntl5 = id%CNTL(5)
      CALL mpi_bcast(cntl5, 1, mpi_double_precision,
     &               master, id%COMM, ierr)
      IF (id%MYID .EQ. master) cntl6 = id%CNTL(6)
      CALL mpi_bcast(cntl6, 1, mpi_double_precision,
     &               master, id%COMM, ierr)
      IF (id%MYID .EQ. master) id%DKEEP(8) = id%CNTL(7)
      CALL mpi_bcast(id%DKEEP(8), 1, mpi_double_precision,
     &               master, id%COMM, ierr)
      id%DKEEP(11) = id%DKEEP(8)/id%KEEP(461)
      id%DKEEP(12) = id%DKEEP(8)/id%KEEP(462)
      IF (keep(486).EQ.0) id%DKEEP(8) = zero
      compute_anorminf = .false.
      IF ( (keep(486) .NE. 0).AND. (id%DKEEP(8).LT.zero)) THEN
        compute_anorminf = .true.
      ENDIF
      IF (keep(19).NE.0) THEN
C       Rank revealing factorisation
        compute_anorminf = .true.
      ENDIF
      IF (keep(110).NE.0) THEN
C       Null pivot detection       
        compute_anorminf = .true.
      ENDIF
      IF (id%DKEEP(8).LT.zero) THEN
C       Experimental setting of CNTL(7)
        IF (compute_anorminf) THEN
         eff_size_schur = 0
         CALL zmumps_anorminf(  id , anorminf, lscal, eff_size_schur )
C        If no schur ANORMINF fine for other cases
        ELSE
         anorminf = zero
        ENDIF
        id%DKEEP(8) = abs(id%DKEEP(8))*anorminf
C       ANORMINF need be recomputed in case of schur
        IF ((keep(60).GT.0).AND.keep(116).GT.0) anorminf=zero
      ENDIF
C     -------------------------------------------------------
C        We compute ANORMINF, when needed, based on
C        the infinite norm of Rowsca *A*Colsca
C        and make it available on all working processes.
      IF (compute_anorminf) THEN
         eff_size_schur = 0
         IF (keep(60).GT.0) eff_size_schur = keep(116)
         CALL zmumps_anorminf(  id , anorminf, lscal, eff_size_schur )
      ELSE
         anorminf = zero
      ENDIF
C
       IF ((keep(19).NE.0).OR.(keep(110).NE.0)) THEN 
         IF (prokg) THEN 
           IF (keep(19).NE.0) THEN
           WRITE(mpg,'(A,1PD16.4)')
     &    ' CNTL(3) for null pivot rows/singularities  =',cntl3      
           ELSE
           WRITE(mpg,'(A,1PD16.4)')
     &    ' CNTL(3) for null pivot row detection       =',cntl3      
           ENDIF
         ENDIF
       ENDIF
       IF (keep(19).EQ.0) THEN 
C        -- RR is off
         seuil = zero
         id%DKEEP(9) = zero
       ELSE
C        -- RR is on
C
C      CNTL(3) is the threshold used in the following to compute  
C      DKEEP(9) the threshold under which the sing val. are considered
C      as null and from which we start to look for a gap between two
C      sing val.
         IF (cntl3 .LT. zero) THEN
           id%DKEEP(9) = abs(cntl(3))
         ELSE IF  (cntl3 .GT. zero) THEN
           id%DKEEP(9) = cntl3*anorminf
         ELSE  !  (CNTL(3) .EQ. ZERO) THEN
         ENDIF
         IF (prokg) THEN 
           WRITE(mpg, '(A,I16)')
     &    ' ICNTL(56) rank revealing effective value   =',keep(19)   
           WRITE(mpg,'(A,1PD16.4)')
     &    ' ...Threshold for singularities on the root =',id%DKEEP(9)
         ENDIF
C        RR postponing considers that pivot rows with norm smaller 
C        than SEUIL should be postponed. 
C        SEUIL should be bigger than DKEEP(9), this means that 
C        DKEEP(13) should be bigger than 1. 
         thresh_seuil = id%DKEEP(13)
         IF (id%DKEEP(13).LT.1)  thresh_seuil = 10
         seuil = id%DKEEP(9)*thresh_seuil
         IF (prokg) WRITE(mpg,'(A,1PD16.4)')
     &   ' ...Threshold for postponing                =',seuil
       ENDIF !end KEEP(19)
       seuil_ldlt_niv2 = seuil
C     -------------------------------
C     -- Null pivot row detection
C     -------------------------------
       IF (keep(110).EQ.0) THEN 
C        -- Null pivot is off
C        Initialize DKEEP(1) to a negative value
C        in order to avoid detection of null pivots
C        (test max(AMAX,RMAX,abs(PIVOT)).LE.PIVNUL
C        in ZMUMPS_FAC_I, where PIVNUL=DKEEP(1))
         id%DKEEP(1) = -1.0d0
         id%DKEEP(2) = zero
       ELSE
C        -- Null pivot is on
        IF (keep(19).NE.0) THEN
C        -- RR is on
C      RR postponing considers that pivot rows of norm smaller that SEUIL 
C      should be postponed, but pivot rows smaller than DKEEP(1) are 
C      directly added to null space and thus considered as null pivot rows. 
         IF ((id%DKEEP(10).LE.0).OR.(id%DKEEP(10).GT.1)) THEN
C           DKEEP(10) is out of range, set to the default value 10-1
           id%DKEEP(1) = id%DKEEP(9)*1d-1
         ELSE
           id%DKEEP(1) = id%DKEEP(9)*id%DKEEP(10) 
         ENDIF    
        ELSE
C        -- RR is off
C        -- only Null pivot detection
C         We keep strategy currently used in MUMPS 4.10.0
          IF (cntl3 .LT. zero) THEN
           id%DKEEP(1)  = abs(cntl(3))
          ELSE IF  (cntl3 .GT. zero) THEN
            id%DKEEP(1)  = cntl3*anorminf
          ELSE !  (CNTL(3) .EQ. ZERO) THEN
c          id%DKEEP(1)  = NPIV_CRITICAL_PATH*EPS*ANORMINF
           CALL mumps_npiv_critical_path(
     &         n, keep(28), id%STEP(1), id%FRERE_STEPS(1), id%FILS(1),
     &         id%NA(1), id%LNA, id%NE_STEPS(1), npiv_critical_path )
           id%DKEEP(1)  = sqrt(dble(npiv_critical_path))*eps*anorminf
          ENDIF
        ENDIF ! fin rank revealing        
        IF ((keep(110).NE.0).AND.(prokg)) THEN
           WRITE(mpg, '(A,I16)')
     &    ' ICNTL(24) null pivot rows detection        =',keep(110)   
           WRITE(mpg,'(A,1PD16.4)')
     &    ' ...Zero pivot detection threshold          =',id%DKEEP(1)
        ENDIF 
        IF (cntl5.GT.zero) THEN
            id%DKEEP(2) = cntl5 * anorminf
            IF (prokg) WRITE(mpg,'(A,1PD10.3)')
     &    ' ...Fixation for null pivots               =',id%DKEEP(2)
         ELSE
            IF (prokg) WRITE(mpg,*) '...Infinite fixation '
            IF (id%KEEP(50).EQ.0) THEN
C             Unsym
            ! the user let us choose a fixation. set in NEGATIVE
            ! to detect during facto when to set row to zero !
             id%DKEEP(2) = -max(1.0d10*anorminf, 
     &                sqrt(huge(anorminf))/1.0d8)
            ELSE
C             Sym
            id%DKEEP(2) = zero
            ENDIF
         ENDIF
       ENDIF ! fin null pivot detection.
C     Find id of root node if RR is on 
      IF (keep(53).NE.0) THEN
        id_root =mumps_procnode(id%PROCNODE_STEPS(id%STEP(keep(20))),
     &                          id%KEEP(199))
        IF ( keep( 46 )  .NE. 1 ) THEN
          id_root = id_root + 1
        END IF
      ENDIF
C Second pass:  set parameters for null pivot detection
C Allocate PIVNUL_LIST in case of null pivot detection
      lpn_list = 1
      IF ( associated( id%PIVNUL_LIST) ) DEALLOCATE(id%PIVNUL_LIST)
      IF(keep(110) .EQ. 1) THEN
         lpn_list = n
      ENDIF
      IF (keep(19).NE.0 .AND.
     &   (id_root.EQ.id%MYID .OR. id%MYID.EQ.master)) THEN
         lpn_list = n
      ENDIF
      ALLOCATE( id%PIVNUL_LIST(lpn_list),stat = ierr )
      IF ( ierr .GT. 0 ) THEN
        id%INFO(1)=-13
        id%INFO(2)=lpn_list
      END IF
      id%PIVNUL_LIST(1:lpn_list) = 0
      keep(109) = 0
C end set parameter for null pivot detection
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF ( id%INFO(1).lt.0 ) GOTO 517
C   --------------------------------------------------------------
C   STATIC PIVOTING
C     -- Static pivoting only when RR and Null pivot detection OFF
C   --------------------------------------------------------------
      keep(97) = 0
      IF ((keep(19).EQ.0).AND.(keep(110).EQ.0)) THEN
        IF (id%MYID .EQ. master) cntl4 = id%CNTL(4)
        CALL mpi_bcast( cntl4, 1, mpi_double_precision,
     &                master, id%COMM, ierr )
C 
        IF ( cntl4 .GE. zero ) THEN
         keep(97) = 1
         IF ( cntl4 .EQ. zero ) THEN
C           -- set seuil to sqrt(eps)*||A||
            IF(anorminf .EQ. zero) THEN
               eff_size_schur = 0
               IF (keep(60).GT.0) eff_size_schur = keep(116)
               CALL zmumps_anorminf( id , anorminf, lscal, 
     &              eff_size_schur )
            ENDIF
            seuil = sqrt(eps) * anorminf
         ELSE
            seuil = cntl4
         ENDIF
         seuil_ldlt_niv2 = seuil
        ELSE 
         seuil = zero
        ENDIF
      ENDIF
C     set number of tiny pivots / 2x2 pivots in types 1 /
C     2x2 pivots in types 2, to zero. This is because the
C     user can call the factorization step several times.
      keep(98)  = 0
      keep(103) = 0
      keep(105) = 0
      maxs      = 1_8
*
*     Start allocations
*     *****************
*
C
C  The slaves can now perform the factorization
C
C
C  Allocate id%S on all nodes
C  or point to user provided data WK_USER when LWK_USER>0
C  =======================
C
C     Compute BLR_STRAT and a first estimation 
C     of MAXS, the size of id%S
      CALL  zmumps_set_blrstrat_and_maxs_k8 (
     &           maxs_base8, maxs_base_relaxed8,
     &           blr_strat,
     &           id%KEEP(1), id%KEEP8(1))
C   
      maxs = maxs_base_relaxed8
      IF (wk_user_provided) THEN
C       -- Set MAXS to size of WK_USER_
        maxs = id%KEEP8(24)
      ENDIF
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF (id%INFO(1) .LT. 0) THEN
        GOTO 517
      ENDIF
C
      id%KEEP8(75) = huge(id%KEEP8(75))
      id%KEEP8(76) = huge(id%KEEP8(76))
      IF (i_am_slave) THEN
C
        IF (id%KEEP8(4) .NE. 0_8) THEN
C
          IF ( .NOT. wk_user_provided ) THEN
C           Set MAXS given BLR_STRAT, KEEP(201) and MAXS_BASE_RELAXED8
            CALL zmumps_mem_allowed_set_maxs (
     &          maxs,
     &          blr_strat, id%KEEP(201), maxs_base_relaxed8,
     &          id%KEEP(1), id%KEEP8(1), id%MYID, id%N, id%NELT, 
     &          id%NA(1), id%LNA, id%NSLAVES, 
     &          keep464copy,  keep465copy,
     &          id%INFO(1), id%INFO(2)
     &           , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &           size(id%I8_L0_OMP,2)
     &           )
C           Given  MAXS and  max memory allowed KEEP8(4)
C           compute in KEEP8(75) the number of real/complex 
C           available for dynamic allocations 
            CALL  zmumps_mem_allowed_set_k75 ( 
     &           maxs, id%MYID, 
     &           .false., ! UNDER_L0_OMP
     &           n, id%NELT, id%NA(1), id%LNA, id%NSLAVES,
     &           blr_strat, id%KEEP(201), 
     &           id%KEEP, id%KEEP8, id%INFO(1), id%INFO(2)
     &           , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &           size(id%I8_L0_OMP,2)
     &           )
          ELSE
C            KEEP8(75) dow not include MAXS, since WK_USER is provided
             CALL  zmumps_mem_allowed_set_k75 ( 
     &           0_8, id%MYID, 
     &           .false., ! UNDER_L0_OMP
     &           n, id%NELT, id%NA(1), id%LNA, id%NSLAVES,
     &           blr_strat, id%KEEP(201), 
     &           id%KEEP, id%KEEP8, id%INFO(1), id%INFO(2)
     &           , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &           size(id%I8_L0_OMP,2)
     &           )
          ENDIF
          IF (keep(400) .GT.0) THEN
C            ------------------------------
C            compute KEEP8(75) under L0_OMP
C            ------------------------------
C            Save KEEP8(75) above L0_OMP to reset KEEP8(75) 
C            when starting FAC_PAR_M
             id%KEEP8(76) = id%KEEP8(75)
             CALL  zmumps_mem_allowed_set_k75 ( 
     &             0_8, ! MAXS=0_8 
     &             id%MYID, 
     &             .true., ! UNDER_L0_OMP
     &             id%N, id%NELT, id%NA(1), id%LNA, id%NSLAVES,
     &             blr_strat, id%KEEP(201), 
     &             id%KEEP, id%KEEP8, id%INFO(1), id%INFO(2)
     &             , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &             size(id%I8_L0_OMP,2)
     &             )
C            KEEP8(75) holds the number of entries that 
C            can be allocated underL0.
C            It will be used during ZMUMPS_FAC_L0_OMP to adjust the 
C            the size of MUMPS_TPS_ARR(ITH)%LA 
          ENDIF
        ENDIF ! MEM_ALLOWED
C
      ENDIF ! I_AM_SLAVE THEN
C
      IF (i_am_slave) THEN
        IF ( (keep(400).GT.0) .AND. (keep(406).EQ.2) ) THEN
C         Compute KEEP8(77) the peak authorized used by 
C         ZMUMPS_PERFORM_COPIES
          CALL zmumps_l0_compute_peak_allowed(
     &          id%MYID, id%N,
     &          id%NELT, id%NA(1), id%LNA, id%NSLAVES,
     &          blr_strat, id%KEEP(201),
     &          id%KEEP(1), id%KEEP8(1), id%INFO(1), id%INFO(2)
     &          , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &          size(id%I8_L0_OMP,2)
     &           )
        ENDIF
      ENDIF ! I_AM_SLAVE)
C
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF (id%INFO(1) .LT. 0) THEN
        GOTO 517
      ENDIF
      CALL mumps_seti8toi4(maxs, id%INFO(39))
      CALL zmumps_avgmax_stat8(prokg, mpg, maxs, id%NSLAVES,
     & print_maxavg,
     & id%COMM, " Effective size of S     (based on INFO(39))=   ")
C
      IF ( i_am_slave ) THEN
C       ------------------
C       Dynamic scheduling
C       ------------------
        CALL zmumps_load_set_inicost( dble(id%COST_SUBTREES),
     &        keep(64), id%DKEEP(15), keep(375), maxs )
        k28=keep(28)
        memory_md_arg = min(int(perlu,8) * ( maxs_base8 / 100_8 + 1_8 ),
C       Restrict freedom from dynamic scheduler when
C       MEM_ALLOWED=ICNTL(23) is small (case where KEEP8(4)-MAXS_BASE8
C       is negative after call to ZMUMPS_MAX_MEM)
     &                      max(0_8, maxs-maxs_base8))
        CALL zmumps_load_init( id, memory_md_arg, maxs )
C
C       Out-Of-Core (OOC) issues. Case where we ran one factorization OOC
C       and the second one is in-core: we try to free OOC
C       related data from previous factorization.
C
        CALL zmumps_clean_ooc_data(id, ierr)
        IF (ierr < 0) THEN
          id%INFO(1) = -90
          id%INFO(2) = 0
          GOTO 112
        ENDIF
        IF (keep(201) .GT. 0) THEN
C          -------------------
C          OOC initializations
C          -------------------
           IF (keep(201).EQ.1 !PANEL Version
     &         .AND.keep(50).EQ.0 ! Unsymmetric
     &         .AND.keep(251).NE.2 ! Store L to disk
     &         ) THEN
              id%OOC_NB_FILE_TYPE=2 ! declared in MUMPS_OOC_COMMON
           ELSE
              id%OOC_NB_FILE_TYPE=1 ! declared in MUMPS_OOC_COMMON
           ENDIF
C          ------------------------------
C          Dimension IO buffer, KEEP(100)
C          ------------------------------
           IF (keep(205) .GT. 0) THEN
             keep(100) = keep(205)
           ELSE
             IF (keep(201).EQ.1) THEN ! PANEL version
               i8tmp = int(id%OOC_NB_FILE_TYPE,8) *
     &               2_8 * int(keep(226),8)
             ELSE
               i8tmp = 2_8 * id%KEEP8(119)
             ENDIF
             i8tmp = i8tmp +  int(max(keep(12),0),8) *
     &               (i8tmp/100_8+1_8)
C            we want to avoid too large IO buffers.
C            12M corresponds to 100Mbytes given to buffers.
             i8tmp = min(i8tmp, 12000000_8)
             keep(100)=int(i8tmp)
           ENDIF
           IF (keep(201).EQ.1) THEN
C            Panel version. Force the use of a buffer. 
             IF ( keep(99) < 3 ) THEN
               keep(99) = keep(99) + 3
             ENDIF
           ENDIF
C          --------------------------
C          Reset KEEP(100) to 0 if no
C          buffer is used for OOC.
C          --------------------------
           IF (keep(99) .LT.3) keep(100)=0
           IF((dble(keep(100))*dble(keep(35))/dble(2)).GT.
     &          (dble(1999999999)))THEN
             IF (prokg) THEN
               WRITE(mpg,*)id%MYID,': Warning: DIM_BUF_IO might be
     &  too big for Filesystem'
             ENDIF
           ENDIF
           ALLOCATE (id%OOC_INODE_SEQUENCE(keep(28),
     &          id%OOC_NB_FILE_TYPE),
     &          stat=ierr)
           IF ( ierr .GT. 0 ) THEN
              id%INFO(1) = -13
              id%INFO(2) = id%OOC_NB_FILE_TYPE*keep(28)
              NULLIFY(id%OOC_INODE_SEQUENCE)
              GOTO 112
           ENDIF
           ALLOCATE (id%OOC_TOTAL_NB_NODES(id%OOC_NB_FILE_TYPE),
     &          stat=ierr)
           IF ( ierr .GT. 0 ) THEN
              id%INFO(1) = -13
              id%INFO(2) = id%OOC_NB_FILE_TYPE
              NULLIFY(id%OOC_TOTAL_NB_NODES)
              GOTO 112
           ENDIF
           ALLOCATE (id%OOC_SIZE_OF_BLOCK(keep(28),
     &          id%OOC_NB_FILE_TYPE),
     &          stat=ierr)
           IF ( ierr .GT. 0 ) THEN
              id%INFO(1) = -13
              id%INFO(2) = id%OOC_NB_FILE_TYPE*keep(28)
              NULLIFY(id%OOC_SIZE_OF_BLOCK)
              GOTO 112
           ENDIF
           ALLOCATE (id%OOC_VADDR(keep(28),id%OOC_NB_FILE_TYPE),
     &          stat=ierr)
           IF ( ierr .GT. 0 ) THEN
              id%INFO(1) = -13
              id%INFO(2) = id%OOC_NB_FILE_TYPE*keep(28)
              NULLIFY(id%OOC_VADDR)
              GOTO 112
           ENDIF
        ENDIF
      ENDIF
 112  CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF (id%INFO(1) < 0) THEN
C       LOAD_END must be done but not OOC_END_FACTO
        GOTO 513
      ENDIF
      IF (i_am_slave) THEN
        IF (keep(201) .GT. 0) THEN
           IF ((keep(201).EQ.1).OR.(keep(201).EQ.2)) THEN
             CALL zmumps_ooc_init_facto(id,maxs)
          ELSE
             WRITE(*,*) "Internal error in ZMUMPS_FAC_DRIVER"
             CALL mumps_abort()
           ENDIF
           IF(id%INFO(1).LT.0)THEN
              GOTO 111
           ENDIF
        ENDIF
C       First increment corresponds to the number of
C       floating-point operations for subtrees allocated
C       to the local processor.
        CALL zmumps_load_update(0,.false.,dble(id%COST_SUBTREES),
     &          id%KEEP(1),id%KEEP8(1))
        IF (id%INFO(1).LT.0) GOTO 111
      END IF
C     -----------------------
C     Manage main workarray S
C     -----------------------
      earlyt3rootins = keep(200) .EQ.0
     &        .OR.  ( keep(200) .LT. 0 .AND. keep(400) .EQ. 0 )
#if defined (LARGEMATRICES)
      IF ( id%MYID .ne. master ) THEN
#endif
      IF (.NOT.wk_user_provided) THEN
        IF ( earlyt3rootins ) THEN
C         Standard allocation strategy
          ALLOCATE (id%S(maxs),stat=ierr)
          id%KEEP8(23) = maxs
          IF ( ierr .GT. 0 ) THEN
            id%INFO(1) = -13
            CALL mumps_set_ierror(maxs, id%INFO(2))
C           On some platforms (IBM for example), an
C           allocation failure returns a non-null pointer.
C           Therefore we nullify S
            NULLIFY(id%S)
            id%KEEP8(23)=0_8
          ENDIF
        ENDIF
      ELSE
       id%S => id%WK_USER(1:id%KEEP8(24))
       id%KEEP8(23) = 0_8
      ENDIF
#if defined (LARGEMATRICES)
      END IF
#endif
C
 111  CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF ( id%INFO(1).LT.0 ) GOTO 514
C     --------------------------
C     Initialization of modules
C     related to data management
C     --------------------------
      nb_active_fronts_estim = 3
      nb_threads = 1
!$    NB_THREADS = OMP_GET_MAX_THREADS()
C     
      nb_active_fronts_estim = 3*nb_threads
      IF (i_am_slave) THEN
C
        CALL mumps_fdm_init('A',nb_active_fronts_estim, id%INFO)
C
        IF ( (keep(486).EQ.2) 
     &     .OR. ((keep(489).NE.0).AND.(keep(400).GT.1))
     &    ) THEN
C         In case of LRSOLVE or CompressCB, 
C         initialize nb of handlers to nb of BLR
C         nodes estimated at analysis
          nb_fronts_f_estim = keep(470)
        ELSE
          IF (keep(489).NE.0) THEN
C          Compress CB and no L0 OMP (or 1 thread under L0): 
C          NB_ACTIVE_FRONTS_ESTIM is too small, 
C          to limit nb of reallocations make it twice larger
           nb_fronts_f_estim = 2*nb_active_fronts_estim
          ELSE
           nb_fronts_f_estim = nb_active_fronts_estim
          ENDIF
        ENDIF
        CALL mumps_fdm_init('F',nb_fronts_f_estim, id%INFO )
        IF (id%INFO(1) .LT. 0 ) GOTO 114
#if ! defined(NO_FDM_DESCBAND)
C         Storage of DESCBAND information
          CALL mumps_fdbd_init( nb_active_fronts_estim, id%INFO )
#endif
#if ! defined(NO_FDM_MAPROW)
C         Storage of MAPROW and ROOT2SON information
          CALL mumps_fmrd_init( nb_active_fronts_estim, id%INFO )
#endif
        CALL zmumps_blr_init_module( nb_fronts_f_estim, id%INFO )
 114    CONTINUE
      ENDIF
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                       id%COMM, id%MYID )
C     GOTO 500: one of the above module initializations failed
      IF ( id%INFO(1).LT.0 ) GOTO 500
C
C
C  Allocate space for matrix in arrowhead 
C  ======================================
C
C  CASE 1 : Matrix is assembled
C  CASE 2 : Matrix is elemental
C
      IF ( keep(55) .eq. 0 ) THEN
C       ------------------------------------
C       Space has been allocated already for
C       the integer part during analysis
C       Only slaves need the arrowheads.
C       ------------------------------------
        IF (associated( id%DBLARR)) THEN
          DEALLOCATE(id%DBLARR)
          NULLIFY(id%DBLARR)
        ENDIF
        IF ( i_am_slave .and. id%KEEP8(26) .ne. 0_8 ) THEN
          ALLOCATE( id%DBLARR( id%KEEP8(26) ), stat = ierr )
        ELSE
          ALLOCATE( id%DBLARR( 1 ), stat =ierr )
        END IF
        IF ( ierr .NE. 0 ) THEN
          IF (lpok) THEN
            WRITE(lp,*) id%MYID,
     &      ': Allocation error for DBLARR(',id%KEEP8(26),')'
          ENDIF
          id%INFO(1)=-13
          CALL mumps_set_ierror(id%KEEP8(26), id%INFO(2))
          NULLIFY(id%DBLARR)
          GOTO 100
        END IF
      ELSE
C        ----------------------------------------
C        Allocate variable lists. Systematically.
C        ----------------------------------------
         IF ( associated( id%INTARR ) ) THEN
           DEALLOCATE( id%INTARR )
           NULLIFY( id%INTARR )
         END IF
         IF ( i_am_slave .and. id%KEEP8(27) .ne. 0_8 ) THEN
           ALLOCATE( id%INTARR( id%KEEP8(27) ), stat = allocok )
           IF ( allocok .GT. 0 ) THEN
             id%INFO(1) = -13
             CALL mumps_set_ierror(id%KEEP8(27), id%INFO(2))
             NULLIFY(id%INTARR)
             GOTO 100
           END IF
         ELSE
           ALLOCATE( id%INTARR(1),stat=allocok )
           IF ( allocok .GT. 0 ) THEN
             id%INFO(1) = -13
             id%INFO(2) = 1
             NULLIFY(id%INTARR)
             GOTO 100
           END IF
         END IF
C        -----------------------------
C        Allocate real values.
C        On master, if hybrid host and
C        no scaling, avoid the copy.
C        -----------------------------
         IF (associated( id%DBLARR)) THEN
           DEALLOCATE(id%DBLARR)
           NULLIFY(id%DBLARR)
         ENDIF
         IF ( i_am_slave ) THEN
           IF (      id%MYID_NODES .eq. master
     &       .AND.   keep(46)   .eq. 1
     &       .AND.   keep(52)   .eq. 0 ) THEN
C            --------------------------
C            Simple pointer association
C            --------------------------
             id%DBLARR => id%A_ELT
           ELSE
C            ----------
C            Allocation
C            ----------
             IF ( id%KEEP8(26) .ne. 0_8 ) THEN
               ALLOCATE( id%DBLARR( id%KEEP8(26) ), stat = allocok )
               IF ( allocok .GT. 0 ) THEN
                 id%INFO(1) = -13
                 CALL mumps_set_ierror(id%KEEP8(26), id%INFO(2))
                 NULLIFY(id%DBLARR)
                 GOTO 100
               END IF
             ELSE
               ALLOCATE( id%DBLARR(1), stat = allocok )
               IF ( allocok .GT. 0 ) THEN
                 id%INFO(1) = -13
                 id%INFO(2) = 1
                 NULLIFY(id%DBLARR)
                 GOTO 100
               END IF 
             END IF
           END IF
         ELSE
           ALLOCATE( id%DBLARR(1), stat = allocok )
           IF ( allocok .GT. 0 ) THEN
             id%INFO(1) = -13
             id%INFO(2) = 1
             NULLIFY(id%DBLARR)
             GOTO 100
           END IF
         END IF
      END IF
C     -----------------
C     Also prepare some
C     data for the root
C     -----------------
      IF ( keep(38).NE.0 .AND.  i_am_slave ) THEN
         CALL zmumps_init_root_fac( id%N,
     &   id%root, id%FILS(1), keep(38), id%KEEP(1), id%INFO(1) )
      END IF
C
C
 100  CONTINUE
C     ----------------
C     Check for errors
C     ----------------
      CALL mumps_propinfo( id%ICNTL(1), id%INFO(1),
     &                        id%COMM, id%MYID )
      IF ( id%INFO(1).LT.0 ) GOTO 500
C
C       -----------------------------------
C
C       DISTRIBUTION OF THE ORIGINAL MATRIX
C
C       -----------------------------------
C
C     TIMINGS: computed (and printed) on the host
C     Next line: global time for distrib(arrowheads,elts)
C     on the host. Synchronization has been performed.
      IF (id%MYID.EQ.master) CALL mumps_secdeb(time)
C     -------------------------------------------
C     S_PTR_ARG / MAXS_ARG will be used for id%S
C     argument to arrowhead/element distribution
C     routines: if id%S is not allocated, we pass
C     S_DUMMY_ARG instead, which is not accessed.
C     -------------------------------------------
      IF (earlyt3rootins) THEN
            s_ptr_arg => id%S
            maxs_arg = maxs
      ELSE
            s_ptr_arg => s_dummy_arg
            maxs_arg = 1
      ENDIF
C
      IF ( keep( 55 ) .eq. 0 ) THEN
C       ----------------------------
C       Original matrix is assembled
C       Arrowhead format to be used.
C       ----------------------------
C       KEEP8(26) and KEEP8(27) hold the number of entries for real/integer
C       for the matrix in arrowhead format. They have been set by the
C       analysis phase (ZMUMPS_ANA_F and ZMUMPS_ANA_G)
C
C       ------------------------------------------------------------------
C       Blocking is used for sending arrowhead records (I,J,VAL)
C              buffer(1) is used to store number of bytes already packed
C              buffer(2) number of records already packed
C       KEEP(39) : Number of records (blocking factor)
C       ------------------------------------------------------------------
C
C     ---------------------------------------------
C     In case of parallel root compute minimum
C     size of workspace to receive arrowheads
C     of root node. Will be used to check that
C     MAXS is large enough for arrowheads (case
C     of EARLYT3ROOTINS (KEEP(200)=0); if .NOT.
C     EARLYT3ROOTINS (KEEP(200)=1), root will
C     be assembled into id%S later and size of
C     id%S will be checked later)
C     ---------------------------------------------
      IF (earlyt3rootins .AND. keep(38).NE.0 .AND.
     &    keep(60) .EQ.0 .AND. i_am_slave) THEN
        lwk = int(numroc( id%root%ROOT_SIZE, id%root%MBLOCK,
     &             id%root%MYROW, 0, id%root%NPROW ),8)
        lwk = max( 1_8, lwk )
        lwk = lwk*
     &        int(numroc( id%root%ROOT_SIZE, id%root%NBLOCK,
     &        id%root%MYCOL, 0, id%root%NPCOL ),8)
        lwk = max( 1_8, lwk )
      ELSE
        lwk = 1_8
      ENDIF
C     MAXS must be at least 1, and in case of
C     parallel root, large enough to receive
C     arrowheads of root.
      IF (maxs .LT. int(lwk,8)) THEN
           id%INFO(1) = -9
           CALL mumps_set_ierror(lwk, id%INFO(2))
      ENDIF
      CALL mumps_propinfo( id%ICNTL(1), id%INFO(1),
     &                        id%COMM, id%MYID )
      IF ( id%INFO(1).LT.0 ) GOTO 500
C
      IF ( keep(54) .eq. 0 ) THEN
C       ================================================
C       FIRST CASE : MATRIX IS NOT INITIALLY DISTRIBUTED
C       ================================================
C       A small integer workspace is needed to
C       send the arrowheads.
        IF ( id%MYID .eq. master ) THEN
          ALLOCATE(iwk(id%N), stat=allocok)
          IF ( allocok .NE. 0 ) THEN
            id%INFO(1)=-13
            id%INFO(2)=id%N
          END IF
#if defined(LARGEMATRICES)
          ALLOCATE (wk(lwk),stat=ierr)
          IF ( ierr .GT. 0 ) THEN
            id%INFO(1) = -13
            CALL mumps_set_ierror(lwk, id%INFO(2))
            write(6,*) ' PB1 ALLOC LARGEMAT'
          ENDIF
#endif
        ENDIF
        CALL mumps_propinfo( id%ICNTL(1), id%INFO(1),
     &                        id%COMM, id%MYID )
        IF ( id%INFO(1).LT.0 ) GOTO 500
        IF ( id%MYID .eq. master ) THEN
C
C         --------------------------------
C         MASTER sends arowheads using the
C         global communicator with ranks
C         also in global communicator
C         IWK is used as temporary
C         workspace of size N.
C         --------------------------------
           IF ( .not. associated( id%INTARR ) ) THEN
              ALLOCATE( id%INTARR( 1 ),stat=ierr)
              IF ( ierr .GT. 0 ) THEN
                 id%INFO(1) = -13
                 id%INFO(2) = 1
                 NULLIFY(id%INTARR)
                 write(6,*) ' PB2 ALLOC INTARR'
                 CALL mumps_abort()
              ENDIF
           ENDIF
          nbrecords = keep(39)
          IF (id%KEEP8(28) .LT. int(nbrecords,8)) THEN
            nbrecords = int(id%KEEP8(28))
          ENDIF
#if defined(LARGEMATRICES)
          CALL zmumps_facto_send_arrowheads(id%N, id%KEEP8(28), id%A(1),
     &      id%IRN(1), id%JCN(1), id%SYM_PERM(1),
     &      lscal, id%COLSCA(1), id%ROWSCA(1),   
     &      id%MYID, id%NSLAVES, id%PROCNODE_STEPS(1),
     &      nbrecords,
     &      lp, id%COMM, id%root, keep,id%KEEP8,
     &      id%FILS(1), iwk(1), ! workspace of size N
     &
     &      id%INTARR(1), id%KEEP8(27), id%DBLARR(1), id%KEEP8(26),
     &      id%PTRAR(1), id%PTRAR(id%N+1),
     &      id%FRERE_STEPS(1), id%STEP(1), wk(1), lwk,
     &      id%ISTEP_TO_INIV2(1), id%I_AM_CAND(1),
     &      id%CANDIDATES(1,1)) 
C         write(6,*) '!!! A,IRN,JCN are freed during factorization '
          DEALLOCATE (id%A)
          NULLIFY(id%A)
          DEALLOCATE (id%IRN)
          NULLIFY (id%IRN)
          DEALLOCATE (id%JCN)
          NULLIFY (id%JCN)
          IF (.NOT.wk_user_provided) THEN
            IF (earlyt3rootins) THEN
              ALLOCATE (id%S(maxs),stat=ierr)
              id%KEEP8(23) = maxs
              IF ( ierr .GT. 0 ) THEN
                id%INFO(1) = -13
                id%INFO(2) = maxs
                NULLIFY(id%S)
                id%KEEP8(23)=0_8
                write(6,*) ' PB2 ALLOC LARGEMAT',maxs
                CALL mumps_abort()
               ENDIF
             ENDIF
            ENDIF
          ELSE
            id%S => id%WK_USER(1:id%KEEP8(24))
          ENDIF
          IF (earlyt3rootins) THEN
            id%S(maxs-lwk+1_8:maxs) = wk(1_8:lwk)
          ENDIF
          DEALLOCATE (wk)
#else
          CALL zmumps_facto_send_arrowheads(id%N, id%KEEP8(28), id%A(1),
     &    id%IRN(1), id%JCN(1), id%SYM_PERM(1),
     &    lscal, id%COLSCA(1), id%ROWSCA(1),   
     &    id%MYID, id%NSLAVES, id%PROCNODE_STEPS(1),
     &    nbrecords,
     &    lp, id%COMM, id%root, keep(1),id%KEEP8(1),
     &    id%FILS(1), iwk(1),
     &
     &    id%INTARR(1), id%KEEP8(27), id%DBLARR(1), id%KEEP8(26),
     &    id%PTRAR(1), id%PTRAR(id%N+1),
     &    id%FRERE_STEPS(1), id%STEP(1), s_ptr_arg(1), maxs_arg,
     &    id%ISTEP_TO_INIV2(1), id%I_AM_CAND(1),
     &    id%CANDIDATES(1,1) ) 
#endif
          DEALLOCATE(iwk)
        ELSE
          nbrecords = keep(39)
          IF (id%KEEP8(28) .LT. int(nbrecords,8)) THEN
            nbrecords = int(id%KEEP8(28))
          ENDIF
          CALL zmumps_facto_recv_arrowhd2( id%N,
     &       id%DBLARR(1), id%KEEP8(26),
     &       id%INTARR(1), id%KEEP8(27),
     &       id%PTRAR( 1 ),
     &       id%PTRAR(id%N+1),
     &       keep( 1 ), id%KEEP8(1), id%MYID, id%COMM,
     &       nbrecords,
     &
     &       s_ptr_arg(1), maxs_arg,
     &       id%root,
     &       id%PROCNODE_STEPS(1), id%NSLAVES,
     &       id%SYM_PERM(1), id%FRERE_STEPS(1), id%STEP(1),
     &       id%INFO(1), id%INFO(2) )
        ENDIF
      ELSE
C
C     =============================================
C     SECOND CASE : MATRIX IS INITIALLY DISTRIBUTED
C     =============================================
C     Timing on master.
      IF (id%MYID.EQ.master) THEN
        CALL mumps_secdeb(time)
      END IF
      IF ( i_am_slave ) THEN
C       ---------------------------------------------------
C       In order to have possibly IRN_loc/JCN_loc/A_loc
C       of size 0, avoid to pass them inside REDISTRIBUTION
C       and pass id instead
C       NZ_locMAX8 gives as a maximum buffer size (send/recv) used
C        an upper bound to limit buffers on small matrices
C       ---------------------------------------------------
       CALL mpi_allreduce(id%KEEP8(29), nz_locmax8, 1, mpi_integer8,
     &                   mpi_max, id%COMM_NODES, ierr)
       nbrecords = keep(39)
       IF (nz_locmax8 .LT. int(nbrecords,8)) THEN
            nbrecords = int(nz_locmax8)
       ENDIF
        CALL zmumps_redistribution( id%N,
     &  id%KEEP8(29),
     &  id,
     &  id%DBLARR(1), id%KEEP8(26), id%INTARR(1),
     &  id%KEEP8(27), id%PTRAR(1), id%PTRAR(id%N+1),
     &  keep(1), id%KEEP8(1), id%MYID_NODES,
     &  id%COMM_NODES, nbrecords,
     &  s_ptr_arg(1), maxs_arg, id%root, id%PROCNODE_STEPS(1),
     &  id%NSLAVES, id%SYM_PERM(1), id%STEP(1),
     &  id%ICNTL(1), id%INFO(1), nsend8, nlocal8,
     &  id%ISTEP_TO_INIV2(1),
     &  id%CANDIDATES(1,1) )
        IF ( ( keep(52).EQ.7 ).OR. (keep(52).EQ.8) ) THEN
C         -------------------------------------------------
C         In that case, scaling arrays have been allocated
C         on all processors. They were useful for matrix
C         distribution. But we now really only need them
C         on the host. In case of distributed solution, we
C         will have to broadcast either ROWSCA or COLSCA
C         (depending on MTYPE) but this is done later.
C
C         In other words, on exit from the factorization,
C         we want to have scaling arrays available only
C         on the host.
C         -------------------------------------------------
          IF ( id%MYID > 0 ) THEN
            IF (associated(id%ROWSCA)) THEN
              DEALLOCATE(id%ROWSCA)
              NULLIFY(id%ROWSCA)
            ENDIF
            IF (associated(id%COLSCA)) THEN
              DEALLOCATE(id%COLSCA)
              NULLIFY(id%COLSCA)
            ENDIF
          ENDIF
        ENDIF
#if defined(LARGEMATRICES)
C      deallocate id%IRN_loc, id%JCN(loc) to free extra space
C      Note that in this case IRN_loc cannot be used
C      anymore during the solve phase for IR and Error analysis.
         IF (associated(id%IRN_loc)) THEN
            DEALLOCATE(id%IRN_loc)
            NULLIFY(id%IRN_loc)
         ENDIF
         IF (associated(id%JCN_loc)) THEN
            DEALLOCATE(id%JCN_loc)
            NULLIFY(id%JCN_loc)
         ENDIF
         IF (associated(id%A_loc)) THEN
            DEALLOCATE(id%A_loc)
            NULLIFY(id%A_loc)
         ENDIF
       write(6,*) ' Warning :', 
     &        ' id%A_loc, IRN_loc, JCN_loc deallocated !!! '
#endif
      IF (prok) THEN
        WRITE(mp,120) nlocal8, nsend8
      END IF
      END IF
      IF ( keep(46) .eq. 0 .AND. id%MYID.eq.master ) THEN
C       ------------------------------
C       The host is not working -> had
C       no data from initial matrix
C       ------------------------------
        nsend8  = 0_8
        nlocal8 = 0_8
      END IF
C     --------------------------
C     Put into some info/infog ?
C     --------------------------
      CALL mpi_reduce( nsend8, nsend_tot8, 1, mpi_integer8,
     &                 mpi_sum, master, id%COMM, ierr )
      CALL mpi_reduce( nlocal8, nlocal_tot8, 1, mpi_integer8,
     &                 mpi_sum, master, id%COMM, ierr )
      IF ( prokg ) THEN
        WRITE(mpg,125) nlocal_tot8, nsend_tot8
      END IF
C
C     -------------------------
C     Check for possible errors
C     -------------------------
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF ( id%INFO( 1 ) .LT. 0 ) GOTO 500
C
      ENDIF
      ELSE
C       -------------------
C       Matrix is elemental,
C       provided on the
C       master only
C       -------------------
        IF ( id%MYID.eq.master)
     &   CALL zmumps_maxelt_size( id%ELTPTR(1),
     &                        id%NELT,
     &                        maxelt_size )
C
C         Perform the distribution of the elements.
C         A this point,
C           PTRAIW/PTRARW have been computed.
C           INTARR/DBLARR have been allocated
C           ELTPROC gives the mapping of elements
C
        CALL zmumps_elt_distrib( id%N, id%NELT, id%KEEP8(30),
     &     id%COMM, id%MYID,
     &     id%NSLAVES, id%PTRAR(1),
     &     id%PTRAR(id%NELT+2),
     &     id%INTARR(1), id%DBLARR(1), id%KEEP8(27), id%KEEP8(26),
     &     id%KEEP(1), id%KEEP8(1), maxelt_size,
     &     id%FRTPTR(1), id%FRTELT(1),
     &     s_ptr_arg(1), maxs_arg, id%FILS(1),
     &     id, id%root )
C       ----------------
C       Broadcast errors
C       ----------------
        CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
        IF ( id%INFO( 1 ) .LT. 0 ) GOTO 500
      END IF ! Element entry
C     ------------------------
C     Time the redistribution:
C     ------------------------
      IF ( id%MYID.EQ.master) THEN
        CALL mumps_secfin(time)
        id%DKEEP(93) = time
        IF (prokg) WRITE(mpg,160) id%DKEEP(93)
      END IF
      IF ( keep(400) .GT. 0 ) THEN
C       L0-OMP was active at analysis and
C       thus will be active at factorization
C       We check the number of threads.
        nomp=1
!$      NOMP = omp_get_max_threads()
          IF ( nomp .NE. keep(400) ) THEN
            id%INFO(1)=-58
            id%INFO(2)=keep(400)
            IF (lpok) WRITE(lp,'(A,A,I5,A,I5)')
     &" FAILURE DETECTED IN FACTORIZATION: #threads for KEEP(401)",
     &" changed from",keep(400)," at analysis to", nomp
          ENDIF
C       error check
        CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
        IF ( id%INFO( 1 ) .LT. 0 ) GOTO 500
      ENDIF
C
C     TIMINGS:
C     Next line: elapsed time for factorization
      IF (id%MYID.EQ.master)  CALL mumps_secdeb(time)
C
C  Allocate buffers on the workers
C  ===============================
C
      IF ( i_am_slave )  THEN
        CALL zmumps_buf_ini_myid(id%MYID_NODES)
C
C  Some buffers are required to pack/unpack data and for
C  receiving MPI messages.
C  For packing/unpacking : the buffer must be large
C  enough to send several messages while receives might not
C  be posted yet.
C  It is assumed that the size of an integer is held in KEEP(34)
C  while the size of a complex is held in KEEP(35).
C  BUFR and LBUFR are declared of type integer, since byte is not
C  a standard datatype.
C  We now use KEEP(43) or KEEP(379) and KEEP(44) or KEEP(380)
C  as estimated at analysis to allocate appropriate buffer sizes
C
C  Reception buffer
C  ----------------
        IF (keep(486).NE.0) THEN
          zmumps_lbufr_bytes8 = int(keep( 380 ),8) * int(keep( 35 ),8)
        ELSE
          zmumps_lbufr_bytes8 = int(keep( 44 ),8) * int(keep( 35 ),8)
        ENDIF
C       ---------------------------------------
C       Ensure a reasonable minimal buffer size
C       ---------------------------------------
        zmumps_lbufr_bytes8 = max( zmumps_lbufr_bytes8,
     &                      100000_8 )
        IF ((keep(50).NE.0).AND.(keep(489).GT.0).AND.
     &       (id%NSLAVES.GE.2)) THEN
C        ----------------------------------------------------------
C        Ensure large enough receive buffer in case of BLR with
C        CompressCB for symmetric matrices. 
C        -----------------------------------------------------------
          ratiok465 = dble(keep465copy)/dble(1000)
          zmumps_lbufr_bytes8 = max(zmumps_lbufr_bytes8, 
     &     int( 
     &       ratiok465*
     &       dble(
     &        int(keep(2)+1,8)*int(keep(142),8)*int(keep(35),8) 
     &           ) 
     &       ,8)
     &                              )
        ENDIF
C
C  If there is pivoting, size of the message might still increase.
C  We use a relaxation (so called PERLU) to increase the estimate.
C
C  Note: PERLU is a global estimate for pivoting. 
C  It may happen that one large contribution block size is increased 
C  by more than that.
C  This is why we use an extra factor 2 relaxation coefficient for 
C  the relaxation of
C  the reception buffer in the case where pivoting is allowed.
C  A more dynamic strategy could be applied: if message to
C  be received is larger than expected, reallocate a larger
C  buffer. (But this won't work with IRECV.)
C  Finally, one may want (as we are currently doing it for 
C  most messages)
C  to cut large messages into a series of smaller ones.
C
        IF (keep(48).EQ.5) THEN
          min_perlu = 2
        ELSE
          min_perlu = 0
        ENDIF
C
        zmumps_lbufr_bytes8 = zmumps_lbufr_bytes8
     &        + int( 2.0d0 * dble(max(perlu,min_perlu))*
     &        dble(zmumps_lbufr_bytes8)/100d0, 8)
        zmumps_lbufr_bytes8 = min(zmumps_lbufr_bytes8,
     &                            int(huge(i4)-100,8))
        zmumps_lbufr_bytes  = int( zmumps_lbufr_bytes8 )
        IF (keep(48)==5) THEN
C          Since the buffer is going to be allocated, use
C          it as the constraint for memory/granularity
C          in hybrid scheduler
C
           id%KEEP8(21) = id%KEEP8(22) +
     &        int( dble(max(perlu,min_perlu))*
     &        dble(id%KEEP8(22))/100d0,8)
        ENDIF
C
C  Now estimate the size for the buffer for asynchronous
C  sends of contribution blocks (so called CB). We want to be able to send at
C  least KEEP(213)/100 (two in general) messages at the
C  same time.
C
C   Send buffer
C   -----------
        IF (keep(486).NE.0) THEN
         zmumps_lbuf8 = int( dble(keep(213)) / 100.0d0 *
     &                      dble(keep(379)) * dble(keep(35)), 8  )
        ELSE
         zmumps_lbuf8 = int( dble(keep(213)) / 100.0d0 *
     &                      dble(keep(43)) * dble(keep(35)), 8  )
        ENDIF
        zmumps_lbuf8 = max( zmumps_lbuf8, 100000_8 )
        zmumps_lbuf8 = zmumps_lbuf8
     &                 + int( 2.0d0 * dble(max(perlu,min_perlu))*
     &                   dble(zmumps_lbuf8)/100d0, 8)
C       Make ZMUMPS_LBUF8 small enough to be stored in a standard integer
        zmumps_lbuf8 = min(zmumps_lbuf8, int(huge(i4)-100,8))
C
C       No reason to have send buffer smaller than receive buffer.
C       This should never occur with the formulas above but just
C       in case:
        zmumps_lbuf8 = max(zmumps_lbuf8, zmumps_lbufr_bytes8+3*keep(34))
        zmumps_lbuf  = int(zmumps_lbuf8)
        IF(id%KEEP(48).EQ.4)THEN
           zmumps_lbufr_bytes=zmumps_lbufr_bytes*5
           zmumps_lbuf=zmumps_lbuf*5
        ENDIF
C
C  Estimate size of buffer for small messages 
C  Each node can send ( NSLAVES - 1 ) messages to (NSLAVES-1) nodes
C
C  KEEP(56) is the number of nodes of level II.
C  Messages will be sent for the symmetric case
C  for synchronisation issues.
C
C  We take an upperbound
C
        zmumps_lbuf_int = ( keep(56) + id%NSLAVES * id%NSLAVES ) * 5
     &               * keep(34)
        IF ( keep( 38 ) .NE. 0 ) THEN
C
C
          kkkk = mumps_procnode( id%PROCNODE_STEPS(id%STEP(keep(38))),
     &                           id%KEEP(199) )
          IF ( kkkk .EQ. id%MYID_NODES ) THEN
             zmumps_lbuf_int = zmumps_lbuf_int + 4 * keep(34) *
     &         ( id%NSLAVES + id%NE_STEPS(id%STEP(keep(38)))
     &      + min(keep(56), id%NE_STEPS(id%STEP(keep(38)))) * id%NSLAVES
     &         )
          END IF
        END IF
C       At this point, ZMUMPS_LBUFR_BYTES, ZMUMPS_LBUF
C       and ZMUMPS_LBUF_INT have been computed (all
C       are in numbers of bytes).
        IF ( prok ) THEN
          WRITE( mp, 9999 ) zmumps_lbufr_bytes,
     &                      zmumps_lbuf, zmumps_lbuf_int
        END IF
 9999   FORMAT( /,' Allocated buffers',/,' ------------------',/,
     &  ' Size of reception buffer in bytes ...... = ', i10,
     &  /,
     &  ' Size of async. emission buffer (bytes).. = ', i10,/,
     &  ' Small emission buffer (bytes) .......... = ', i10)
C       --------------------------
C       Allocate small send buffer
C       required for ZMUMPS_FAC_B
C       --------------------------
        CALL zmumps_buf_alloc_small_buf( zmumps_lbuf_int, ierr )
        IF ( ierr .NE. 0 ) THEN
          id%INFO(1)= -13
C         convert to size in integer  id%INFO(2)= ZMUMPS_LBUF_INT
          id%INFO(2)= (zmumps_lbuf_int+keep(34)-1)/keep(34)
          IF (lpok) THEN
            WRITE(lp,*) id%MYID,
     &     ':Allocation error in ZMUMPS_BUF_ALLOC_SMALL_BUF'
     &     ,id%INFO(2)
          ENDIF
          GO TO 110
        END IF
C
C       --------------------------------------
C       Allocate reception buffer on all procs
C       This is done now.
C       --------------------------------------
        zmumps_lbufr = (zmumps_lbufr_bytes+keep(34)-1)/keep(34)
        ALLOCATE( bufr( zmumps_lbufr ),stat=ierr )
        IF ( ierr .NE. 0 ) THEN
          id%INFO(1) = -13
          id%INFO(2) = zmumps_lbufr
          IF (lpok) THEN
            WRITE(lp,*)
     &     ': Allocation error for BUFR(', zmumps_lbufr,
     &     ') on MPI process',id%MYID
          ENDIF
          GO TO 110
        END IF
C       -----------------------------------------
C       Estimate MAXIS. IS will be allocated in
C       ZMUMPS_FAC_B. It will contain factors and
C       contribution blocks integer information
C       -----------------------------------------
C       Relax integer workspace based on PERLU
        perlu          = keep( 12 )
        IF (keep(201).GT.0) THEN
C         OOC panel or non panel (note that
C         KEEP(15)=KEEP(225) if non panel)
          maxis_estim   = keep(225)
        ELSE
C         In-core or reals for factors not stored
          maxis_estim   = keep(15)
        ENDIF
        maxis = max( 1, int( min(  int(huge(maxis),8),
     &          int(maxis_estim,8) + 3_8 * max(int(perlu,8),10_8) *
     &          ( int(maxis_estim,8) / 100_8 + 1_8 )
     &                          )   ! min
     &                     )  ! int
     &             ) !max
C       ----------------------------
C       Allocate PTLUST_S and PTRFAC
C       They will be used to access
C       factors in the solve phase.
C       They are also needed for
C       ZMUMPS_FAC_L0_OMP.
C       ----------------------------
        ALLOCATE( id%PTLUST_S( id%KEEP(28) ), stat = ierr )
        IF ( ierr .NE. 0 ) THEN
          id%INFO(1)=-13
          id%INFO(2)=id%KEEP(28)
          IF (lpok) THEN
            WRITE(lp,*) id%MYID,
     &      ': Allocation error for id%PTLUST_S(', id%KEEP(28),')'
          ENDIF
          NULLIFY(id%PTLUST_S)
          GOTO 110
        END IF
        ALLOCATE( id%PTRFAC( id%KEEP(28) ), stat = ierr )
        IF ( ierr .NE. 0 ) THEN
          id%INFO(1)=-13
          id%INFO(2)=id%KEEP(28)
          NULLIFY(id%PTRFAC)
          IF (lpok) THEN
            WRITE(lp,*) id%MYID,
     &      ': Allocation error for id%PTRFAC(', id%KEEP(28),')'
          ENDIF
          GOTO 110
        END IF
C       -----------------------------
C       Reserve temporary workspace :
C       IPOOL, PTRWB, ITLOC, PTRIST
C       PTRWB will be subdivided again
C       in routine ZMUMPS_FAC_B
C       -----------------------------
        ptrist = 1
        ptrwb  = ptrist + id%KEEP(28)
        itloc  = ptrwb  + 2 * id%KEEP(28)
C Fwd in facto: ITLOC of size id%N + id%KEEP(253)
        ipool  = itloc  + id%N + id%KEEP(253)
C
C       --------------------------------
C       NA(1) is an upperbound for LPOOL
C       --------------------------------
C       Structure of the pool:
C     ____________________________________________________
C    | Subtrees   |         | Top nodes           | 1 2 3 |
C     ----------------------------------------------------
        lpool = mumps_get_pool_length(id%NA(1), id%KEEP(1),id%KEEP8(1))
        ALLOCATE( iwk(  ipool + lpool - 1 ), stat = ierr )
        IF ( ierr .NE. 0 ) THEN
          id%INFO(1)=-13
          id%INFO(2)=ipool + lpool - 1
          IF (lpok) THEN
            WRITE(lp,*) id%MYID,
     &      ': Allocation error for IWK(',ipool+lpool-1,')'
          ENDIF
          GOTO 110
        END IF
        ALLOCATE(iwk8( 2 * id%KEEP(28)), stat = ierr)
        IF ( ierr .NE. 0 ) THEN
          id%INFO(1)=-13
          id%INFO(2)=2 * id%KEEP(28)
          IF (lpok) THEN
            WRITE(lp,*) id%MYID,
     &      ': Allocation error for IWKB(', 2*id%KEEP(28),')'
          ENDIF
          GOTO 110
        END IF
C
C  Return to SPMD
C
      ENDIF
C
 110  CONTINUE
C     ----------------
C     Broadcast errors
C     ----------------
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF ( id%INFO( 1 ) .LT. 0 ) GOTO 500
C
      IF ( i_am_slave )  THEN
C       Store size of receive buffers in ZMUMPS_LBUF module
        CALL zmumps_buf_dist_irecv_size( zmumps_lbufr_bytes )
        IF (prok) THEN
          WRITE( mp, 170 ) maxs, maxis, id%KEEP8(12), keep(15),
     &    id%KEEP8(26), id%KEEP8(27), id%KEEP8(11), keep(26), keep(27)
        ENDIF
      END IF
C     ===============================================================
C     Before calling the main driver, ZMUMPS_FAC_B,
C     some statistics should be initialized to 0,
C     even on the host node because they will be
C     used in REDUCE operations afterwards.
C     --------------------------------------------
C     Size of factors written. It will be set to POSFAC in
C     IC, otherwise we accumulate written factors in it.
      id%KEEP8(31)= 0_8
C     Size of factors under L0 will be returned
C     in id%KEEP8(64), not included in KEEP8(31))
C     Number of entries in factors
      id%KEEP8(10) = 0_8
C     KEEP8(8) will hold the volume of extra copies due to
C              in-place stacking in fac_mem_stack.F
      id%KEEP8(8)=0_8
      id%INFO(9:14)=0
      rinfo(2:3)=zero
      IF ( i_am_slave ) THEN
C       ------------------------------------
C       Call effective factorization routine
C       ------------------------------------
        IF ( keep(55) .eq. 0 ) THEN
          ldptrar = id%N
        ELSE
          ldptrar = id%NELT + 1
        END IF
        IF ( id%KEEP(55) .NE. 0 ) THEN
          nelt_arg = id%NELT
        ELSE
C         ------------------------------
C         Use size 1 to avoid complaints
C         when using check bound options
C         ------------------------------
          nelt_arg = 1
        END IF
      ENDIF
      IF (i_am_slave) THEN
        IF (associated(id%L0_OMP_MAPPING))
     &    DEALLOCATE(id%L0_OMP_MAPPING)
        IF (keep(400) .GT. 0) THEN
          id%LL0_OMP_MAPPING = keep(28)
        ELSE
          id%LL0_OMP_MAPPING = 1
        ENDIF
        ALLOCATE(id%L0_OMP_MAPPING(id%LL0_OMP_MAPPING), stat=allocok)
        IF ( allocok > 0) THEN
            write(*,*) "Problem allocating L0_OMP_MAPPING",
     &              ierr, keep(28)
            GOTO 115
        ENDIF
        IF (associated(id%L0_OMP_FACTORS)) THEN
          CALL zmumps_free_l0_omp_factors(id%L0_OMP_FACTORS)
        ENDIF
        IF (keep(400) .GT. 0) THEN
          id%LL0_OMP_FACTORS = keep(400)
        ELSE
          id%LL0_OMP_FACTORS = 1
        ENDIF
        ALLOCATE(id%L0_OMP_FACTORS(id%LL0_OMP_FACTORS),stat = allocok)
        IF (allocok > 0) THEN
           id%INFO(1)=-7
           id%INFO(2)=nb_threads
           GOTO 111
        ENDIF
        CALL zmumps_init_l0_omp_factors(id%L0_OMP_FACTORS)
      ENDIF
 115  CONTINUE
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
      IF ( id%INFO( 1 ) .LT. 0 ) GOTO 500
C     Compute DKEEP(17) 
      avg_flops    =  rinfog(1)/(dble(id%NSLAVES))
      id%DKEEP(17) = max( id%DKEEP(18), avg_flops/dble(50) )
     &               
      IF (prok.AND.id%MYID.EQ.master) THEN
        IF (id%NSLAVES.LE.1) THEN
         WRITE(mp,'(/A,A,1PD10.3)')
     &' Start factorization with total', 
     &' estimated flops (RINFOG(1))                         = ', 
     &   rinfog(1)
        ELSE
         WRITE(mp,'(/A,A,1PD10.3,A,1PD10.3)')
     &' Start factorization with total', 
     &' estimated flops RINFOG(1) / Average per MPI proc    = ', 
     &      rinfog(1), ' / ', avg_flops
        ENDIF
      ENDIF
      IF (i_am_slave) THEN
C       IS/S pointers passed to ZMUMPS_FAC_B with
C       implicit interface through intermediate
C       structure S_IS_POINTERS. IS will be allocated
C       during ZMUMPS_FAC_B.
C       In case of L0OMP, id%IS and id%S are allocated during
C       ZMUMPS_FAC_B, and only after L0OMP nodes are processed,
C       in order to limit the global memory peak. 
        s_is_pointers%IW => id%IS; NULLIFY(id%IS)
        s_is_pointers%A  => id%S ; NULLIFY(id%S)
        CALL zmumps_fac_b(id%N,s_is_pointers,maxs,maxis,id%SYM_PERM(1),
     &  id%NA(1),id%LNA,id%NE_STEPS(1),id%ND_STEPS(1), id%FILS(1),
     &  id%STEP(1),id%FRERE_STEPS(1),id%DAD_STEPS(1),id%CANDIDATES(1,1),
     &  id%ISTEP_TO_INIV2(1),id%TAB_POS_IN_PERE(1,1), id%PTRAR(1),
     &  ldptrar,iwk(ptrist),id%PTLUST_S(1),id%PTRFAC(1),iwk(ptrwb),iwk8,
     &  iwk(itloc),rhs_mumps(1),iwk(ipool),lpool,cntl1,icntl(1),
     &  id%INFO(1), rinfo(1),keep(1),id%KEEP8(1),id%PROCNODE_STEPS(1),
     &  id%NSLAVES,id%COMM_NODES,id%MYID,id%MYID_NODES,bufr,zmumps_lbufr
     &  , zmumps_lbufr_bytes, zmumps_lbuf, id%INTARR(1),id%DBLARR(1),
     &  id%root, nelt_arg, id%FRTPTR(1), id%FRTELT(1),id%COMM_LOAD,
     &  id%ASS_IRECV, seuil, seuil_ldlt_niv2, id%MEM_DIST(0),
     &  id%DKEEP(1), id%PIVNUL_LIST(1), lpn_list, id%LRGROUPS(1)
     &  ,id%IPOOL_B_L0_OMP(1),id%LPOOL_B_L0_OMP,
     &  id%IPOOL_A_L0_OMP(1),id%LPOOL_A_L0_OMP,id%L_VIRT_L0_OMP,
     &  id%VIRT_L0_OMP(1), id%VIRT_L0_OMP_MAPPING(1),id%L_PHYS_L0_OMP,
     &  id%PHYS_L0_OMP(1), id%PERM_L0_OMP(1), id%PTR_LEAFS_L0_OMP(1),
     &  id%L0_OMP_MAPPING(1),id%LL0_OMP_MAPPING,
     &  id%THREAD_LA, id%L0_OMP_FACTORS(1), id%LL0_OMP_FACTORS,
     &  id%I4_L0_OMP(1,1), size(id%I4_L0_OMP,1), size(id%I4_L0_OMP,2),
     &  id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1), size(id%I8_L0_OMP,2)
     &     )
        id%IS => s_is_pointers%IW; NULLIFY(s_is_pointers%IW)
        id%S  => s_is_pointers%A ; NULLIFY(s_is_pointers%A)
C
C       ------------------------------
C       Deallocate temporary workspace
C       ------------------------------
        DEALLOCATE( iwk  )
        DEALLOCATE( iwk8 )
      ENDIF
C     ---------------------------------
C     Free some workspace corresponding
C     to the original matrix in
C     arrowhead or elemental format.
C                  -----
C     Note : INTARR was not allocated
C     during factorization in the case
C     of an assembled matrix.
C     ---------------------------------
        IF ( keep(55) .eq. 0 ) THEN
C
C         ----------------
C         Assembled matrix
C         ----------------
          IF (associated( id%DBLARR)) THEN
            DEALLOCATE(id%DBLARR)
            NULLIFY(id%DBLARR)
          ENDIF
C
        ELSE
C
C         ----------------
C         Elemental matrix
C         ----------------
          IF (associated(id%INTARR)) THEN
            DEALLOCATE( id%INTARR)
            NULLIFY( id%INTARR )
          ENDIF
C         ------------------------------------
C         For the master from an hybrid host
C         execution without scaling, then real
C         values have not been copied !
C         -------------------------------------
          IF (      id%MYID_NODES .eq. master
     &      .AND.   keep(46)   .eq. 1
     &      .AND.   keep(52)   .eq. 0 ) THEN
            NULLIFY( id%DBLARR )
          ELSE
            IF (associated( id%DBLARR)) THEN
              DEALLOCATE(id%DBLARR)
              NULLIFY(id%DBLARR)
            ENDIF
          END IF
        END IF
C     Memroy statistics
C     -----------------------------------
C     If QR (Keep(19)) is not zero, and if
C     the host does not have the information
C     (ie is not slave), send information
C     computed on the slaves during facto
C     to the host.
C     -----------------------------------
      IF ( keep(19) .NE. 0 ) THEN
        IF ( keep(46) .NE. 1 ) THEN
C         Host was not working during facto_root
C         Send him the information
          IF ( id%MYID .eq. master ) THEN
            CALL mpi_recv( keep(17), 1, mpi_integer, 1, defic_tag,
     &                   id%COMM, status, ierr )
            CALL mpi_recv( keep(143), 1, mpi_integer, 1, defic_tag,
     &                   id%COMM, status, ierr )
          ELSE IF ( id%MYID .EQ. 1 ) THEN
            CALL mpi_send( keep(17), 1, mpi_integer, 0, defic_tag,
     &                   id%COMM, ierr )
            CALL mpi_send( keep(143), 1, mpi_integer, 0, defic_tag,
     &                   id%COMM, ierr )
          END IF
        END IF
      END IF
C     --------------------------------
C     Deallocate communication buffers
C     They will be reallocated
C     in the solve.
C     --------------------------------
      IF (allocated(bufr)) DEALLOCATE(bufr)
      CALL zmumps_buf_deall_small_buf( ierr )
C//PIV
      IF (keep(219).NE.0) THEN
      CALL zmumps_buf_deall_max_array()
      ENDIF
C
C     Check for errors.
C     After ZMUMPS_FAC_B every slave is aware of an error. 
C     If master is included in computations, the call below should
C     not be necessary.
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &                    id%COMM, id%MYID )
C
      CALL zmumps_extract_schur_redrhs(id)
      IF (keep(201) .GT. 0) THEN
         IF ((keep(201).EQ.1) .OR. (keep(201).EQ.2)) THEN
            IF ( i_am_slave ) THEN
               CALL zmumps_ooc_clean_pending(ierr)
               IF(ierr.LT.0)THEN
                  id%INFO(1)=ierr
                  id%INFO(2)=0
               ENDIF
            ENDIF
            CALL mumps_propinfo( id%ICNTL(1), id%INFO(1),
     &           id%COMM, id%MYID )
C             We want to collect statistics even in case of 
C             error to understand if it is due to numerical
C             issues
CC            IF ( id%INFO(1) < 0 ) GOTO 500
         END IF
      END IF
      IF (id%MYID.EQ.master) THEN
        CALL mumps_secfin(time)
        id%DKEEP(94)=time
        IF (keep(400).GT.0) THEN
C         Facto time above L0_OMP = total time - facto time under L0_OMP
          id%DKEEP(96)=id%DKEEP(94)-id%DKEEP(95)
        ENDIF
      ENDIF
C     =====================================================================
C     COMPUTE MEMORY ALLOCATED BY MUMPS, INFO(16)
C     ---------------------------------------------
      mem_eff_allocated = .true.
      CALL zmumps_max_mem( id%KEEP(1),id%KEEP8(1),
     &     id%MYID, n, id%NELT, id%NA(1), id%LNA, id%KEEP8(28),
     &     id%KEEP8(30),
     &     id%NSLAVES, total_mbytes, .true., id%KEEP(201),
     &     blr_strat, .true., total_bytes,
     &     idummy, bdummy, mem_eff_allocated
     &     , .false.   ! UNDER_L0_OMP
     &    , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &     size(id%I8_L0_OMP,2)
     &     )
      IF (keep(400) .GT. 0 ) THEN ! L0 activated
        CALL zmumps_max_mem( id%KEEP(1),id%KEEP8(1),
     &     id%MYID, n, id%NELT, id%NA(1), id%LNA, id%KEEP8(28),
     &     id%KEEP8(30),
     &     id%NSLAVES, total_mbytes_under_l0, .true., id%KEEP(201),
     &     blr_strat, .true., total_bytes_under_l0,
     &     idummy, bdummy, mem_eff_allocated
     &     , .true.   ! UNDER_L0_OMP
     &     , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &       size(id%I8_L0_OMP,2)
     &        )
        total_mbytes = max(total_mbytes,total_mbytes_under_l0) 
        total_bytes  = max(total_bytes, total_bytes_under_l0)
      ENDIF
       IF (id%KEEP8(24).NE.0) THEN
C         WK_USER is not part of memory allocated by MUMPS
C         and is not counted, id%KEEP8(23) should be zero
          id%INFO(16) = total_mbytes
       ELSE
C         Note that even for the case of ICNTL(23)>0 
C         we report here the memory effectively allocated
C         that can be smaller than ICNTL(23) !
          id%INFO(16) = total_mbytes
       ENDIF
C       ----------------------------------------------------
C       Centralize memory statistics on the host
C       id%INFOG(18) = size of mem in Mbytes for facto,
C                   for the processor using largest memory
C       id%INFOG(19) = size of mem in Mbytes for facto,
C                   sum over all processors
C       ----------------------------------------------------
        CALL mumps_mem_centralize( id%MYID, id%COMM,
     &                           id%INFO(16), id%INFOG(18), irank )
        CALL zmumps_print_allocated_mem( prok, prokg, print_maxavg,
     &       mp, mpg, id%INFO(16), id%INFOG(18), id%INFOG(19),
     &       id%NSLAVES, irank,
     &       id%KEEP(1) )
C     If WK_USER is provided, this excludes WK_USER
      IF (prok ) THEN
          WRITE(mp,'(A,I12) ')
     & ' ** Eff. min. Space MBYTES for facto             (INFO(16)):',
     &                total_mbytes
      ENDIF
C     ========================(INFO(16) RELATED)======================
C     ---------------------------------------
C     COMPUTE EFFECTIVE MEMORY USED INFO(22)
C     ---------------------------------------
      perlu_on = .true.
      mem_eff_allocated = .false.
      CALL zmumps_max_mem( id%KEEP(1),id%KEEP8(1),
     &     id%MYID, n, id%NELT, id%NA(1), id%LNA, id%KEEP8(28),
     &     id%KEEP8(30),
     &     id%NSLAVES, total_mbytes, .true., id%KEEP(201),
     &     blr_strat, perlu_on, total_bytes,
     &     idummy, bdummy, mem_eff_allocated
     &     , .false.   ! UNDER_L0_OMP
     &     , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &       size(id%I8_L0_OMP,2)
     &     )
      IF (keep(400) .GT. 0 ) THEN ! L0 activated
        CALL zmumps_max_mem( id%KEEP(1),id%KEEP8(1),
     &    id%MYID, n, id%NELT, id%NA(1), id%LNA, id%KEEP8(28),
     &    id%KEEP8(30),
     &    id%NSLAVES, total_mbytes_under_l0, .true., id%KEEP(201),
     &    blr_strat, perlu_on, total_bytes_under_l0,
     &    idummy, bdummy, mem_eff_allocated
     &    , .true. ! UNDER_L0_OMP
     &    , id%I8_L0_OMP(1,1), size(id%I8_L0_OMP,1),
     &      size(id%I8_L0_OMP,2)
     &  )
        total_mbytes = max(total_mbytes,total_mbytes_under_l0) 
        total_bytes  = max(total_bytes, total_bytes_under_l0)
      ENDIF
C     -- TOTAL_BYTES and TOTAL_MBYTES includes both static 
C     -- (MAXS) and BLR structures computed as the SUM of the PEAKS 
C     -- (KEEP8(67) + KEEP8(70))
      id%KEEP8(7) = total_bytes
C     -- INFO(22) holds the effective space (in Mbytes) used by MUMPS
C     -- (it includes part of WK_USER used if provided by user)
      id%INFO(22) = total_mbytes
C     ----------------------------------------------------
C     Centralize memory statistics on the host
C       INFOG(21) = size of effective mem (Mbytes) for facto,
C                   for the processor using largest memory
C       INFOG(22) = size of effective mem (Mbytes) for facto,
C                   sum over all processors
C     ----------------------------------------------------
      CALL mumps_mem_centralize( id%MYID, id%COMM,
     &                    id%INFO(22), id%INFOG(21), irank )
      IF ( prokg ) THEN
       IF (print_maxavg) THEN
        WRITE( mpg,'(A,I12) ')
     &  ' ** Memory effectively used, max in  Mbytes     (INFOG(21)):',
     &  id%INFOG(21)
       ENDIF
       WRITE( mpg,'(A,I12) ')
     &  ' ** Memory effectively used, total in Mbytes    (INFOG(22)):',
     &  id%INFOG(22)
      END IF
      sum_info22_this_node=0
      CALL mpi_reduce( id%INFO(22), sum_info22_this_node, 1,
     &                 mpi_integer,
     &                 mpi_sum, 0, id%KEEP(411), ierr )
      CALL mpi_reduce( sum_info22_this_node, max_sum_info22_this_node,
     &                 1, mpi_integer, mpi_max, 0, id%COMM, ierr )
      IF (prokg .AND. print_nodeinfo) THEN
        WRITE(mpg,'(A,I12)')
     &  ' ** Max. effective space per compute node, in MBytes       :',
     &  max_sum_info22_this_node
      ENDIF
C
      IF (i_am_slave) THEN
       k67 = id%KEEP8(67)
       k68 = id%KEEP8(68)
       k70 = id%KEEP8(70)
       k74 = id%KEEP8(74)
       k75 = id%KEEP8(75)
      ELSE
       k67 = 0_8
       k68 = 0_8
       k70 = 0_8
       k74 = 0_8
       k75 = 0_8
      ENDIF
C     -- Save the number of entries effectively used
C        in main working array S
      CALL mumps_seti8toi4(k67,id%INFO(21))
C
C
      IF (keep(400) .GT.0 ) THEN
        IF (.NOT. i_am_slave) THEN
          id%DKEEP(95) = 0.0d0
          id%DKEEP(16) = 0.0d0
        ENDIF
        IF (id%NPROCS .GT. 1) THEN
C         Compute average and max (across MPI's)
          CALL mpi_reduce(id%DKEEP(95), tmptime, 1,
     &    mpi_double_precision, mpi_sum, master, id%COMM, ierr)
          IF (id%MYID.EQ.master) timeavg = tmptime
          CALL mpi_reduce(id%DKEEP(16), tmpflop, 1,
     &    mpi_double_precision, mpi_sum, master, id%COMM, ierr)
          IF (id%MYID.EQ.master) flopavg = tmpflop
          IF (id%MYID.EQ.master) THEN
            timeavg = timeavg / id%NSLAVES
            flopavg = flopavg / id%NSLAVES
          ENDIF
          CALL mpi_reduce(id%DKEEP(95), timemax, 1, 
     &    mpi_double_precision, mpi_max, master, id%COMM, ierr)
          CALL mpi_reduce(id%DKEEP(16), flopmax, 1, 
     &    mpi_double_precision, mpi_max, master, id%COMM, ierr)
C         (PROKG may only be true on master)
          IF ( prokg ) THEN
            WRITE(mpg,190) flopavg, flopmax
            WRITE(mpg,188) timeavg, timemax
          ENDIF
        ELSE
C          Print DKEEP(95) directly without reduction
           IF ( prokg ) THEN
             WRITE(mpg,189) id%DKEEP(16)
             WRITE(mpg,187) id%DKEEP(95)
           ENDIF
        ENDIF
      ENDIF
      IF ( prokg ) THEN
          IF (id%INFO(1) .GE.0) THEN
            WRITE(mpg,180) id%DKEEP(94)
          ELSE
            WRITE(mpg,185) id%DKEEP(94)
          ENDIF
      ENDIF
C
C  Sum RINFO(2) : total number of flops for assemblies
C  Sum RINFO(3) : total number of flops for eliminations
C     Initialize RINFO(4) in case BLR was not activated
      rinfo(4) = rinfo(3)
C 
C  Should work even if the master does some work
C
      CALL mpi_reduce( rinfo(2), rinfog(2), 2,
     &                 mpi_double_precision,
     &                 mpi_sum, master, id%COMM, ierr)
C     Reduce needed to dimension small working array
C     on all procs during ZMUMPS_GATHER_SOLUTION
      keep(247) = 0
      CALL mpi_reduce( keep(246), keep(247), 1, mpi_integer, 
     &                 mpi_max, master, id%COMM, ierr)
C
C     Reduce compression times: get max compression times
      CALL mpi_reduce( id%DKEEP(97), id%DKEEP(98), 1,
     &     mpi_double_precision,
     &     mpi_max, master, id%COMM, ierr)
C
      CALL mpi_reduce( rinfo(2), rinfog(2), 2,
     &                 mpi_double_precision,
     &                 mpi_sum, master, id%COMM, ierr)
      CALL mumps_reducei8( id%KEEP8(31)+id%KEEP8(64),id%KEEP8(6),
     &                     mpi_sum, master, id%COMM )
C
      IF (id%MYID.EQ.0) THEN
C      In MegaBytes
       rinfog(16) = dble(id%KEEP8(6)*int(keep(35),8))/dble(1d6)
       IF (keep(201).LE.0) THEN
        rinfog(16) = zero
       ENDIF
      ENDIF
      CALL mumps_reducei8( id%KEEP8(48),id%KEEP8(148), mpi_sum,
     &                     master, id%COMM )
      CALL mumps_seti8toi4(id%KEEP8(148), infog(9))
C
      CALL mpi_reduce( int(id%INFO(10),8), id%KEEP8(128),
     &                 1, mpi_integer8,
     &                 mpi_sum, master, id%COMM, ierr)
      IF (id%MYID.EQ.master) THEN
        CALL mumps_seti8toi4(id%KEEP8(128), id%INFOG(10))
      ENDIF
C     Use MPI_MAX for this one to get largest front size
      CALL mpi_allreduce( id%INFO(11), infog(11), 1, mpi_integer,
     &                 mpi_max, id%COMM, ierr)
C     make maximum effective frontal size available on all procs
C     for solve phase
C     (Note that INFO(11) includes root size on root master)
      keep(133) = infog(11)
      CALL mpi_reduce( id%INFO(12), infog(12), 3, mpi_integer,
     &                 mpi_sum, master, id%COMM, ierr)
      CALL mpi_reduce( keep(103), infog(25), 1, mpi_integer,
     &                 mpi_sum, master, id%COMM, ierr)
      keep(229) = infog(25)
      CALL mpi_reduce( keep(105), infog(25), 1, mpi_integer,
     &                 mpi_sum, master, id%COMM, ierr)
      keep(230) = infog(25)
C
      id%INFO(25) = keep(98)
      CALL mpi_allreduce( id%INFO(25), infog(25), 1, mpi_integer,
     &                 mpi_sum, id%COMM, ierr)
C     Extra copies due to in-place stacking
      CALL mumps_reducei8( id%KEEP8(8), id%KEEP8(108), mpi_sum,
     &                     master, id%COMM )
C     Entries in factors
      CALL mumps_seti8toi4(id%KEEP8(10), id%INFO(27))
      CALL mumps_reducei8( id%KEEP8(10),id%KEEP8(110), mpi_sum,
     &                     master, id%COMM )
      CALL mumps_seti8toi4(id%KEEP8(110), infog(29))
C     Initialize INFO(28)/INFOG(35) in case BLR not activated
      id%INFO(28)  = id%INFO(27)
      infog(35)    = infog(29)
C     ==============================
C     LOW-RANK
C     ==============================
      IF ( keep(486) .NE. 0 ) THEN  !LR is activated
C Compute and Save local amount of flops in case of BLR
            rinfo(4) = dble(flop_frfronts + flop_facto_fr - flop_lrgain
     &                 + flop_compress  + flop_frfronts)
C
C Compute and Save local number of entries in compressed factors
C 
            itmp8 =  id%KEEP8(10) - int(mry_lu_lrgain,8)
            CALL mumps_seti8toi4( itmp8, id%INFO(28))
C
            CALL mpi_reduce( mry_lu_lrgain, tmp_mry_lu_lrgain
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( mry_lu_fr, tmp_mry_lu_fr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( mry_cb_fr, tmp_mry_cb_fr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( mry_cb_lrgain, tmp_mry_cb_lrgain
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_lrgain, tmp_flop_lrgain
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_trsm_fr, tmp_flop_trsm_fr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_trsm_lr, tmp_flop_trsm_lr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_update_fr, tmp_flop_update_fr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_update_lr, tmp_flop_update_lr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_frswap_compress,
     &                       tmp_flop_frswap_compress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_midblk_compress,
     &                       tmp_flop_midblk_compress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_update_lrlr3, tmp_flop_update_lrlr3
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce(flop_accum_compress, tmp_flop_accum_compress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_trsm, tmp_flop_trsm
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_panel, tmp_flop_panel
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_frfronts, tmp_flop_frfronts
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_compress, tmp_flop_compress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_decompress, tmp_flop_decompress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_cb_compress, tmp_flop_cb_compress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_cb_decompress,tmp_flop_cb_decompress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( flop_facto_fr, tmp_flop_facto_fr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( cnt_nodes,tmp_cnt_nodes
     &                      , 1, mpi_integer,
     &                       mpi_sum, master, id%COMM, ierr)
            IF (id%NPROCS.GT.1) THEN
              flop_facto_lr = flop_facto_fr - flop_lrgain
     &                          + flop_compress + flop_frfronts
              CALL mpi_reduce( flop_facto_lr, avg_flop_facto_lr
     &                        , 1, mpi_double_precision,
     &                         mpi_sum, master, id%COMM, ierr)
              IF (id%MYID.EQ.master) THEN
                avg_flop_facto_lr = avg_flop_facto_lr/id%NPROCS
              ENDIF
              CALL mpi_reduce( flop_facto_lr, min_flop_facto_lr
     &                        , 1, mpi_double_precision,
     &                         mpi_min, master, id%COMM, ierr)
              CALL mpi_reduce( flop_facto_lr, max_flop_facto_lr
     &                        , 1, mpi_double_precision,
     &                         mpi_max, master, id%COMM, ierr)
            ENDIF ! NPROCS > 1
            CALL mpi_reduce( time_update, tmp_time_update
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_update_lrlr1, tmp_time_update_lrlr1
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_update_lrlr2, tmp_time_update_lrlr2
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_update_lrlr3, tmp_time_update_lrlr3
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_update_frlr, tmp_time_update_frlr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_update_frfr, tmp_time_update_frfr
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_diagcopy, tmp_time_diagcopy
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_compress,tmp_time_compress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_midblk_compress, 
     &                       tmp_time_midblk_compress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_frswap_compress,
     &                       tmp_time_frswap_compress
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_cb_compress, tmp_time_cb_compress
     &                      , 1, mpi_double_precision, 
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_decomp, tmp_time_decomp
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_decomp_ucfs, tmp_time_decomp_ucfs
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_decomp_asm1, tmp_time_decomp_asm1
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce(time_decomp_locasm2, tmp_time_decomp_locasm2
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce(time_decomp_maplig1, tmp_time_decomp_maplig1
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_decomp_asms2s, tmp_time_decomp_asms2s
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_decomp_asms2m, tmp_time_decomp_asms2m
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_panel, tmp_time_panel
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_fac_i, tmp_time_fac_i
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_fac_mq, tmp_time_fac_mq
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_fac_sq, tmp_time_fac_sq
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_lrtrsm, tmp_time_lrtrsm
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_frtrsm, tmp_time_frtrsm
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_frfronts, tmp_time_frfronts
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
            CALL mpi_reduce( time_lr_module, tmp_time_lr_module
     &                      , 1, mpi_double_precision,
     &                       mpi_sum, master, id%COMM, ierr)
        IF (id%MYID.EQ.master) THEN
         IF (id%NPROCS.GT.1) THEN
C rename the stat variable so that COMPUTE_GLOBAL_GAINS can work for any
C number of procs
            mry_lu_fr            = tmp_mry_lu_fr
            mry_lu_lrgain        = tmp_mry_lu_lrgain
            mry_cb_fr            = tmp_mry_cb_fr
            mry_cb_lrgain        = tmp_mry_cb_lrgain
            flop_lrgain          = tmp_flop_lrgain
            flop_panel           = tmp_flop_panel
            flop_trsm            = tmp_flop_trsm
            flop_trsm_fr         = tmp_flop_trsm_fr
            flop_trsm_lr         = tmp_flop_trsm_lr
            flop_update_fr       = tmp_flop_update_fr
            flop_update_lr       = tmp_flop_update_lr
            flop_update_lrlr3    = tmp_flop_update_lrlr3
            flop_compress        = tmp_flop_compress
            flop_midblk_compress = tmp_flop_midblk_compress
            flop_frswap_compress = tmp_flop_frswap_compress
            flop_accum_compress  = tmp_flop_accum_compress
            flop_cb_compress     = tmp_flop_cb_compress
            flop_decompress      = tmp_flop_decompress
            flop_cb_decompress   = tmp_flop_cb_decompress
            flop_frfronts        = tmp_flop_frfronts
            flop_facto_fr        = tmp_flop_facto_fr
            cnt_nodes            = tmp_cnt_nodes
            time_update          = tmp_time_update         /id%NPROCS
            time_update_lrlr1    = tmp_time_update_lrlr1   /id%NPROCS
            time_update_lrlr2    = tmp_time_update_lrlr2   /id%NPROCS
            time_update_lrlr3    = tmp_time_update_lrlr3   /id%NPROCS
            time_update_frlr     = tmp_time_update_frlr    /id%NPROCS
            time_update_frfr     = tmp_time_update_frfr    /id%NPROCS
            time_compress        = tmp_time_compress       /id%NPROCS
            time_midblk_compress = tmp_time_midblk_compress/id%NPROCS
            time_frswap_compress = tmp_time_frswap_compress/id%NPROCS
            time_diagcopy        = tmp_time_diagcopy       /id%NPROCS
            time_cb_compress     = tmp_time_cb_compress    /id%NPROCS
            time_panel           = tmp_time_panel          /id%NPROCS
            time_fac_i           = tmp_time_fac_i          /id%NPROCS
            time_fac_mq          = tmp_time_fac_mq         /id%NPROCS
            time_fac_sq          = tmp_time_fac_sq         /id%NPROCS
            time_lrtrsm          = tmp_time_lrtrsm         /id%NPROCS
            time_frtrsm          = tmp_time_frtrsm         /id%NPROCS
            time_frfronts        = tmp_time_frfronts       /id%NPROCS
            time_lr_module       = tmp_time_lr_module      /id%NPROCS
            time_decomp          = tmp_time_decomp         /id%NPROCS
            time_decomp_ucfs     = tmp_time_decomp_ucfs    /id%NPROCS
            time_decomp_asm1     = tmp_time_decomp_asm1    /id%NPROCS
            time_decomp_locasm2  = tmp_time_decomp_locasm2 /id%NPROCS
            time_decomp_maplig1  = tmp_time_decomp_maplig1 /id%NPROCS
            time_decomp_asms2s   = tmp_time_decomp_asms2s  /id%NPROCS
            time_decomp_asms2m   = tmp_time_decomp_asms2m  /id%NPROCS
         ENDIF
         CALL compute_global_gains(id%KEEP8(110),id%RINFOG(3),
     &        id%KEEP8(49), prokg, mpg)
C         Number of entries in factor  INFOG(35) in 
C         compressed form is updated as long as 
C         BLR is activated,  this independently of the 
C         fact that factors are saved in LR.
          CALL mumps_seti8toi4(id%KEEP8(49),  id%INFOG(35))
          frontwise = 0
C         WRITE gains also compute stats stored in DKEEP array
          IF (lpok) THEN
            IF (cntl(7) < 0.0d0) THEN
C           Warning : using negative values is an experimental and 
C            non recommended setting.
             WRITE(lp,'(/A/,A/,A/,A,A)') 
     &  ' WARNING in BLR input setting',
     &  '          CNTL(7) < 0 is experimental: ',
     &  '          RRQR precision = |CNTL(7| x ||A_pre||, ',
     &  '          where A_pre is the preprocessed matrix as defined',
     &  ' in the Users guide '
            ENDIF
          ENDIF
          CALL saveandwrite_gains(frontwise,
     &                keep(489), id%DKEEP, n,  id%ICNTL(36),
     &                keep(487), keep(488), keep(490),
     &                keep(491), keep(50), keep(486), 
     &                keep(472), keep(475), keep(478), keep(480), 
     &                keep(481),
     &                keep(483), keep(484), 
     &                id%KEEP8(110), id%KEEP8(49),
     &                keep(28), id%NPROCS, mpg, prokg)
C           flops when BLR activated
            rinfog(14) = id%DKEEP(56)
          ELSE
            rinfog(14) = 0.0d00
          ENDIF
      ENDIF
C     ==============================
C     NULL PIVOTS AND RANK-REVEALING
C     ==============================
      IF(keep(110) .EQ. 1) THEN
C        -- make available to users the local number of null pivots detected 
C        -- with ICNTL(24) = 1.
         id%INFO(18) = keep(109)
         CALL mpi_allreduce( keep(109), keep(112), 1, mpi_integer,
     &        mpi_sum, id%COMM, ierr)
      ELSE
         id%INFO(18)  = 0
         keep(109) = 0
         keep(112) = 0
      ENDIF
      IF (id%MYID.EQ.master) THEN
C      INFOG(28) deficiency resulting from ICNTL(24) and ICNTL(56).
       infog(28)=keep(112)+keep(17)
      ENDIF
C     ========================================
C     We now provide to the host the part of
C     PIVNUL_LIST resulting from the processing
C     of the root node and we update id%INFO(18)
C     on the processor holding the root to
C     include null pivots relative to the root
C     ========================================
      IF (keep(17) .NE. 0) THEN
        IF (id%MYID .EQ. id_root) THEN
C         Include in id%INFO(18) null pivots resulting
C         from deficiency on the root. In this way,
C         the sum of all id%INFO(18) is equal to INFOG(28).
          id%INFO(18)=id%INFO(18)+keep(17)
        ENDIF
        IF (id_root .EQ. master) THEN
          IF (id%MYID.EQ.master) THEN
C           --------------------------------------------------
C           Null pivots of root have been stored in
C           PIVNUL_LIST(KEEP(109)+1:KEEP(109)+KEEP(17).
C           Shift them at the end of the list because:
C           * this is what we need to build the null space
C           * we would otherwise overwrite them on the host
C             when gathering null pivots from other processors
C           --------------------------------------------------
            DO i= keep(17), 1, -1
c           DO I=1, KEEP(17) % incorrect 
C                            when KEEP(112) < KEEP(109)+ KEEP(17)
              id%PIVNUL_LIST(keep(112)+i)=id%PIVNUL_LIST(keep(109)+i)
            ENDDO
          ENDIF
        ELSE
C         ---------------------------------
C         Null pivots of root must be sent
C         from the processor responsible of
C         the root to the host (or MASTER).
C         ---------------------------------
          IF (id%MYID .EQ. id_root) THEN
            CALL mpi_send(id%PIVNUL_LIST(keep(109)+1), keep(17),
     &                    mpi_integer, master, zero_piv,
     &                    id%COMM, ierr)
          ELSE IF (id%MYID .EQ. master) THEN
            CALL mpi_recv(id%PIVNUL_LIST(keep(112)+1), keep(17),
     &                    mpi_integer, id_root, zero_piv,
     &                    id%COMM, status, ierr )
          ENDIF
        ENDIF
      ENDIF
C     ===========================
C     gather zero pivots indices
C     on the host node
C     ===========================
C     In case of non working host, the following code also
C     works considering that KEEP(109) is equal to 0 on
C     the non-working host
      IF(keep(110) .EQ. 1) THEN
         ALLOCATE(itmp2(id%NPROCS),stat = ierr )  ! deallocated in 490
         IF ( ierr .GT. 0 ) THEN
            id%INFO(1)=-13
            id%INFO(2)=id%NPROCS
         END IF
         CALL mumps_propinfo( icntl(1), id%INFO(1),
     &     id%COMM, id%MYID )
         IF (id%INFO(1).LT.0) GOTO 490
         CALL mpi_gather ( keep(109),1, mpi_integer, 
     &        itmp2(1), 1, mpi_integer, 
     &        master, id%COMM, ierr)
         IF(id%MYID .EQ. master) THEN
            posbuf = itmp2(1)+1
C           First null pivot of master is in
C           position 1 of global list
            keep(220)=1
            DO i = 1,id%NPROCS-1
               CALL mpi_recv(id%PIVNUL_LIST(posbuf), itmp2(i+1), 
     &              mpi_integer,i, 
     &              zero_piv, id%COMM, status, ierr)
C              Send position POSBUF of first null pivot of proc I
C              in global list. Will allow to quickly identify during
C              the solve step if one is concerned by a global position
C              K, 0 <= K <= INFOG(28).
               CALL mpi_send(posbuf, 1, mpi_integer, i, zero_piv,
     &              id%COMM, ierr)
               posbuf = posbuf + itmp2(i+1)
            ENDDO
         ELSE
            CALL mpi_send( id%PIVNUL_LIST(1), keep(109), mpi_integer,
     &           master,zero_piv, id%COMM, ierr)
            CALL mpi_recv( keep(220), 1, mpi_integer, master, zero_piv,
     &           id%COMM, status, ierr )
         ENDIF
      ENDIF
C     =====================================
C     Statistics relative to min/max pivots
C     =====================================
      CALL mpi_reduce( id%DKEEP(19), rinfog(19), 1, 
     &                 mpi_double_precision,
     &                 mpi_min, master, id%COMM, ierr )
      CALL mpi_reduce( id%DKEEP(20), rinfog(20), 1, 
     &                 mpi_double_precision,
     &                 mpi_min, master, id%COMM, ierr )
      CALL mpi_reduce( id%DKEEP(21), rinfog(21), 1, 
     &                 mpi_double_precision,
     &                 mpi_max, master, id%COMM, ierr )
C     =========================================
C     Centralized number of swaps for pivoting
C     =========================================
      CALL mpi_reduce( id%KEEP8(80), itemp8, 1, mpi_integer8,
     &                 mpi_sum, master, id%COMM, ierr )
      IF (id%MYID .EQ. master) THEN
        CALL mumps_seti8toi4(itemp8,id%INFOG(48))
      ENDIF
C     ==========================================
C     Centralized largest increase of panel size
C     ==========================================
      CALL mpi_reduce( id%KEEP(425), id%INFOG(49), 1, mpi_integer,
     &                 mpi_max, master, id%COMM, ierr )
C     =====================================
C     Statistics concerning the determinant
C     =====================================
C
C     1/ on the host better take into account null pivots if scaling:
C
C     Since null pivots are excluded from the computation
C     of the determinant, we also exclude the corresponding
C     scaling entries. Since those entries have already been
C     taken into account before the factorization, we multiply
C     the determinant on the host by the scaling values corresponding
C     to pivots in PIVNUL_LIST.
      IF (id%MYID.EQ.master .AND. lscal. and. keep(258).NE.0) THEN
        k = min(keep(143), keep(17))
        k = max(k, 0)
        DO i = 1, keep(112)+ k
c       DO I = 1, id%INFOG(28)  ! all null pivots + singular values
          CALL zmumps_updatedeter_scaling(
     &                            id%ROWSCA(id%PIVNUL_LIST(i)),
     &                            id%DKEEP(6), keep(259))
          CALL zmumps_updatedeter_scaling(
     &                            id%COLSCA(id%PIVNUL_LIST(i)),
     &                            id%DKEEP(6), keep(259))
        ENDDO
      ENDIF
C
C     2/ Swap signs depending on pivoting on each proc
C
      IF (keep(258).NE.0) THEN
C       Return the determinant in INFOG(34) and RINFOG(12/13)
        IF (keep(260).EQ.-1) THEN ! Local to each processor
          id%DKEEP(6)=-id%DKEEP(6)
          id%DKEEP(7)=-id%DKEEP(7)
        ENDIF
C
C       3/ Perform a reduction
C
        CALL zmumps_deter_reduction(
     &           id%COMM, id%DKEEP(6), keep(259),
     &           rinfog(12), infog(34), id%NPROCS)
C
C       4/ Swap sign if needed
C
        IF (id%KEEP(50).EQ.0 .AND. id%MYID.EQ. master) THEN
C         Modify sign of determinant according
C         to unsymmetric permutation (max-trans
C         of max-weighted matching)
          IF (id%KEEP(23).NE.0) THEN
            CALL zmumps_deter_sign_perm(
     &           rinfog(12), id%N,
C           id%STEP: used as workspace of size N still
C                    allocated on master; restored on exit
     &           id%STEP(1),
     &           id%UNS_PERM(1) )
C           Remark that RINFOG(12/13) are modified only
C           on the host but will be broadcast on exit
C           from MUMPS (see ZMUMPS_DRIVER)
          ENDIF
        ENDIF
      ENDIF
 490  IF (allocated(itmp2)) DEALLOCATE(itmp2)
      IF ( prokg ) THEN
C     -----------------------------
C     PRINT STATISTICS  (on master)
C     -----------------------------
          WRITE(mpg,99984) rinfog(2),rinfog(3),keep(52), 
     &                    id%KEEP8(148),
     &                    id%KEEP8(128), infog(11), id%KEEP8(110)
          IF (id%KEEP(50) == 0) THEN
            ! off diag pivots
            WRITE(mpg, 99985) infog(12)
          END IF
          IF (id%KEEP(50) .NE. 1) THEN
            ! delayed pivots
            WRITE(mpg, 99982) infog(13)
          END IF
          IF (keep(97) .NE. 0) THEN
            ! tiny pivots
            WRITE(mpg, '(A,D16.4)') 
     &     ' Effective static pivoting thresh., CNTL(4) =', seuil
            WRITE(mpg, 99986) infog(25)
          ENDIF
          IF (id%KEEP(50) == 2) THEN
            !number of 2x2 pivots in type 1 nodes
             WRITE(mpg, 99988) keep(229)
            !number of 2x2 pivots in type 2 nodes
             WRITE(mpg, 99989) keep(230)
          ENDIF
          !number of zero pivots
          IF (keep(110) .NE.0) THEN
              WRITE(mpg, 99991) keep(112)
          ENDIF
          !Deficiency on root
          IF ( keep(19) .ne. 0 )
c         IF ( KEEP(17) .ne. 0 )
     &    WRITE(mpg, 99983) keep(17)
          !Total deficiency
          IF (keep(110).NE.0.OR.keep(19).NE.0)
c          IF (KEEP(110).NE.0.OR.KEEP(17).NE.0)
     &    WRITE(mpg, 99992) keep(17)+keep(112)
          ! Memory compress
          WRITE(mpg, 99981) infog(14)
          ! Extra copies due to ip stack in unsym case
          ! in core case (or OLD_OOC_PANEL)
          IF (id%KEEP8(108) .GT. 0_8) THEN
            WRITE(mpg, 99980) id%KEEP8(108)
          ENDIF
          IF  ((keep(60).NE.0) .AND. infog(25).GT.0) THEN
          !  Schur on and tiny pivots set in last level 
          ! before the Schur if KEEP(114)=0
           WRITE(mpg, '(A)') 
     & " ** Warning Static pivoting was necessary"
           WRITE(mpg, '(A)') 
     & " ** to factor interior variables with Schur ON"
          ENDIF
          IF (keep(258).NE.0) THEN
            WRITE(mpg,99978) rinfog(12)
            WRITE(mpg,99979) rinfog(13)
            WRITE(mpg,99977) infog(34)
          ENDIF
      END IF
* ==========================================
*
*  End of Factorization Phase
*
* ==========================================
C
C  Goto 500 is done when
C  LOAD_INIT
C  OOC_INIT_FACTO
C  MUMPS_FDM_INIT
#if ! defined(NO_FDM_DESCBAND)
C  MUMPS_FDBD_INIT
#endif
#if ! defined(NO_FDM_MAPROW)
C  MUMPS_FMRD_INIT
#endif
C  are all called.
C
 500  CONTINUE
C     Redo free DBLARR (as in end_driver.F)
C     in case an error occurred after allocating
C     DBLARR and before freeing it above.
      IF (id%KEEP(46).EQ.1 .AND.
     &    id%KEEP(55).NE.0 .AND.
     &    id%MYID.EQ.master .AND.
     &    id%KEEP(52) .EQ. 0) THEN
        NULLIFY(id%DBLARR)
      ELSE
        IF (associated(id%DBLARR)) THEN
          DEALLOCATE(id%DBLARR)
          NULLIFY(id%DBLARR)
        ENDIF
      ENDIF
#if ! defined(NO_FDM_DESCBAND)
      IF (i_am_slave) THEN
        CALL mumps_fdbd_end(id%INFO(1))  ! INFO(1): input only
      ENDIF
#endif
#if ! defined(NO_FDM_MAPROW)
      IF (i_am_slave) THEN
        CALL mumps_fmrd_end(id%INFO(1))  ! INFO(1): input only
      ENDIF
#endif
      IF (i_am_slave) THEN
C       Terminate BLR module except if it is still needed for solve
        IF ( 
     &       (
     &         (keep(486).EQ.2) 
     &       )
     &      .AND. id%INFO(1).GE.0
     &     ) THEN
C         Store pointer to BLR_ARRAY in MUMPS structure
C         (requires successful factorization otherwise module is freed)
          CALL zmumps_blr_mod_to_struc(id%BLRARRAY_ENCODING)
        ELSE
C         INFO(1) positive or negative
          CALL zmumps_blr_end_module(id%INFO(1), id%KEEP8, id%KEEP(34))
        ENDIF
      ENDIF
      IF (i_am_slave) THEN
        CALL mumps_fdm_end('A')
C       Terminate BLR module except if it is still needed for solve
        IF ( 
     &       (
     &         (keep(486).EQ.2) 
     &       )
     &      .AND. id%INFO(1).GE.0
     &     ) THEN
           CALL mumps_fdm_mod_to_struc('F', id%FDM_F_ENCODING,
     &        id%INFO(1))
           IF (.NOT. associated(id%FDM_F_ENCODING)) THEN
             WRITE(*,*) "Internal error 2 in ZMUMPS_FAC_DRIVER"
           ENDIF
        ELSE
           CALL mumps_fdm_end('F')
        ENDIF
      ENDIF
C
C  Goto 514 is done when an
C  error occurred in MUMPS_FDM_INIT
C  or (after FDM_INIT but before
C  OOC_INIT)
C
 514  CONTINUE
      IF ( i_am_slave ) THEN
         IF ((keep(201).EQ.1).OR.(keep(201).EQ.2)) THEN
            CALL zmumps_ooc_end_facto(id,ierr)
            IF (id%ASSOCIATED_OOC_FILES) THEN
              id%ASSOCIATED_OOC_FILES = .false.
            ENDIF
            IF (ierr.LT.0 .AND. id%INFO(1) .GE. 0) id%INFO(1) = ierr
         ENDIF
         IF (wk_user_provided) THEN
C     at the end of a phase S is always freed when WK_USER provided
            NULLIFY(id%S)
         ELSE IF (keep(201).NE.0) THEN
C           ----------------------------------------
C           In OOC or if KEEP(201).EQ.-1 we always
C           free S at end of factorization. As id%S
C           may be unassociated in case of error
C           during or before the allocation of id%S,
C           we only free S when it was associated.
C           ----------------------------------------
            IF (associated(id%S))  DEALLOCATE(id%S)
            NULLIFY(id%S)   ! in all cases
            id%KEEP8(23)=0_8
         ENDIF
      ELSE  ! host not working
         IF (wk_user_provided) THEN
C     at the end of a phase S is always freed when WK_USER provided
            NULLIFY(id%S)
         ELSE
            IF (associated(id%S))  DEALLOCATE(id%S)
            NULLIFY(id%S)   ! in all cases
            id%KEEP8(23)=0_8
         END IF
      END IF
C
C     Goto 513 is done in case of error where LOAD_INIT was
C     called but not OOC_INIT_FACTO.
 513  CONTINUE
      IF ( i_am_slave ) THEN
         CALL zmumps_load_end( id%INFO(1), id%NSLAVES, ierr )
         IF (ierr.LT.0 .AND. id%INFO(1) .GE. 0) id%INFO(1) = ierr
      ENDIF
      CALL mumps_propinfo( icntl(1), id%INFO(1),
     &     id%COMM, id%MYID )
C
C     Goto 517 is done when an error occurs when GPU initialization
C     has been performed but not LOAD_INIT or OOC_INIT_FACTO
C
 517  CONTINUE
C
C     Goto 530 is done when an error occurs before
C     the calls to GPU_INIT, LOAD_INIT and OOC_INIT_FACTO
 530  CONTINUE
C  Fwd in facto: free RHS_MUMPS in case
C  it was allocated.
      IF (rhs_mumps_allocated) DEALLOCATE(rhs_mumps)
      NULLIFY(rhs_mumps)
C
      id%KEEP8(26) = keep826_save
      RETURN
 120  FORMAT(/' Local redistrib: data local/sent           =',i16,i16)
 125  FORMAT(/' Redistrib: total data local/sent           =',i16,i16)
 130  FORMAT(//'****** FACTORIZATION STEP ********'/)
 160  FORMAT(
     & /' Elapsed time to reformat/distribute matrix =',f12.4)
 166  FORMAT(' Max difference from 1 after scaling the entries',
     &       ' for ONE-NORM (option 7/8)   =',d9.2)
 170  FORMAT(' STATISTICS PRIOR NUMERICAL FACTORIZATION ...'/
     &        ' Size of internal working array S           =',i16/
     &        ' Size of internal working array IS          =',i16/
     &        ' Minimum (ICNTL(14)=0) size of S            =',i16/
     &        ' Minimum (ICNTL(14)=0) size of IS           =',i16/
     &        ' Real space for original matrix             =',i16/
     &        ' Integer space for original matrix          =',i16/
     &        ' INFO(3) Real space for factors (estimated) =',i16/ 
     &        ' INFO(4) Integer space for factors (estim.) =',i16/
     &        ' Maximum frontal size (estimated)           =',i16)
 172  FORMAT(' GLOBAL STATISTICS PRIOR NUMERICAL FACTORIZATION ...'/
     &        ' Number of working processes                =',i16/
     &        ' ICNTL(22) Out-of-core option               =',i16/
     &        ' ICNTL(35) BLR activation (eff. choice)     =',i16/
     &        ' ICNTL(14) Memory relaxation                =',i16/
     &        ' INFOG(3) Real space for factors (estimated)=',i16/
     &        ' INFOG(4) Integer space for factors (estim.)=',i16/
     &        ' Maximum frontal size (estimated)           =',i16/
     &        ' Number of nodes in the tree                =',i16/
     &        ' ICNTL(23) Memory allowed (value on host)   =',i16/
     &        '           Sum over all procs               =',i16/
     &        ' Memory provided by user, sum of LWK_USER   =',i16/
     &        ' Effective threshold for pivoting, CNTL(1)  =',d16.4)
 173  FORMAT( ' Perform forward during facto, NRHS         =',i16)
 174  FORMAT( ' KEEP(268) Relaxed pivoting effective value =',i16)
 180  FORMAT(/' Elapsed time for factorization             =',f12.4)
 185  FORMAT(/' Elapsed time for (failed) factorization    =',f12.4)
 187  FORMAT( ' Elapsed time under L0                      =',f12.4)
 188  FORMAT( ' Elapsed time under L0 (avg/max across MPI) =',
     &          f12.4,f12.4)
 189  FORMAT(/' Flops under L0 layer                       =',1pd12.3)
 190  FORMAT(/' Flops under L0 Layer  (avg/max across MPI) =',
     &          1pd12.3,1pd12.3)
99977 FORMAT( ' infog(34)  determinant(base 2 exponent)   =',I16)
99978 FORMAT( ' rinfog(12) determinant(real part)         =',F16.8)
99979 FORMAT( ' rinfog(12) determinant(imaginary part)    =',F16.8)
99980 FORMAT( ' extra copies due to in-place stacking      =',I16)
99981 FORMAT( ' infog(14)  number of memory compress       =',I16)
99982 FORMAT( ' infog(13)  number of delayed pivots        =',I16)
99983 FORMAT( ' nb of singularities detected by icntl(56)  =',I16)
99991 FORMAT( ' nb of null pivots detected by icntl(24)    =',I16)
99992 FORMAT( ' infog(28)  estimated deficiency            =',I16)
99984 FORMAT(/'leaving factorization with ...'/
     &        ' rinfog(2)  operations in node assembly     =',1PD10.3/
     &        ' ------(3)  operations in node elimination  =',1PD10.3/
     &        ' icntl(8)  scaling effectively used        =',I16/
     &        ' infog(9)  real space for factors          =',I16/
     &        ' infog(10)  Integer space for factors       =',I16/
     &        ' infog(11)  maximum front size              =',I16/
     &        ' infog(29)  number of entries in factors    =',I16)
99985 FORMAT( ' infog(12)  number of off diagonal pivots   =',I16)
99986 FORMAT( ' infog(25)  number of tiny pivots(static)   =',I16)
99988 FORMAT( ' number of 2x2 pivots in type 1 nodes       =',I16)
99989 FORMAT( ' number of 2x2 pivots in type 2 nodes       =',I16)

zmumps_print_allocated_mem()

subroutine zmumps_print_allocated_mem	(	logical, intent(in)	prok,
		logical, intent(in)	prokg,
		logical, intent(in)	print_maxavg,
		integer, intent(in)	mp,
		integer, intent(in)	mpg,
		integer, intent(in)	info16,
		integer, intent(in)	infog18,
		integer, intent(in)	infog19,
		integer, intent(in)	nslaves,
		integer, intent(in)	irank,
		integer, dimension(500), intent(in)	keep )

Definition at line 3849 of file zfac_driver.F.

      IMPLICIT NONE
C
C  Purpose:
C  =======
C     Print memory allocated during factorization
C     - called at beginning of factorization in full-rank
C     - called at end of factorization in low-rank (because
C       of dynamic allocations)
C
      LOGICAL, INTENT(IN) :: PROK, PROKG, PRINT_MAXAVG
      INTEGER, INTENT(IN) :: MP, MPG, INFO16, INFOG18, INFOG19
      INTEGER, INTENT(IN) :: IRANK, NSLAVES
      INTEGER, INTENT(IN) :: KEEP(500)
C
      IF ( prokg ) THEN
        IF (print_maxavg) THEN
         WRITE( mpg,'(A,I12) ')
     &   ' ** Memory allocated, max in Mbytes             (INFOG(18)):',
     &   infog18
        ENDIF
        WRITE( mpg,'(/A,I12) ')
     &   ' ** Memory allocated, total in Mbytes           (INFOG(19)):',
     &    infog19
      END IF
      RETURN