ana_omp_m.F

Go to the documentation of this file.

C
C  This file is part of MUMPS 5.5.1, released
C  on Tue Jul 12 13:17:24 UTC 2022
C
C
C  Copyright 1991-2022 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
C  Mumps Technologies, University of Bordeaux.
C
C  This version of MUMPS is provided to you free of charge. It is
C  released under the CeCILL-C license 
C  (see doc/CeCILL-C_V1-en.txt, doc/CeCILL-C_V1-fr.txt, and
C  https://cecill.info/licences/Licence_CeCILL-C_V1-en.html)
C
      MODULE mumps_ana_omp_m
      CONTAINS
      SUBROUTINE mumps_ana_l0_omp( NB_THREADS, N, NSTEPS, SYM, SLAVEF,
     &  DAD, FRERE, FILS, NSTK_STEPS, ND, STEP, PROCNODE_STEPS, KEEP,
     &  KEEP8, MYID_NODES, NA, LNA, ARITH, LPOOL_B_L0_OMP,
     &  IPOOL_B_L0_OMP, LPOOL_A_L0_OMP, IPOOL_A_L0_OMP,
     &  L_VIRT_L0_OMP, VIRT_L0_OMP, VIRT_L0_OMP_MAPPING,
     &  L_PHYS_L0_OMP, PHYS_L0_OMP, PERM_L0_OMP, PTR_LEAFS_L0_OMP,
     &  THREAD_LA, INFO, ICNTL )
      USE mumps_idll
      USE mumps_ddll
      IMPLICIT NONE
      include 'mpif.h'
      INTEGER, INTENT ( IN ) :: NB_THREADS, N, NSTEPS, SYM
      INTEGER, INTENT ( IN ) :: SLAVEF, MYID_NODES
      INTEGER, INTENT ( IN ) :: LNA
      INTEGER, INTENT ( IN ) :: DAD (:), FRERE (:) 
      INTEGER, INTENT ( IN ) :: FILS (:)  
      INTEGER, INTENT ( IN ) :: NSTK_STEPS (:) 
      INTEGER, INTENT ( IN ) :: ND (:), STEP (:) 
      INTEGER, INTENT ( IN ) :: PROCNODE_STEPS(:) 
      INTEGER, INTENT ( IN ) :: KEEP ( : ) 
      INTEGER(8), INTENT ( IN ) :: KEEP8(:)  
      INTEGER, INTENT ( IN ) :: NA ( : ) 
      CHARACTER(1), INTENT(IN) :: ARITH 
      INTEGER, INTENT ( OUT ) :: LPOOL_B_L0_OMP
      INTEGER, INTENT ( OUT ) :: LPOOL_A_L0_OMP
      INTEGER, INTENT ( OUT ) :: L_PHYS_L0_OMP
      INTEGER, INTENT ( OUT ) :: L_VIRT_L0_OMP
      INTEGER, DIMENSION(:), POINTER :: IPOOL_B_L0_OMP
      INTEGER, DIMENSION(:), POINTER :: IPOOL_A_L0_OMP
      INTEGER, DIMENSION(:), POINTER :: PHYS_L0_OMP
      INTEGER, DIMENSION(:), POINTER :: VIRT_L0_OMP, VIRT_L0_OMP_MAPPING
      INTEGER, DIMENSION(:), POINTER :: PERM_L0_OMP
      INTEGER, DIMENSION(:), POINTER :: PTR_LEAFS_L0_OMP
      INTEGER(8), INTENT ( OUT ) :: THREAD_LA
      INTEGER, INTENT(INOUT) :: INFO(80)
      INTEGER, INTENT(IN)    :: ICNTL(60)
      LOGICAL :: LPOK
      INTEGER :: LP
      INTEGER :: NB_REPEAT_ACCEPTL0, NB_MAX_IN_L0_ACCEPTL0
      INTEGER :: THRESH_MEM, SLAVEF_DURING_MAPPING
      REAL    :: THRESH_EQUILIB
      DOUBLE PRECISION, DIMENSION(1,1,1) :: BENCH 
      INTEGER :: INODE
      INTEGER :: NBLEAF_MYID
      DOUBLE PRECISION :: COST_UNDER, COST_ABOVE, COST_TOTAL_BEST
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: THREADS_CHARGE
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: COSTS_MONO_THREAD
      DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: COSTS_MULTI_THREAD
      INTEGER(8), DIMENSION(:), ALLOCATABLE :: SCHUR_MEMORY
      INTEGER(8), DIMENSION(:), ALLOCATABLE :: SUBTREE_FACTOR_MEMORY
      INTEGER(8), DIMENSION(:), ALLOCATABLE :: SUBTREE_MEMORY
      INTEGER(8) :: FACTOR_SIZE_UNDER_L0, FACTOR_SIZE_PER_MPI
      INTEGER, DIMENSION(:), ALLOCATABLE :: CP_NSTK_STEPS
      TYPE ( IDLL_T ), POINTER :: L0_OMP_DLL
      TYPE ( IDLL_T ), POINTER :: LEAFS_ABOVE_L0_OMP_DLL
      INTEGER :: I
      thresh_equilib = real(keep(408))/real(100)
      IF ((thresh_equilib.GT..99).OR.(thresh_equilib.LT.0.01)) THEN
         thresh_equilib = 0.9
      ENDIF
      thresh_mem = keep(397)
      IF ((thresh_mem.LT.-1).OR.(thresh_mem.GT.100)) thresh_mem=100
      IF (thresh_mem.EQ.-1) THEN
        IF (nb_threads.EQ.2) THEN
         thresh_mem = 50
        ELSEIF (nb_threads.LE.4) THEN
         thresh_mem = 60
        ELSEIF (nb_threads.LT.8) THEN
         thresh_mem = 70
        ELSEIF (nb_threads.LE.12) THEN
         thresh_mem = 80
        ELSEIF (nb_threads.LE.20) THEN
         thresh_mem = 85
        ELSEIF (nb_threads.LE.36) THEN
         thresh_mem = 90
        ELSE
         thresh_mem = 95
        ENDIF
      ENDIF
        slavef_during_mapping = slavef
      factor_size_per_mpi = keep8(101) / slavef_during_mapping
      IF ( keep(261) .EQ. 0) THEN
          WRITE(*,*)"KEEP(261) MUST BE SET TO 1 IN ORDER TO USE
     & MULTITHREADED TREE PARALLELISM"
          CALL mumps_abort()
      END IF
      lp = icntl(1)
      lpok = ( lp .GT. 0 .AND. icntl(4) .GE.1 )
      nb_repeat_acceptl0 = -1
      nb_max_in_l0_acceptl0 = -1
      CALL mumps_ana_initialize_l0_omp ()
      IF (info(1) .LT. 0) GOTO 500
10    CONTINUE
      DO WHILE ( .NOT. mumps_ana_accept_l0_omp() )
        IF (info(1) .LT. 0) GOTO 500
        CALL l0_remove_node ( inode )
          IF (inode .LT. 0) THEN 
              DO i = 1, lpool_b_l0_omp
                  IF (ipool_b_l0_omp( i ) .EQ. (- inode)) THEN
                      ipool_b_l0_omp( i ) = inode 
                      GOTO 10
                  END IF
              END DO
          END IF
          CALL l0_insert_children ( inode )
      END DO
      CALL mumps_ana_finalize_l0_omp ()
 500  CONTINUE
      CALL mumps_ana_free_l0_workspace()
      RETURN
      CONTAINS
      SUBROUTINE mumps_ana_initialize_l0_omp ( )
      IMPLICIT NONE
      INTEGER :: INODE, IFATH, IGRANDFATH, SPECIAL_ROOT,
     &           nfront, npiv, leaf, varnum, ierr
      LOGICAL :: INODE_IS_A_LEAF
      INTEGER(8) :: NFRONT8, NPIV8
      INTEGER(8) :: SUM_CB, MAX_MEM
      DOUBLE PRECISION :: COST_NODE
      LOGICAL :: IN_L0INIT, SKIP_ABOVE
      LOGICAL, EXTERNAL :: MUMPS_ROOTSSARBR, MUMPS_IN_OR_ROOT_SSARBR
      INTEGER, EXTERNAL :: MUMPS_GET_POOL_LENGTH, MUMPS_TYPENODE
      IF (associated(ipool_b_l0_omp)) THEN
        WRITE(*,*) " Internal error 1 MUMPS_ANA_INITIALIZE_L0_OMP" 
        CALL mumps_abort()
      ENDIF
      IF (associated(ipool_a_l0_omp)) THEN
        WRITE(*,*) " Internal error 2 MUMPS_ANA_INITIALIZE_L0_OMP" 
        CALL mumps_abort()
      ENDIF
      IF (associated(virt_l0_omp)) THEN
        WRITE(*,*) " Internal error 3 MUMPS_ANA_INITIALIZE_L0_OMP" 
        CALL mumps_abort()
      ENDIF
      IF (associated(virt_l0_omp_mapping)) THEN
        WRITE(*,*) " Internal error 4 MUMPS_ANA_INITIALIZE_L0_OMP" 
        CALL mumps_abort()
      ENDIF
      IF (associated(perm_l0_omp)) THEN
        WRITE(*,*) " Internal error 5 MUMPS_ANA_INITIALIZE_L0_OMP" 
        CALL mumps_abort()
      ENDIF
      IF (associated(ptr_leafs_l0_omp)) THEN
        WRITE(*,*) " Internal error 6 MUMPS_ANA_INITIALIZE_L0_OMP" 
        CALL mumps_abort()
      ENDIF
      ierr = idll_create( l0_omp_dll )
      ierr = idll_create( leafs_above_l0_omp_dll )
      ALLOCATE( threads_charge( nb_threads ), stat=ierr )
      IF (ierr .GT. 0) THEN
         info(1) = -7
         info(2) = nb_threads
         IF (lpok) WRITE(lp,150) 'THREADS_CHARGE'
         GOTO 500
      ENDIF
      ALLOCATE( costs_mono_thread( nsteps ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1) = -7
         info(2) = nsteps
         IF (lpok) WRITE(lp, 150) ' COSTS_MONO_THREAD'
         GOTO 500
      ENDIF
      ALLOCATE( costs_multi_thread( nsteps ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1) = -7
         info(2) = nsteps
         IF (lpok) WRITE(lp, 150) ' COSTS_MULTI_THREAD'
         GOTO 500
      ENDIF
      ALLOCATE( schur_memory( nsteps ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1) = -7
         info(2) = nsteps
         IF (lpok) WRITE(lp, 150) ' SCHUR_MEMORY'
         GOTO 500
      ENDIF
      ALLOCATE( subtree_factor_memory( nsteps ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1) = -7
         info(2) = nsteps
         IF (lpok) WRITE(lp, 150) ' SCHUR_FACTOR_MEMORY'
         GOTO 500
      ENDIF
      ALLOCATE( subtree_memory( nsteps ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1) = -7
         info(2) = nsteps
         IF (lpok) WRITE(lp, 150) ' SUBTREE_MEMORY'
         GOTO 500
      ENDIF
      ALLOCATE( cp_nstk_steps( nsteps ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1) = -7
         info(2) = nsteps
         IF (lpok) WRITE(lp, 150) ' CP_NSTK_STEPS'
         GOTO 500
      ENDIF
      lpool_b_l0_omp=mumps_get_pool_length(na(1),keep(1),keep8(1))
      ALLOCATE( ipool_b_l0_omp( lpool_b_l0_omp) , stat=ierr )
      IF(ierr.GT.0) THEN
         info(1) = -7
         info(2) = nsteps
         IF (lpok) WRITE(lp, 150) ' id%IPOOL_B_L0_OMP'
         GOTO 500
      ENDIF
      costs_mono_thread = 0.0d0
      costs_multi_thread = 0.0d0
      cost_under = 0.0d0
      cost_above = 0.0d0
      cost_total_best = huge(cost_total_best)
      schur_memory = 0_8
      subtree_factor_memory = 0_8
      subtree_memory = 0_8
      factor_size_under_l0 = 0_8
      cp_nstk_steps(:) = nstk_steps(:)
      IF (keep(403).NE.0) THEN
      CALL read_bench( arith, keep(50) )
      ENDIF
      CALL mumps_init_pool_dist(n, leaf,
     &                     myid_nodes,
     &                     keep(199), na(1), lna, 
     &                     keep(1), keep8(1), step(1),
     &                     procnode_steps(1),
     &                     ipool_b_l0_omp(1), lpool_b_l0_omp)
      leaf = leaf - 1
      nbleaf_myid = leaf
      IF (nbleaf_myid .EQ. 0) THEN
        RETURN
      ENDIF
 90   CONTINUE
      inode = ipool_b_l0_omp( leaf )
      leaf = leaf - 1
      inode_is_a_leaf=.true.
 95   CONTINUE
      nfront = nd( step( inode ) )
      nfront8= int(nfront,8)
      npiv = 0
      varnum = inode
      DO WHILE (varnum .GT. 0 )
          npiv = npiv + 1
          varnum = fils( varnum )
      END DO
      npiv8=int(npiv,8)
      varnum = - varnum
      IF (keep(50) .EQ. 0) THEN
        schur_memory( step( inode ) ) =
     &           (nfront8 - npiv8)*(nfront8 - npiv8)
        IF (keep(251) .EQ. 0) THEN
          subtree_factor_memory( step( inode ) ) = nfront8 * nfront8
     &           - schur_memory( step( inode ) )
        ELSE
          subtree_factor_memory( step( inode ) ) = 0_8
        END IF
      ELSE
        schur_memory( step( inode ) ) =
     &           (nfront8 - npiv8)*(nfront8 + 1_8 - npiv8)/2_8
        IF (keep(251) .EQ. 0) THEN
          subtree_factor_memory( step( inode ) ) = nfront8 * npiv8
        ELSE
          subtree_factor_memory( step( inode ) ) = 0_8
        END IF
      END IF
      sum_cb  = 0_8
      max_mem = 0_8
      IF (keep(403) .EQ. 0) THEN
        CALL mumps_get_flops_cost ( nfront, npiv, npiv,
     &                              sym, 1, cost_node )
        costs_mono_thread( step( inode ) ) = cost_node
      ELSE
        CALL cost_bench (npiv, nfront-npiv, 1, keep(50), cost_node)
        costs_mono_thread( step( inode ) ) = cost_node
        CALL cost_bench (npiv,nfront-npiv,nb_threads,keep(50),cost_node)
        costs_multi_thread( step( inode ) ) = cost_node
      END IF
      DO WHILE (varnum .GT. 0 )
          costs_mono_thread( step( inode ) ) =
     &                      costs_mono_thread( step( inode ) )
     &                      +
     &                      costs_mono_thread( step( varnum ) )
          max_mem = max(max_mem,
     &                   subtree_memory( step( varnum ) ) + sum_cb )
          sum_cb = sum_cb + schur_memory( step( varnum ) ) +
     &             subtree_factor_memory( step( varnum ) )
          subtree_factor_memory( step( inode ) ) =
     &                         subtree_factor_memory( step( inode ) )
     &                       + subtree_factor_memory( step( varnum ) )
          varnum = frere( step( varnum ) )
      END DO
      subtree_memory( step( inode ) ) =
     &                        max( max_mem, nfront8*nfront8 + sum_cb )
      ifath = dad( step( inode ) )
      IF (ifath .NE. 0) THEN
        igrandfath = dad( step( ifath ) )
      ELSE
        igrandfath = 0
      ENDIF
      special_root = max(keep(38), keep(20))
      skip_above = .false.
      in_l0init  = .false.
      IF ( inode .EQ. special_root ) THEN
        in_l0init  = .false.
        IF (inode_is_a_leaf) THEN
          skip_above = .true.
          GOTO 80
        ELSE
          WRITE(*,*) " Internal error 1 in MUMPS_ANA_INITIALIZE_L0_OMP",
     &    inode, special_root
          CALL mumps_abort()
        ENDIF
      ENDIF
      IF ( ifath .NE. 0 .AND. ifath .EQ. keep(38) ) THEN
        in_l0init  = .false.
        IF (inode_is_a_leaf) THEN
          skip_above = .true.
          GOTO 80
        ELSE
          WRITE(*,*) " Internal error 2 in MUMPS_ANA_INITIALIZE_L0_OMP",
     &    inode, ifath, keep(38)
          CALL mumps_abort()
        ENDIF
      ENDIF
      IF ( slavef_during_mapping > 1 ) THEN
        IF (mumps_rootssarbr(
     &  procnode_steps( step( inode ) ), keep(199) )
     &  .OR. .NOT. mumps_in_or_root_ssarbr(
     &  procnode_steps( step ( inode ) ), keep(199) )
     &) THEN
          in_l0init = .false.
          IF (inode_is_a_leaf) THEN
            skip_above = .true.
            GOTO 80
          ELSE
            WRITE(*,*)
     &      " Internal error 3 in MUMPS_ANA_INITIALIZE_L0_OMP",
     &      inode
            CALL mumps_abort()
          ENDIF
        ENDIF
      ENDIF
      IF (ifath.NE.0) THEN
        IF ( mumps_typenode(step(ifath),keep(199)).EQ.2) THEN
          in_l0init = .false.
          IF (inode_is_a_leaf) THEN
            skip_above = .true.
            GOTO 80
          ELSE
            WRITE(*,*)
     &      " Internal error 5 in MUMPS_ANA_INITIALIZE_L0_OMP",
     &      inode, ifath
            CALL mumps_abort()
          ENDIF
        ENDIF
      ENDIF
      IF ( mumps_typenode(step(inode),keep(199)).EQ.2) THEN
        in_l0init = .false.
        IF (inode_is_a_leaf) THEN
          skip_above = .true.
          GOTO 80
        ELSE
          WRITE(*,*)
     &    " Internal error 6 in MUMPS_ANA_INITIALIZE_L0_OMP",
     &    inode
          CALL mumps_abort()
        ENDIF
      ENDIF
      IF ( ifath .EQ. 0 ) THEN
        in_l0init = .true.
        GOTO 80
      ELSE
        IF ( ifath .EQ. keep(20) ) THEN
          in_l0init = .true.
          GOTO 80
        ENDIF
        IF ( igrandfath .EQ. keep(38) .AND. keep(38) .NE. 0 ) THEN
          in_l0init = .true.
          GOTO 80
        ENDIF
        IF ( slavef_during_mapping > 1 ) THEN
          IF (mumps_rootssarbr(
     &    procnode_steps( step( ifath ) ), keep(199) )) THEN
            in_l0init = .true.
            GOTO 80
          ENDIF
        ENDIF
      ENDIF
 80   CONTINUE
      IF (.NOT. skip_above) THEN
        IF (keep(50).EQ.0) THEN
          factor_size_under_l0 = factor_size_under_l0 +
     &    npiv8 * ( nfront8 + nfront8 - npiv8 )
        ELSE
          factor_size_under_l0  = factor_size_under_l0  +
     &    nfront8 * npiv8
        ENDIF
      ENDIF
      IF ( in_l0init ) THEN
          CALL l0_insert_node ( l0_omp_dll, inode )
      ELSE IF ( skip_above ) THEN
          ierr = idll_push_back( leafs_above_l0_omp_dll, inode )
          IF ( .NOT. inode_is_a_leaf ) THEN
            WRITE(*,*)
     &      " Internal error 7 in MUMPS_ANA_INITIALIZE_L0_OMP",
     &      inode
            CALL mumps_abort()
          ENDIF
          ipool_b_l0_omp(leaf+1) = -inode
      ELSE
          cp_nstk_steps( step( ifath ) ) =
     &                          cp_nstk_steps( step( ifath ) ) - 1
          IF ( cp_nstk_steps( step( ifath ) ) .EQ. 0 ) THEN
              inode = ifath
              inode_is_a_leaf = .false.
              GOTO 95
          ENDIF
      END IF
      IF ( leaf .GT. 0 ) THEN
        GOTO 90
      END IF
 500  CONTINUE
      RETURN
 150  FORMAT(
     & /' ** ALLOC FAILURE IN MUMPS_ANA_INITIALIZE_L0_OMP FOR ',
     & a30)
      END SUBROUTINE mumps_ana_initialize_l0_omp
      SUBROUTINE l0_insert_node ( DLL, INODE )
      IMPLICIT NONE
      INTEGER, INTENT ( IN ) :: INODE
      TYPE ( IDLL_T ), POINTER :: DLL
      INTEGER :: IERR
      TYPE ( IDLL_NODE_T ), POINTER :: IDLL_NODE
      ierr = idll_iterator_begin( dll, idll_node )
      DO WHILE ( associated ( idll_node ) )
          IF ( costs_mono_thread( step( idll_node%ELMT ) )
     &         .GT.
     &         costs_mono_thread( step( inode ) ) ) THEN
              idll_node => idll_node%NEXT
          ELSE
              EXIT
          END IF
      END DO
      IF ( .NOT. associated ( idll_node ) ) THEN
          ierr = idll_push_back(dll, inode)
      ELSE
          ierr = idll_insert_before(dll, idll_node, inode)
      ENDIF
      RETURN
      END SUBROUTINE l0_insert_node
      SUBROUTINE l0_insert_children ( I_FATHER )
      IMPLICIT NONE
      INTEGER, INTENT ( IN ) :: I_FATHER
      INTEGER :: I_SON, IERR
      TYPE ( IDLL_T ), POINTER :: SON_DLL
      TYPE ( IDLL_NODE_T ), POINTER :: IDLL_NODE
      ierr = idll_create( son_dll )
      i_son = i_father
      DO WHILE ( i_son .GT. 0 )
          i_son = fils( i_son )
      END DO
      i_son = - i_son
      IF ( i_son .EQ. 0 ) THEN
          RETURN
      END IF
      DO WHILE ( i_son .GT. 0 )
          CALL l0_insert_node ( son_dll, i_son )
           i_son = frere( step( i_son ) )
      END DO
      ierr = idll_iterator_begin( l0_omp_dll, idll_node )
      ierr = idll_pop_front( son_dll, i_son )
      IF ( ierr .NE. 0 ) THEN
          GOTO 190
      END IF
      IF ( .NOT. associated( idll_node ) ) THEN
              DO
                  ierr = idll_push_back( l0_omp_dll, i_son )
                  ierr = idll_pop_front( son_dll, i_son )
                  IF ( ierr .NE. 0 ) THEN
                      GOTO 190
                  END IF
              END DO
      ELSE
          DO
              IF ( costs_mono_thread( step( i_son )) .LE.
     &             costs_mono_thread( step( idll_node%ELMT ) ) ) THEN
                  IF ( associated ( idll_node%NEXT ) ) THEN
                      idll_node => idll_node%NEXT
                  ELSE
                      ierr = idll_push_back(l0_omp_dll, i_son)
                      ierr = idll_pop_front( son_dll, i_son )
                      IF ( ierr .NE. 0 ) THEN
                          GOTO 190
                      END IF
                  END IF
              ELSE
                  ierr = idll_insert_before(l0_omp_dll, idll_node,i_son)
                  ierr = idll_pop_front( son_dll, i_son )
                  IF ( ierr .NE. 0 ) THEN
                      GOTO 190
                  END IF
              END IF
          END DO
      END IF
190   CONTINUE
      ierr = idll_destroy( son_dll )
      RETURN
      END SUBROUTINE l0_insert_children
      SUBROUTINE l0_remove_node ( INODE )
      IMPLICIT NONE
      INTEGER, INTENT ( OUT ) :: INODE
      INTEGER :: I_SON, IERR, NPIV
      ierr = idll_pop_front( l0_omp_dll, inode )
      i_son = inode
      npiv = 0
      DO WHILE ( i_son .GT. 0 )
          npiv = npiv + 1
          i_son = fils( i_son )
      END DO
      i_son = - i_son
      IF (keep(50) .EQ. 0) THEN
        factor_size_under_l0 = factor_size_under_l0 -
     &  int(npiv, 8) * int(2 *  nd(step(inode)) - npiv, 8)
      ELSE
        factor_size_under_l0 = factor_size_under_l0 -
     &  int(npiv, 8) * int(nd(step(inode)), 8)
      ENDIF
      IF ( i_son .EQ. 0 ) THEN 
          ierr = idll_push_back( leafs_above_l0_omp_dll, inode )
          inode = -inode
      ELSE IF (inode .GT. 0) THEN
          cost_above = cost_above + costs_multi_thread(step( inode ))
      END IF
      RETURN
      END SUBROUTINE l0_remove_node
      FUNCTION mumps_ana_accept_l0_omp ()
      LOGICAL :: mumps_ana_accept_l0_omp
      INTEGER :: i, i_less_charged, IERR, nb_in_l0
      DOUBLE PRECISION :: lightest_charge, heaviest_charge
      TYPE ( idll_node_t ), POINTER :: idll_node
      threads_charge = 0.0d0
      nb_in_l0 = 0 
      ierr = idll_iterator_begin( l0_omp_dll, idll_node )
      DO WHILE ( associated ( idll_node ) )
          nb_in_l0 = nb_in_l0 + 1
          i_less_charged = 1
          lightest_charge = threads_charge( 1 )
          DO i = 2, nb_threads
              IF ( threads_charge( i ) .LT. lightest_charge ) THEN
                  i_less_charged = i
                  lightest_charge = threads_charge( i )
              END IF
          END DO
          threads_charge( i_less_charged ) =
     &                   threads_charge( i_less_charged )
     &                   +
     &                   costs_mono_thread( step( idll_node%ELMT ) )
          idll_node => idll_node%NEXT
      END DO
      nb_max_in_l0_acceptl0 = max(nb_max_in_l0_acceptl0, nb_in_l0)
      lightest_charge = threads_charge( 1 )
      heaviest_charge = threads_charge( 1 )
      DO i = 2, nb_threads
          IF ( threads_charge( i ) .LT. lightest_charge ) THEN
              lightest_charge = threads_charge( i )
          ELSEIF ( threads_charge( i ) .GT. heaviest_charge ) THEN
              heaviest_charge = threads_charge( i )
          END IF
      END DO
      cost_under = heaviest_charge
      IF (keep(403) .EQ. 0) THEN
          mumps_ana_accept_l0_omp =
     &   (
     &     dble(lightest_charge)/(dble(heaviest_charge)+1.d-12) 
     &     .GT.thresh_equilib .AND.
     &
     &     factor_size_under_l0 .LE.
     &     factor_size_per_mpi * int(thresh_mem,8) / 100_8
     &
     &      )
     &  .OR.
     &   ( nb_in_l0 .LT. nb_max_in_l0_acceptl0 .AND.
     &     lightest_charge .EQ. 0.0d0 ) 
     &  .OR. ( nb_in_l0 .EQ. 0 ) 
        IF (mumps_ana_accept_l0_omp) THEN
          IF (associated(phys_l0_omp)) THEN
            DEALLOCATE(phys_l0_omp)
            nullify(phys_l0_omp)
          ENDIF
          ierr = idll_2_array( l0_omp_dll, phys_l0_omp, l_phys_l0_omp )
          IF (ierr .EQ. -2) THEN
            info(1) = -7
            info(2) = l_phys_l0_omp
            RETURN
          ENDIF
        END IF
      ELSE
        IF (cost_under + cost_above .LT. cost_total_best) THEN
          IF (associated(phys_l0_omp)) THEN
            DEALLOCATE(phys_l0_omp)
            nullify(phys_l0_omp)
          ENDIF
          ierr = idll_2_array( l0_omp_dll, phys_l0_omp, l_phys_l0_omp )
          cost_total_best = cost_under + cost_above
          nb_repeat_acceptl0 = 100
        END IF
        nb_repeat_acceptl0 = nb_repeat_acceptl0- 1
        mumps_ana_accept_l0_omp = (nb_repeat_acceptl0 .EQ. 0)
      END IF
      RETURN
      END FUNCTION mumps_ana_accept_l0_omp
      SUBROUTINE mumps_ana_finalize_l0_omp ()
      IMPLICIT NONE
      INTEGER :: INODE, OLD_INODE, I, J, K, LEAF, IERR
      DOUBLE PRECISION :: LIGHTEST_CHARGE
      INTEGER :: I_LESS_CHARGED
      INTEGER(8) :: SUM_CB, MAX_MEM, MAX_MEM_ALL_THREADS
      INTEGER :: MAX_TASK_PER_THREAD
      TYPE ( IDLL_NODE_T ), POINTER :: IDLL_NODE
      INTEGER, DIMENSION(:,:), ALLOCATABLE :: THREADS_TASKS
      INTEGER, DIMENSION(:), ALLOCATABLE :: NB_TASK_PER_THREAD
      INTEGER, DIMENSION(:), ALLOCATABLE :: INV_PERM_L0_OMP
      EXTERNAL :: mumps_get_pool_length
      INTEGER :: MUMPS_GET_POOL_LENGTH
      IF (keep(402) .EQ. 0) THEN
        l_virt_l0_omp = nb_threads + 1
      ELSE
        l_virt_l0_omp = l_phys_l0_omp + 1
      END IF
      lpool_a_l0_omp = mumps_get_pool_length(na(1),keep(1),keep8(1))
      ALLOCATE ( virt_l0_omp( max(l_virt_l0_omp,1) ),
     &           virt_l0_omp_mapping( max(l_virt_l0_omp,1) ),
     &           stat=ierr )
      IF(ierr.GT.0) THEN
         info(1)=-7
         info(2)=2*max(l_virt_l0_omp,1)
         IF (lpok) WRITE(lp,150) 'id%VIRT_L0_OMP[_MAPPING]'
         GOTO 300
      ENDIF
      ALLOCATE ( perm_l0_omp( max(l_phys_l0_omp,1) ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1)=-7
         info(2)=max(l_phys_l0_omp,1)
         IF (lpok) WRITE(lp,150) 'id%PERM_L0_OMP'
         GOTO 300
      ENDIF
      ALLOCATE ( ptr_leafs_l0_omp( l_phys_l0_omp + 1 ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1)=-7
         info(2)=max(l_phys_l0_omp,1)
         IF (lpok) WRITE(lp,150) 'id%PTR_LEAFS_L0_OMP'
         GOTO 300
      ENDIF
      ALLOCATE ( ipool_a_l0_omp( lpool_a_l0_omp ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1)=-7
         info(2)=lpool_a_l0_omp
         IF (lpok) WRITE(lp,150) 'id%IPOOL_A_L0_OMP'
         GOTO 300
      ENDIF
      ALLOCATE ( nb_task_per_thread( nb_threads ), stat=ierr )
      IF(ierr.GT.0) THEN
         info(1)=-7
         info(2)=nb_threads
         IF (lpok) WRITE(lp,150) 'NB_TASK_PER_THREAD'
         GOTO 300
      ENDIF
      ALLOCATE ( inv_perm_l0_omp( l_phys_l0_omp ), stat=ierr )
      IF(ierr.GT.0) THEN
         WRITE(*,*) "Allocation Error in MUMPS_ANA_FINALIZE_L0_OMP"
         CALL mumps_abort()
      ENDIF
      nb_task_per_thread = 0
      threads_charge = 0.0d0
      DO i = 1, l_phys_l0_omp
          i_less_charged = 1
          lightest_charge = threads_charge( 1 )
          DO j = 2, nb_threads
              IF ( threads_charge( j ) .LT. lightest_charge ) THEN
                  i_less_charged = j
                  lightest_charge = threads_charge( j )
                  IF (threads_charge( j ) .EQ. 0) THEN
                    EXIT
                  ENDIF
              END IF
          END DO
          nb_task_per_thread( i_less_charged ) =
     &                        nb_task_per_thread( i_less_charged ) + 1
          IF (keep(402) .NE. 0) THEN
            virt_l0_omp_mapping(i) = i_less_charged
          ENDIF
          threads_charge( i_less_charged ) =
     &                  threads_charge( i_less_charged )
     &                  +
     &                  costs_mono_thread( step( phys_l0_omp( i ) ) )
      END DO
      IF (keep(402) .EQ. 0) THEN
        DO i = 1, nb_threads
          virt_l0_omp_mapping(i) = i
        ENDDO
      ENDIF
      virt_l0_omp_mapping(l_virt_l0_omp) = -999999
      max_task_per_thread = 0
      DO i = 1, nb_threads 
          max_task_per_thread = max(max_task_per_thread,
     &                               nb_task_per_thread( i ) )
      END DO
      ALLOCATE ( threads_tasks( nb_threads, max_task_per_thread ),
     &     stat=ierr )
      IF(ierr.GT.0) THEN
         info(1)=-7
         info(2)=nb_threads*max_task_per_thread
         IF (lpok) WRITE(lp,150) 'THREADS_TASK'
         GOTO 300
      ENDIF
      nb_task_per_thread = 0
      threads_charge = 0.0d0
      threads_tasks = 0 
      DO i = 1, l_phys_l0_omp
          i_less_charged = 1
          lightest_charge = threads_charge( 1 )
          DO j = 2, nb_threads
              IF ( threads_charge( j ) .LT. lightest_charge ) THEN
                  i_less_charged = j
                  lightest_charge = threads_charge( j )
              END IF
          END DO
          nb_task_per_thread( i_less_charged ) =
     &                         nb_task_per_thread( i_less_charged ) + 1
          threads_tasks( i_less_charged, nb_task_per_thread
     &                    ( i_less_charged ) ) = phys_l0_omp( i )
          threads_charge( i_less_charged ) =
     &                  threads_charge( i_less_charged )
     &                  +
     &                  costs_mono_thread( step( phys_l0_omp( i ) ) )
      END DO
      max_mem_all_threads = 0_8
      DO i = 1, nb_threads
          sum_cb = 0_8
          max_mem = 0_8
          DO j = 1, nb_task_per_thread( i )
                  max_mem = max( max_mem, subtree_memory( step(
     &                threads_tasks(i,j) ) ) + sum_cb )
                  sum_cb = sum_cb
     &              +schur_memory(step(threads_tasks(i,j)))
     &              +subtree_factor_memory(
     &                            step(threads_tasks(i,j)))
          END DO
          max_mem = max( max_mem, sum_cb )
          IF (keep(402) .EQ. 0) THEN
              threads_charge( i ) = dble(max_mem)
          END IF
          max_mem_all_threads =  max( max_mem_all_threads, max_mem )
      END DO
      max_mem_all_threads =  ( max_mem_all_threads 
     &                         * int(100 + keep(12),8) ) / 100_8
      thread_la = max(max_mem_all_threads,6_8)
      IF (keep(402) .EQ. 0) THEN
        k = 1
        DO i = 1, nb_threads
          virt_l0_omp(i) = k
          DO j = 1, nb_task_per_thread( i )
              phys_l0_omp(k) = threads_tasks(i,j)
              k = k + 1
          END DO
        END DO
        virt_l0_omp(nb_threads+1) = k
      ELSE
        DO i = 1, l_virt_l0_omp 
            virt_l0_omp(i) = i
        END DO
      END IF
      DO i = 1, l_phys_l0_omp
          inv_perm_l0_omp( i ) = i
      END DO
      IF ( l_phys_l0_omp .GT. 1 ) THEN
        CALL mumps_quick_sort_phys_l0( n, step(1), phys_l0_omp(1),
     &     inv_perm_l0_omp, l_phys_l0_omp, 1, l_phys_l0_omp )
      ENDIF
      DO i = 1, l_phys_l0_omp
          perm_l0_omp( inv_perm_l0_omp( i ) ) = i
      END DO
      j = nbleaf_myid 
      ptr_leafs_l0_omp( 1 ) = j
      DO i = 1, l_phys_l0_omp 
          old_inode = 0
          inode = phys_l0_omp( i )
          DO WHILE ( inode .NE. 0 )
              old_inode = inode
              DO WHILE ( inode .GT. 0 )
                  inode = fils( inode )
              END DO
              inode = - inode
          END DO
          DO WHILE ( ipool_b_l0_omp( j ) .NE. old_inode )
              j = j - 1
          END DO
          j = j - 1
          ptr_leafs_l0_omp( i + 1 ) = j
      END DO
      cp_nstk_steps(:) = nstk_steps(:)
      ipool_a_l0_omp = 0
      leaf = 1
      ierr = idll_iterator_begin( leafs_above_l0_omp_dll, idll_node )
      DO WHILE ( associated( idll_node ) )
          ipool_a_l0_omp( leaf ) = idll_node%ELMT
          leaf = leaf + 1
          idll_node => idll_node%NEXT
      END DO
      DO i = 1 , l_phys_l0_omp 
        IF ( dad( step( phys_l0_omp(i) ) ) .NE. 0 ) THEN
          cp_nstk_steps( step( dad( step( phys_l0_omp(i) ) ) ) ) =
     &     cp_nstk_steps( step( dad( step( phys_l0_omp(i) ) ) ) )-1
          IF (cp_nstk_steps(step(dad(step(phys_l0_omp(i))))) .EQ. 0)THEN
            ipool_a_l0_omp( leaf ) = dad(step(phys_l0_omp( i )))
            leaf = leaf + 1
          END IF
        END IF
      END DO
      leaf = leaf - 1 
      ipool_a_l0_omp(lpool_a_l0_omp) = leaf
      ipool_a_l0_omp(lpool_a_l0_omp-1) = 0
      ipool_a_l0_omp(lpool_a_l0_omp-2) = 0
      IF (leaf .GT. 1) THEN
        CALL mumps_quick_sort_ipool_po( n, step(1), 
     &            ipool_a_l0_omp(1), leaf, 1, leaf )
      ENDIF
 300  CONTINUE
      IF (allocated(nb_task_per_thread)) DEALLOCATE (nb_task_per_thread)
      IF (allocated(inv_perm_l0_omp   )) DEALLOCATE ( inv_perm_l0_omp )
      IF (allocated(threads_tasks     )) DEALLOCATE (threads_tasks )
      RETURN
 150  FORMAT(
     & /' ** ALLOC FAILURE IN MUMPS_ANA_FINALIZE_L0_OMP FOR ',
     & a30)
      END SUBROUTINE mumps_ana_finalize_l0_omp 
      SUBROUTINE mumps_ana_free_l0_workspace()
      INTEGER :: IERR
      IF (allocated(threads_charge))     DEALLOCATE(threads_charge    )
      IF (allocated(cp_nstk_steps ))     DEALLOCATE(cp_nstk_steps     )
      IF (allocated(costs_mono_thread))  DEALLOCATE(costs_mono_thread )
      IF (allocated(costs_multi_thread)) DEALLOCATE(costs_multi_thread)
      IF (allocated(schur_memory))       DEALLOCATE(schur_memory      )
      IF (allocated(subtree_factor_memory))
     &                                DEALLOCATE(subtree_factor_memory)
      IF (allocated(subtree_memory))     DEALLOCATE(subtree_memory    )
      ierr = idll_destroy( leafs_above_l0_omp_dll )
      ierr = idll_destroy( l0_omp_dll )
      RETURN
      END SUBROUTINE mumps_ana_free_l0_workspace
      SUBROUTINE read_bench(ARITH, K50)
      IMPLICIT NONE
      INTEGER,      INTENT(in) :: K50
      CHARACTER(1), INTENT(in) :: ARITH
      INTEGER NLINES, INDEX_NPIV, INDEX_NSCHUR, NB_CORE
      INTEGER V, S, OLD_V, OLD_S, I
      parameter(nlines=2812)
      DOUBLE PRECISION :: AUX
      CHARACTER(1)             :: K50_STR
      index_npiv = 0
      index_nschur = 0
      old_v = -1
      old_s = -1
      WRITE(k50_str,'(I1)') k50
      OPEN(1,file=arith//'benchmark_sym_'//k50_str//'.csv')
      DO I=1,NLINES
        READ(1,*) V, S, NB_CORE, AUX
.NE.        IF (V  OLD_V) THEN
            INDEX_NPIV = INDEX_NPIV + 1
            OLD_V = V
        END IF
.GT.        IF (S  OLD_S) THEN
            INDEX_NSCHUR = INDEX_NSCHUR + 1
            OLD_S = S
.LT.        ELSEIF (S  OLD_S) THEN
            INDEX_NSCHUR = 1
            OLD_S = S
        END IF
        BENCH (INDEX_NPIV, INDEX_NSCHUR, NB_CORE) = AUX
      END DO
      CLOSE(1)
      RETURN
      END SUBROUTINE READ_BENCH
      SUBROUTINE COST_BENCH (NPIV, NSCHUR, NB_CORE, SYM, COST)
      IMPLICIT NONE
      INTEGER, INTENT(IN) ::  NPIV, NSCHUR, NB_CORE, SYM
      DOUBLE PRECISION, INTENT(OUT) :: COST
      INTEGER V, VV, S, SS
      INTEGER LOW_INDEX_NPIV, LOW_INDEX_NSCHUR
      INTEGER HIGH_INDEX_NPIV, HIGH_INDEX_NSCHUR
      DOUBLE PRECISION :: APROX_COST_FLOPS, REAL_COST_FLOPS
.LE.      IF (NPIV  10) THEN
        LOW_INDEX_NPIV = NPIV
        V = NPIV
        VV = NPIV + 1
.LE.      ELSEIF (NPIV  100) THEN
        LOW_INDEX_NPIV = 9 + NPIV/10
        V = (NPIV/10)*10
        VV = (NPIV/10+1)*10
.LE.      ELSEIF (NPIV  1000) THEN
        LOW_INDEX_NPIV = 18 + NPIV/100
        V = (NPIV/100)*100
        VV = (NPIV/100+1)*100
.LE.      ELSEIF (NPIV  10000) THEN
        LOW_INDEX_NPIV = 27 + NPIV/1000
        V = (NPIV/1000)*1000
        VV = (NPIV/1000+1)*1000
      ELSE
        LOW_INDEX_NPIV = 37
        V = (NPIV/10000)*10000
        VV = (NPIV/10000+1)*10000
      END IF
.LE.      IF (NSCHUR  10) THEN
        LOW_INDEX_NSCHUR = NSCHUR + 1
        S = NSCHUR
        SS = NSCHUR + 1
.LE.      ELSEIF (NSCHUR  100) THEN
        LOW_INDEX_NSCHUR = 10 + NSCHUR/10
        S = (NSCHUR/10)*10
        SS = (NSCHUR/10+1)*10
.LE.      ELSEIF (NSCHUR  1000) THEN
        LOW_INDEX_NSCHUR = 19 + NSCHUR/100
        S = (NSCHUR/100)*100
        SS = (NSCHUR/100+1)*100
.LE.      ELSEIF (NSCHUR  10000) THEN
        LOW_INDEX_NSCHUR = 28 + NSCHUR/1000
        S = (NSCHUR/1000)*1000
        SS = (NSCHUR/1000+1)*1000
      ELSE
        LOW_INDEX_NSCHUR = 38
        S = (NSCHUR/10000)*10000
        SS = (NSCHUR/10000+1)*10000
      END IF
.LT.      IF (V  10000) THEN
.LT.        IF (S  10000) THEN
          HIGH_INDEX_NPIV = LOW_INDEX_NPIV + 1
          HIGH_INDEX_NSCHUR = LOW_INDEX_NSCHUR + 1
          COST = (BENCH(LOW_INDEX_NPIV, LOW_INDEX_NSCHUR, NB_CORE)
     &             *(VV - NPIV)*(SS - NSCHUR)
     &            +BENCH(LOW_INDEX_NPIV, HIGH_INDEX_NSCHUR, NB_CORE)
     &             *(VV - NPIV)*(NSCHUR - S)
     &            +BENCH(HIGH_INDEX_NPIV, LOW_INDEX_NSCHUR, NB_CORE)
     &             *(NPIV - V)*(SS - NSCHUR)
     &            +BENCH(HIGH_INDEX_NPIV, HIGH_INDEX_NSCHUR, NB_CORE)
     &             *(NPIV - V)*(NSCHUR - S))
     &           /((VV - V)*(SS - S))
        ELSE
          HIGH_INDEX_NPIV = LOW_INDEX_NPIV + 1
          HIGH_INDEX_NSCHUR = LOW_INDEX_NSCHUR
          COST = (BENCH(LOW_INDEX_NPIV, LOW_INDEX_NSCHUR, NB_CORE)
     &             *(VV - NPIV)
     &            +BENCH(HIGH_INDEX_NPIV, LOW_INDEX_NSCHUR, NB_CORE)
     &             *(NPIV - V))
     &           /(VV - V)
          CALL MUMPS_GET_FLOPS_COST ( NPIV+NSCHUR, NPIV, NPIV,
     &                                SYM, 1, REAL_COST_FLOPS )
          CALL MUMPS_GET_FLOPS_COST ( V+S, V, V,
     &                                SYM, 1, APROX_COST_FLOPS )
          COST = COST * (REAL_COST_FLOPS/APROX_COST_FLOPS)
        END IF
      ELSE
        IF (NSCHUR < 10000) THEN
          HIGH_INDEX_NPIV = LOW_INDEX_NPIV
          HIGH_INDEX_NSCHUR = LOW_INDEX_NSCHUR + 1
          COST = (BENCH(LOW_INDEX_NPIV, LOW_INDEX_NSCHUR, NB_CORE)
     &             *(SS - NSCHUR)
     &            +BENCH(LOW_INDEX_NPIV, HIGH_INDEX_NSCHUR, NB_CORE)
     &             *(NSCHUR - S))
     &           /(SS - S)
          CALL MUMPS_GET_FLOPS_COST ( NPIV+NSCHUR, NPIV, NPIV,
     &                                SYM, 1, REAL_COST_FLOPS )
          CALL MUMPS_GET_FLOPS_COST ( V+S, V, V,
     &                                SYM, 1, APROX_COST_FLOPS )
          COST = COST * (REAL_COST_FLOPS/APROX_COST_FLOPS)
        ELSE
          HIGH_INDEX_NPIV = LOW_INDEX_NPIV
          HIGH_INDEX_NSCHUR = LOW_INDEX_NSCHUR
          COST = (BENCH(LOW_INDEX_NPIV, LOW_INDEX_NSCHUR, NB_CORE))
          CALL MUMPS_GET_FLOPS_COST ( NPIV+NSCHUR, NPIV, NPIV,
     &                                SYM, 1, REAL_COST_FLOPS )
          CALL MUMPS_GET_FLOPS_COST ( V+S, V, V,
     &                                SYM, 1, APROX_COST_FLOPS )
          COST = COST * (REAL_COST_FLOPS/APROX_COST_FLOPS)
        END IF
      END IF
      END SUBROUTINE COST_BENCH
      END SUBROUTINE MUMPS_ANA_L0_OMP
      END MODULE MUMPS_ANA_OMP_M
      RECURSIVE SUBROUTINE MUMPS_QUICK_SORT_IPOOL_PO( N, STEP, 
     &            INTLIST, TAILLE, LO, HI )
      IMPLICIT NONE
      INTEGER N, TAILLE
      INTEGER STEP( N )
      INTEGER INTLIST( TAILLE )
      INTEGER LO, HI
      INTEGER I,J
      INTEGER ISWAP, PIVOT
      I = LO
      J = HI
      PIVOT = STEP(INTLIST((I+J)/2))
 10   IF (STEP(INTLIST(I)) > PIVOT) THEN
        I=I+1
        GOTO 10
      ENDIF
 20   IF (STEP(INTLIST(J)) < PIVOT) THEN
        J=J-1
        GOTO 20
      ENDIF
      IF (I < J) THEN
        ISWAP = INTLIST(I)
        INTLIST(I) = INTLIST(J)
        INTLIST(J)=ISWAP
      ENDIF
      IF ( I <= J) THEN
        I = I+1
        J = J-1
      ENDIF
      IF ( I <= J ) GOTO 10
      IF ( LO < J ) CALL MUMPS_QUICK_SORT_IPOOL_PO(N, STEP,
     &              INTLIST, TAILLE, LO, J)
      IF ( I < HI ) CALL MUMPS_QUICK_SORT_IPOOL_PO(N, STEP,
     &              INTLIST, TAILLE, I, HI)
      RETURN
      END SUBROUTINE MUMPS_QUICK_SORT_IPOOL_PO
      RECURSIVE SUBROUTINE MUMPS_QUICK_SORT_PHYS_L0( N, STEP, 
     &            INTLIST, INVPERM, TAILLE, LO, HI )
      IMPLICIT NONE
      INTEGER N, TAILLE
      INTEGER STEP( N )
      INTEGER INTLIST( TAILLE )
      INTEGER INVPERM( TAILLE )
      INTEGER LO, HI
      INTEGER I,J
      INTEGER ISWAP, PIVOT
      INTEGER dswap
      I = LO
      J = HI
      PIVOT = STEP(INTLIST((I+J)/2))
 10   IF (STEP(INTLIST(I)) < PIVOT) THEN
        I=I+1
        GOTO 10
      ENDIF
 20   IF (STEP(INTLIST(J)) > PIVOT) THEN
        J=J-1
        GOTO 20
      ENDIF
      IF (I < J) THEN
        ISWAP = INTLIST(I)
        INTLIST(I) = INTLIST(J)
        INTLIST(J)=ISWAP
        dswap = INVPERM(I)
        INVPERM(I) = INVPERM(J)
        INVPERM(J) = dswap
      ENDIF
      IF ( I <= J) THEN
        I = I+1
        J = J-1
      ENDIF
      IF ( I <= J ) GOTO 10
      IF ( LO < J ) CALL MUMPS_QUICK_SORT_PHYS_L0(N, STEP,
     &              INTLIST, INVPERM, TAILLE, LO, J)
      IF ( I < HI ) CALL MUMPS_QUICK_SORT_PHYS_L0(N, STEP,
     &              INTLIST, INVPERM, TAILLE, I, HI)
      RETURN
      END SUBROUTINE MUMPS_QUICK_SORT_PHYS_L0
      SUBROUTINE MUMPS_ANA_OMP_RETURN()
#if defined(BLR_MT)
#if ! defined(_OPENMP)
      COMPILATION FAILURE: -DBLR_MT requires compilation with openmp
      Please modify Makefile.inc and do 'make clean; make'
#endif
#endif
      RETURN
      END SUBROUTINE MUMPS_ANA_OMP_RETURN