cmumps_sol_es.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 cmumps_sol_es
      PRIVATE
      PUBLIC:: pruned_size_loaded
      PUBLIC:: cmumps_chain_prun_nodes
      PUBLIC:: cmumps_chain_prun_nodes_stats
      PUBLIC:: cmumps_initialize_rhs_bounds
      PUBLIC:: cmumps_propagate_rhs_bounds
      PUBLIC:: cmumps_tree_prun_nodes
      PUBLIC:: cmumps_tree_prun_nodes_stats
      PUBLIC:: cmumps_sol_es_init
      INTEGER(8), POINTER, DIMENSION(:,:) :: size_of_block
      INTEGER(8) :: PRUNED_SIZE_LOADED
      include 'mumps_headers.h'
      CONTAINS
      SUBROUTINE cmumps_sol_es_init(SIZE_OF_BLOCK_ARG, KEEP201)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: keep201
      INTEGER(8), POINTER, DIMENSION(:,:) :: size_of_block_arg
      IF (keep201 > 0) THEN
        size_of_block => size_of_block_arg
      ELSE
        NULLIFY(size_of_block)
      ENDIF
      RETURN
      END SUBROUTINE cmumps_sol_es_init
      SUBROUTINE cmumps_tree_prun_nodes( 
     &     fill,
     &     DAD, NE_STEPS, FRERE, KEEP28,
     &     FILS, STEP, N,
     &     nodes_RHS, nb_nodes_RHS,
     &     TO_PROCESS,
     &     nb_prun_nodes, nb_prun_roots, nb_prun_leaves,
     &     Pruned_List, Pruned_Roots, Pruned_Leaves
     &     )
      IMPLICIT NONE
      LOGICAL, INTENT(IN) :: fill
      INTEGER, INTENT(IN) :: n, keep28
      INTEGER, INTENT(IN) :: DAD(keep28),ne_steps(keep28),frere(keep28)
      INTEGER, INTENT(IN) :: fils(n), step(n)
      INTEGER, INTENT(IN) :: nodes_rhs(keep28),  nb_nodes_rhs
      INTEGER :: nb_prun_nodes
      INTEGER, OPTIONAL, INTENT(INOUT):: pruned_list(nb_prun_nodes)
      INTEGER :: nb_prun_roots
      INTEGER, OPTIONAL, INTENT(INOUT):: pruned_roots(nb_prun_roots)
      INTEGER :: nb_prun_leaves
      INTEGER, OPTIONAL, INTENT(INOUT):: pruned_leaves(nb_prun_leaves)
      LOGICAL :: to_process(keep28) 
      INTEGER :: in, i, ISTEP, tmp, tmpsave
      LOGICAL :: FILS_VISITED
      nb_prun_nodes = 0
      nb_prun_leaves = 0
      to_process(:) = .false.
      DO i = 1, nb_nodes_rhs
         tmp = nodes_rhs(i)
         tmpsave = tmp
         istep = step(tmp)
         DO WHILE(.NOT.to_process(istep))
            to_process(istep) = .true.
            nb_prun_nodes = nb_prun_nodes + 1
            IF(fill) THEN
               pruned_list(nb_prun_nodes) = tmp
            END IF
            in = fils(tmp) 
            DO WHILE(in.GT.0) 
               in = fils(in)
            END DO
            fils_visited = .false.
            IF (in.LT.0) THEN 
             fils_visited = to_process(step(-in))
            ENDIF
            IF ( in.LT.0.and..NOT.fils_visited)
     &            THEN 
               tmp = -in
               istep = step(tmp)
            ELSE 
               IF (in.EQ.0) THEN
                 nb_prun_leaves = nb_prun_leaves + 1
                 IF (fill) THEN
                    pruned_leaves(nb_prun_leaves) = tmp
                 END IF
               ELSE 
                 tmp = -in
                 istep = step(tmp)
               ENDIF
               DO WHILE (tmp.NE.tmpsave) 
                  tmp = abs(frere(istep))
                  IF(tmp.NE.0) THEN 
                     istep = step(tmp) 
                  ELSE 
                     exit
                  END IF
                  IF (.NOT.to_process(istep)) exit
               END DO
            END IF
         END DO
      END DO
      nb_prun_roots = 0
      DO i=1,nb_nodes_rhs
         tmp = nodes_rhs(i)
         istep = step(tmp)
         IF(dad(istep).NE.0) THEN 
            IF(.NOT.to_process(step(dad(istep)))) THEN
               nb_prun_roots = nb_prun_roots + 1
               IF(fill) THEN
                  pruned_roots(nb_prun_roots) = tmp
               END IF
            END IF
         ELSE 
            nb_prun_roots = nb_prun_roots + 1
            IF(fill) THEN
               pruned_roots(nb_prun_roots) = tmp
            END IF          
         END IF
      END DO
      RETURN
      END SUBROUTINE cmumps_tree_prun_nodes
      SUBROUTINE cmumps_chain_prun_nodes(
     &     fill,
     &     DAD, KEEP28,
     &     STEP, N,
     &     nodes_RHS, nb_nodes_RHS,
     &     Pruned_SONS, TO_PROCESS,
     &     nb_prun_nodes,nb_prun_roots, nb_prun_leaves,
     &     Pruned_List, Pruned_Roots, Pruned_Leaves
     &     )
      IMPLICIT NONE
      LOGICAL, INTENT(IN) :: fill
      INTEGER, INTENT(IN) :: n
      INTEGER, INTENT(IN) :: step(n)
      INTEGER, INTENT(IN) :: keep28
      INTEGER, INTENT(IN) :: dad(keep28)
      INTEGER, INTENT(IN) :: nb_nodes_rhs
      INTEGER, INTENT(IN) :: nodes_RHS(nb_nodes_rhs)
      INTEGER :: nb_prun_nodes
      INTEGER, OPTIONAL, INTENT(INOUT):: pruned_list(nb_prun_nodes)
      INTEGER :: nb_prun_roots
      INTEGER, OPTIONAL, INTENT(INOUT):: pruned_roots(nb_prun_roots)
      INTEGER :: nb_prun_leaves
      INTEGER, OPTIONAL, INTENT(INOUT):: pruned_leaves(nb_prun_leaves)
      INTEGER :: pruned_sons(keep28)
      LOGICAL :: to_process(keep28)
      INTEGER :: in, i, istep, tmp
      nb_prun_nodes = 0
      nb_prun_roots = 0
      to_process(:) = .false.
      pruned_sons(:) = -1
      DO i = 1, nb_nodes_rhs
         tmp = nodes_rhs(i)
         istep = step(tmp)
         to_process(istep) = .true.
         IF (pruned_sons(istep) .eq. -1) THEN
            pruned_sons(istep) = 0
            nb_prun_nodes = nb_prun_nodes + 1
            IF(fill) THEN
               pruned_list(nb_prun_nodes) = nodes_rhs(i)
            END IF
            in = nodes_rhs(i)
            in = dad(step(in))
            DO WHILE (in.NE.0)
               to_process(step(in)) = .true.
               IF (pruned_sons(step(in)).eq.-1) THEN 
                  nb_prun_nodes = nb_prun_nodes + 1
                  IF(fill) THEN
                     pruned_list(nb_prun_nodes) = in
                  END IF
                  pruned_sons(step(in)) = 1
                  tmp = in
                  in = dad(step(in))
               ELSE 
                  pruned_sons(step(in)) = pruned_sons(step(in)) + 1
                  GOTO 201
               ENDIF
            ENDDO
            nb_prun_roots = nb_prun_roots +1
            IF(fill) THEN
               pruned_roots(nb_prun_roots) = tmp
            END IF
         ENDIF
  201    CONTINUE
      ENDDO
      nb_prun_leaves = 0
      DO i = 1, nb_nodes_rhs
         tmp = nodes_rhs(i)
         istep = step(tmp)
         IF (pruned_sons(istep).EQ.0) THEN
            nb_prun_leaves = nb_prun_leaves +1
            IF(fill) THEN
              pruned_leaves(nb_prun_leaves) = tmp
            END IF
         END IF
      ENDDO
      RETURN
      END SUBROUTINE cmumps_chain_prun_nodes
      SUBROUTINE cmumps_initialize_rhs_bounds(
     & STEP, N,
     & IRHS_PTR, NBCOL, IRHS_SPARSE, NZ_RHS,
     & JBEG_RHS, PERM_RHS, SIZE_PERM_RHS, K242, K243,
     & UNS_PERM_INV, SIZE_UNS_PERM_INV, K23,
     & RHS_BOUNDS, NSTEPS,
     & nb_sparse, MYID,
     & mode)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: myid, n, nsteps, k242, k243, k23
      INTEGER, INTENT(IN) :: jbeg_rhs, size_perm_rhs, nb_sparse
      INTEGER, INTENT(IN) :: nbcol, nz_rhs, size_uns_perm_inv
      INTEGER, INTENT(IN) :: step(n), perm_rhs(size_perm_rhs)
      INTEGER, INTENT(IN) :: irhs_ptr(nbcol+1),irhs_sparse(nz_rhs)
      INTEGER, INTENT(IN) :: uns_perm_inv(size_uns_perm_inv)
      INTEGER, INTENT(INOUT):: rhs_bounds(2*nsteps)
      INTEGER, INTENT(IN) :: mode 
      INTEGER :: i, icol, jptr, j, jam1, node, bound
      rhs_bounds = 0
      icol = 0 
      DO i = 1, nbcol
        IF ( (irhs_ptr(i+1)-irhs_ptr(i)).EQ.0) cycle
        icol = icol + 1
        bound = icol - mod(icol, nb_sparse) + 1
        IF(mod(icol, nb_sparse).EQ.0) bound = bound - nb_sparse
        IF(mode.EQ.0) THEN 
          IF ((k242.NE.0).OR.(k243.NE.0)) THEN
            jam1 = perm_rhs(jbeg_rhs+i-1)
          ELSE
            jam1 = jbeg_rhs+i-1
          ENDIF
          node = abs(step(jam1))
          IF(rhs_bounds(2*node - 1).EQ.0) THEN 
            rhs_bounds(2*node - 1) = bound                 
            rhs_bounds(2*node)     = bound + nb_sparse - 1 
          ELSE
            rhs_bounds(2*node) = bound + nb_sparse - 1
          END IF
        ELSE  
          DO jptr = irhs_ptr(i), irhs_ptr(i+1)-1
            j = irhs_sparse(jptr)
            IF ( mode .EQ. 1 ) THEN
              IF (k23.NE.0) j = uns_perm_inv(j)
            ENDIF
            node = abs(step(j))
            IF(rhs_bounds(2*node - 1).EQ.0) THEN
              rhs_bounds(2*node - 1) = bound
              rhs_bounds(2*node)     = bound + nb_sparse - 1
            ELSE
              rhs_bounds(2*node) = bound + nb_sparse - 1
            END IF
          END DO
        END IF
      END DO
      RETURN
      END SUBROUTINE cmumps_initialize_rhs_bounds
      SUBROUTINE cmumps_propagate_rhs_bounds(
     & pruned_leaves, nb_pruned_leaves,
     & STEP, N, Pruned_SONS,
     & DAD, RHS_BOUNDS, NSTEPS,
     & MYID, COMM, KEEP485,
     & IW, LIW, PTRIST, KIXSZ,OOC_FCT_LOC, PHASE, LDLT, K38)
      IMPLICIT NONE
      include 'mpif.h'
      include 'mumps_headers.h'
      INTEGER, INTENT(IN) :: nb_pruned_leaves, n, nsteps
      INTEGER, INTENT(IN) :: step(n), DAD(nsteps), pruned_sons(nsteps)
      INTEGER, INTENT(IN) :: myid, comm, keep485
      INTEGER, INTENT(IN) :: pruned_leaves(nb_pruned_leaves)
      INTEGER, INTENT(IN) :: liw, iw(liw), ptrist(nsteps)
      INTEGER, INTENT(IN) :: kixsz, ooc_fct_loc, phase, ldlt, k38
      INTEGER, INTENT(INOUT):: rhs_bounds(2*nsteps)
      INTEGER :: i, node, father, size_pool, next_size_pool
      INTEGER :: ierr
      INTEGER, ALLOCATABLE, DIMENSION(:) :: pool, nbsons
      ALLOCATE(pool(nb_pruned_leaves),
     &         nbsons(nsteps),
     &         stat=ierr)
      IF (ierr.NE.0) THEN
         WRITE(6,*)'Allocation problem in CMUMPS_PROPAGATE_RHS_BOUNDS'
         CALL mumps_abort()
      END IF
      size_pool = nb_pruned_leaves
      pool = pruned_leaves
      nbsons = pruned_sons
      DO WHILE (size_pool.ne.0)
        next_size_pool =0 
        DO i=1, size_pool 
          node = step(pool(i)) 
          IF (dad(node).NE.0) THEN
            father = step(dad(node))
            nbsons(father) = nbsons(father)-1
            IF (rhs_bounds(2*father-1).EQ.0) THEN
              rhs_bounds(2*father-1) = rhs_bounds(2*node-1)
              rhs_bounds(2*father)   = rhs_bounds(2*node)
            ELSE
              rhs_bounds(2*father-1) = min(rhs_bounds(2*father-1),
     &                                     rhs_bounds(2*node-1))
              rhs_bounds(2*father) = max(rhs_bounds(2*father),
     &                                     rhs_bounds(2*node))
            END IF
            IF(nbsons(father).EQ.0) THEN 
              next_size_pool = next_size_pool+1
              pool(next_size_pool) = dad(node)
            END IF
          END IF
        END DO
        size_pool = next_size_pool 
      END DO
      DEALLOCATE(pool, nbsons)
      RETURN
      END SUBROUTINE cmumps_propagate_rhs_bounds
      INTEGER(8) FUNCTION cmumps_local_factor_size(IW,LIW,PTR,
     &                                 PHASE, LDLT, IS_ROOT)
        INTEGER, INTENT(IN) :: liw, PTR, phase, ldlt
        INTEGER, INTENT(IN) :: iw(liw)
        LOGICAL, INTENT(IN) :: is_root
        INTEGER(8) :: ncb, nelim, liell, npiv, nrow
        ncb   = int(iw(ptr),8)     
        nelim = int(iw(ptr+1),8)   
        nrow  = int(iw(ptr+2),8)
        npiv  = int(iw(ptr+3),8)
        liell = npiv + ncb
        IF (is_root) THEN
          cmumps_local_factor_size = int(iw(ptr+1),8) *     
     &                               int(iw(ptr+2),8) / 2_8 
          RETURN
        ENDIF
        IF (ncb.GE.0_8) THEN 
          IF (phase.EQ.0   
     &      .OR. (phase.EQ.1.AND.ldlt.NE.0) 
     &  ) THEN
            cmumps_local_factor_size =
     &            npiv*(npiv-1_8)/2_8 + (nrow-npiv)*npiv
          ELSE
            cmumps_local_factor_size =
     &      npiv*(npiv+1_8)/2_8 + (liell-npiv)*npiv
          ENDIF
        ELSE
          cmumps_local_factor_size =
     &      -ncb*nelim
        END IF
      RETURN
      END FUNCTION cmumps_local_factor_size
      INTEGER(8) FUNCTION cmumps_local_factor_size_blr(IW,LIW,PTR,
     &                                LRSTATUS, IWHANDLER,
     &                                PHASE, LDLT, IS_ROOT)
      USE cmumps_lr_data_m
      USE cmumps_lr_type
        INTEGER, INTENT(IN) :: liw, ptr, phase, ldlt
        INTEGER, INTENT(IN) :: lrstatus, iwhandler
        INTEGER, INTENT(IN) :: iw(liw)
        LOGICAL, INTENT(IN) :: is_root
        INTEGER(8) :: ncb, nelim, liell, npiv, nrow, factor_size
        INTEGER :: nb_panels, ipanel, loru, iblock
        LOGICAL :: lr_activated
        TYPE(lrb_type), POINTER, DIMENSION(:) :: lrb_panel
        ncb   = int(iw(ptr),8)     
        nelim = int(iw(ptr+1),8)   
        nrow  = int(iw(ptr+2),8)
        npiv  = int(iw(ptr+3),8)
        liell = npiv + ncb
        lr_activated=(lrstatus.GE.2)
        IF (lr_activated) THEN
          factor_size = 0_8
          CALL cmumps_blr_retrieve_nb_panels(iwhandler, nb_panels)
          IF (ldlt.EQ.0) THEN 
            loru = phase
          ELSE
            loru = 0
          ENDIF
          DO ipanel=1,nb_panels
            IF (is_root.AND.ipanel.EQ.nb_panels) THEN
              cycle
            ENDIF
            IF (cmumps_blr_empty_panel_loru(iwhandler, loru, ipanel))
     &        THEN
              cycle
            ENDIF
             CALL cmumps_blr_retrieve_panel_loru(iwhandler, loru,
     &                       ipanel, lrb_panel)
            IF (size(lrb_panel).GT.0) THEN
              IF (phase.EQ.0) THEN
                factor_size = factor_size + 
     &            int(lrb_panel(1)%N,8)*(int(lrb_panel(1)%N,8)-1_8)/2_8
              ELSE
                factor_size = factor_size + 
     &           int(lrb_panel(1)%N,8)*(int(lrb_panel(1)%N,8)+1_8)/2_8
              ENDIF
            ENDIF
            DO iblock=1,size(lrb_panel)
              IF (lrb_panel(iblock)%ISLR) THEN
                factor_size = factor_size + int(lrb_panel(iblock)%K,8)*
     &            int(lrb_panel(iblock)%M+lrb_panel(iblock)%M,8)
              ELSE
                factor_size = factor_size +
     &             int(lrb_panel(iblock)%M*lrb_panel(iblock)%N,8)
              ENDIF
            ENDDO
          ENDDO
          cmumps_local_factor_size_blr = factor_size
        ELSE
          cmumps_local_factor_size_blr = 
     &      cmumps_local_factor_size(iw, liw, ptr, phase, ldlt, is_root)
        ENDIF
      RETURN
      END FUNCTION cmumps_local_factor_size_blr
      SUBROUTINE cmumps_tree_prun_nodes_stats(MYID, N, KEEP28, KEEP201,
     &           FR_FACT,
     &           STEP, Pruned_List, nb_prun_nodes, OOC_FCT_TYPE_LOC)
      INTEGER, intent(in) :: keep28, keep201, ooc_fct_type_loc, myid, n
      INTEGER(8), intent(in) :: fr_fact
      INTEGER, intent(in) :: nb_prun_nodes
      INTEGER, intent(in) :: pruned_list(nb_prun_nodes)
      INTEGER, intent(in) :: step(n)
      INTEGER i, istep
      INTEGER(8) :: pruned_size
      if (keep201 .GT. 0) then
        pruned_size = 0_8
        DO i = 1, nb_prun_nodes
          istep = step(pruned_list(i))
          pruned_size = pruned_size + size_of_block
     &                  (istep, ooc_fct_type_loc)
        ENDDO
        pruned_size_loaded = pruned_size_loaded +pruned_size
      ENDIF
      RETURN
      END SUBROUTINE cmumps_tree_prun_nodes_stats
      SUBROUTINE cmumps_chain_prun_nodes_stats
     &                (myid, n, keep28, keep201, keep485, fr_fact,
     &                step, pruned_list, nb_prun_nodes, ooc_fct_type_loc
     & )
      IMPLICIT NONE
      INTEGER, intent(in) :: keep28, keep201, ooc_fct_type_loc, n,
     &                       keep485
      INTEGER(8), intent(in) :: fr_fact
      INTEGER, intent(in) :: nb_prun_nodes, myid
      INTEGER, intent(in) :: pruned_list(nb_prun_nodes)
      INTEGER, intent(in) :: step(n)
      include 'mpif.h'
      INTEGER I, istep
      INTEGER(8) :: pruned_size
      pruned_size = 0_8
      DO i = 1, nb_prun_nodes
        istep = step(pruned_list(i))
        IF (keep201 .GT. 0) THEN
            pruned_size = pruned_size + size_of_block
     &                    (istep, ooc_fct_type_loc)
        ENDIF
      ENDDO
      IF (keep201.GT.0) THEN
        IF (fr_fact .NE. 0_8) THEN
          pruned_size_loaded = pruned_size_loaded +pruned_size
        ENDIF
      ENDIF
      RETURN
      END SUBROUTINE cmumps_chain_prun_nodes_stats
      END MODULE cmumps_sol_es
      SUBROUTINE cmumps_permute_rhs_gs
     &          (lp, lpok, prokg, mpg, perm_strat, 
     &           sym_perm, n, nrhs,
     &           irhs_ptr, size_irhs_ptr, 
     &           irhs_sparse, nzrhs, 
     &           perm_rhs, ierr
     &         )
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: LP, MPG, PERM_STRAT, N, NRHS, 
     &                       SIZE_IRHS_PTR,
     &                       nzrhs
      LOGICAL, INTENT(IN) :: LPOK, PROKG
      INTEGER, INTENT(IN) :: SYM_PERM(N)
      INTEGER, INTENT(IN) :: IRHS_PTR(SIZE_IRHS_PTR)
      INTEGER, INTENT(IN) :: IRHS_SPARSE(NZRHS)
      INTEGER, INTENT(OUT) :: PERM_RHS(NRHS)
      INTEGER, INTENT(OUT) :: IERR
      INTEGER :: I,J,K, POSINPERMRHS, JJ,
     &           kpos
      INTEGER, ALLOCATABLE :: ROW_REFINDEX(:)
      ierr = 0
      IF ((perm_strat.NE.-1).AND.(perm_strat.NE.1)) THEN
       ierr=-1
       IF (lpok)
     & WRITE(lp,*) " INTERNAL ERROR -1 in ",
     &       " CMUMPS_PERMUTE_RHS_GS, PERM_STRAT =", perm_strat, 
     &       " is out of range "
       RETURN
      ENDIF
      IF (perm_strat.EQ.-1) THEN
       DO i=1,nrhs
        perm_rhs(i) = i
       END DO
       GOTO 490
      ENDIF
      ALLOCATE(row_refindex(nrhs), stat=ierr)
      IF (ierr.GT.0) THEN
       ierr=-1
       IF (lpok) THEN
          WRITE(lp,*) " ERROR -2 : ", 
     &         " ALLOCATE IN CMUMPS_PERMUTE_RHS_GS OF SIZE :",
     &         nrhs
       ENDIF
       RETURN
      ENDIF
      DO i=1,nrhs
        IF (irhs_ptr(i+1)-irhs_ptr(i).LE.0) THEN
          ierr =  1
          IF (i.EQ.1) THEN
            row_refindex(i) = irhs_sparse(irhs_ptr(i))
          ELSE
            row_refindex(i) = row_refindex(i-1)
          ENDIF
        ELSE
          row_refindex(i) = irhs_sparse(irhs_ptr(i))
        ENDIF
      END DO
      posinpermrhs = 0
      DO i=1,nrhs
       kpos = n+1 
       jj   = 0   
       DO j=1,nrhs
        k = row_refindex(j)
        IF (k.LE.0) cycle 
        IF (sym_perm(k).LT.kpos) THEN
         kpos = sym_perm(k)
         jj   = j
        ENDIF
       END DO
       IF (jj.EQ.0) THEN
         ierr = -3 
         IF (lpok)
     &   WRITE(lp,*) " INTERNAL ERROR -3 in ",
     &       " CMUMPS_PERMUTE_RHS_GS "
         GOTO 500
       ENDIF
       posinpermrhs           = posinpermrhs + 1
       perm_rhs(posinpermrhs) = jj
       row_refindex(jj)       = -row_refindex(jj)
      END DO
      IF (posinpermrhs.NE.nrhs) THEN
         IF (lpok)
     &   WRITE(lp,*) " INTERNAL ERROR -4 in ",
     &       " CMUMPS_PERMUTE_RHS_GS ", maxval(row_refindex)
         ierr = -4
         GOTO 500
      ENDIF
  490 CONTINUE
 500  CONTINUE
      IF (allocated(row_refindex)) DEALLOCATE(row_refindex)
      END SUBROUTINE cmumps_permute_rhs_gs
      SUBROUTINE cmumps_permute_rhs_am1
     &          (perm_strat, sym_perm,
     &           irhs_ptr, nhrs,
     &           perm_rhs, sizeperm, ierr
     &         )
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: PERM_STRAT, NHRS, SIZEPERM
      INTEGER, INTENT(IN) :: SYM_PERM(SIZEPERM)
      INTEGER, INTENT(IN) :: IRHS_PTR(NHRS)
      INTEGER, INTENT(OUT):: IERR   
      INTEGER, INTENT(OUT):: PERM_RHS(SIZEPERM)
      DOUBLE PRECISION :: RAND_NUM
      INTEGER  I, J, STRAT
      ierr = 0
      strat = perm_strat 
      IF( (strat.NE.-3).AND.
     &    (strat.NE.-2).AND.
     &    (strat.NE.-1).AND.
     &    (strat.NE. 1).AND.
     &    (strat.NE. 2).AND.
     &    (strat.NE. 6) ) THEN
        WRITE(*,*)"Warning: incorrect value for the RHS permutation; ",
     &            "defaulting to post-order"
        strat = 1
      END IF
      IF (strat .EQ. -3) THEN
         perm_rhs(1:sizeperm)=0
         DO i=1, sizeperm 
           CALL random_number(rand_num) 
           rand_num = rand_num*dble(sizeperm) 
           j = ceiling(rand_num) 
           DO WHILE (perm_rhs(j).NE.0) 
             CALL random_number(rand_num)
             rand_num = rand_num*dble(sizeperm)
             j = ceiling(rand_num)
           ENDDO
           perm_rhs(j)=i
         ENDDO
      ELSEIF (strat .EQ. -2) THEN
         DO i=1, sizeperm
            perm_rhs(sizeperm -i +1) = i
         ENDDO
      ELSEIF (strat .EQ. -1) THEN
         DO i=1, sizeperm
            perm_rhs(i) = i
         ENDDO
      ELSEIF (strat .EQ.  1) THEN
         DO i=1, sizeperm
            perm_rhs(sym_perm(i)) = i
         ENDDO
      ELSEIF (strat .EQ.  2) THEN
         DO i=1, sizeperm
            perm_rhs(sizeperm-sym_perm(i)+1) = i
         ENDDO
      ENDIF
      END SUBROUTINE cmumps_permute_rhs_am1
      SUBROUTINE cmumps_interleave_rhs_am1(
     &  PERM_RHS, SIZE_PERM,
     &  IPTR_WORKING, SIZE_IPTR_WORKING, WORKING, SIZE_WORKING,
     &  IRHS_PTR,
     &  STEP, SYM_PERM, N, NBRHS,
     &  PROCNODE, NSTEPS, SLAVEF, KEEP199,
     &  behaviour_L0, reorder, n_select, PROKG, MPG 
     &  )
      IMPLICIT NONE
      INTEGER, INTENT(IN) ::  SIZE_PERM,
     &                        size_iptr_working,
     &                        iptr_working(size_iptr_working),
     &                        size_working,
     &                        working(size_working),
     &                        n,
     &                        irhs_ptr(n+1),
     &                        step(n),
     &                        sym_perm(n),
     &                        nbrhs,
     &                        nsteps,
     &                        procnode(nsteps),
     &                        slavef, keep199,
     &                        n_select, mpg
      LOGICAL, INTENT(IN) :: behaviour_L0,
     &                        reorder, prokg
      INTEGER, INTENT(INOUT) :: PERM_RHS(SIZE_PERM)
      INTEGER :: I, J, K,
     &           entry,            
     &           node,             
     &           size_perm_working,
     &           nb_non_empty,     
     &           to_be_found,      
     &           posintmprhs,      
     &           selected,         
     &           local_selected,   
     &           current_proc,     
     &           nprocs,           
     &           n_pass,           
     &           pass,             
     &           nblocks,          
     &           n_select_loc,     
     &           ierr
      INTEGER, ALLOCATABLE, DIMENSION(:) :: TMP_RHS,    
     &                                      PTR_PROCS,  
     &                                      LOAD_PROCS,  
     &                                      IPTR_PERM_WORKING,
     &                                      PERM_WORKING,
     &                                      MYTYPENODE,
     &                                      perm_po
      LOGICAL, ALLOCATABLE, DIMENSION(:) :: USED
      LOGICAL :: allow_above_L0
      INTEGER, EXTERNAL :: MUMPS_TYPENODE_ROUGH
      NPROCS = size_iptr_working - 1
      ALLOCATE(tmp_rhs(size_perm),
     &         ptr_procs(nprocs),
     &         load_procs(nprocs),
     &         used(size_perm),
     &         iptr_perm_working(nprocs+1),
     &         mytypenode(nsteps),
     &         stat=ierr)
      IF(ierr.GT.0) THEN
        WRITE(*,*)'Allocation error in CMUMPS_INTERLEAVE_RHS_AM1'
        CALL mumps_abort()
      END IF
      DO i=1, nsteps
        mytypenode(i) = mumps_typenode_rough( procnode(i), keep199 )
      ENDDO
      nb_non_empty = 0
      DO i=1,size_perm
        IF(irhs_ptr(i+1)-irhs_ptr(i).NE.0) THEN 
          nb_non_empty = nb_non_empty + 1
        END IF
      END DO
      k = 0
      iptr_perm_working(1)=1
      DO i=1,nprocs
        used = .false.
        DO j=iptr_working(i),iptr_working(i+1)-1
          used(working(j)) = .true.        
        END DO
        DO j=1,n
          IF (used(abs(step(perm_rhs(j)))).AND.
     &      ((irhs_ptr(perm_rhs(j)+1)-irhs_ptr(perm_rhs(j))).NE.0))
     &    THEN
            k = k + 1
          END IF
        END DO
        iptr_perm_working(i+1) = k+1
      END DO
      size_perm_working = k
      ALLOCATE(perm_working(size_perm_working),
     &         stat=ierr)
      IF(ierr.GT.0) THEN
        WRITE(*,*)'Allocation error in CMUMPS_INTERLEAVE_RHS_AM1'
        CALL mumps_abort()
      END IF
      k = 0
      DO i=1,nprocs
        used = .false.
        DO j=iptr_working(i),iptr_working(i+1)-1
          used(working(j)) = .true.        
        END DO
        DO j=1,n
          IF (used(abs(step(perm_rhs(j)))).AND.
     &      ((irhs_ptr(perm_rhs(j)+1)-irhs_ptr(perm_rhs(j))).NE.0)) 
     &    THEN
            k = k + 1
            perm_working(k) = perm_rhs(j)
          END IF
        END DO
      END DO      
      IF(behaviour_l0) THEN
        n_pass = 2
        allow_above_l0 = .false.
        to_be_found = 0
        DO i=1,size_perm
          IF((mytypenode(abs(step(i))).LE.1).AND. 
     &    (irhs_ptr(i+1)-irhs_ptr(i).NE.0))       
     &    THEN
            to_be_found = to_be_found + 1
          END IF
        END DO
      ELSE
        n_pass = 1
        allow_above_l0 = .true.
        to_be_found = nb_non_empty
      END IF
      ptr_procs(1:nprocs) = iptr_perm_working(1:nprocs)
      load_procs = 0
      used = .false.
      current_proc = 1
      n_select_loc = n_select
      IF (n_select_loc.LE.0) THEN
       n_select_loc = 1
      ENDIF
      posintmprhs = 0
      DO pass=1,n_pass
        selected = 0
        DO WHILE(selected.LT.to_be_found)
          local_selected = 0
          DO WHILE(local_selected.LT.n_select_loc)
            IF(ptr_procs(current_proc).EQ.
     &        iptr_perm_working(current_proc+1))
     &      THEN 
              EXIT
            ELSE 
              entry = perm_working(ptr_procs(current_proc))
              node  = abs(step(entry))
              IF(.NOT.used(entry)) THEN
                IF(allow_above_l0.OR.(mytypenode(node).LE.1)) THEN
                  used(entry) = .true.
                  selected = selected + 1
                  local_selected = local_selected + 1
                  posintmprhs = posintmprhs + 1
                  tmp_rhs(posintmprhs) = entry
                  IF(selected.EQ.to_be_found) EXIT
                END IF
              END IF
              ptr_procs(current_proc) = ptr_procs(current_proc) + 1
            END IF
          END DO
          current_proc = mod(current_proc,nprocs)+1
        END DO
        to_be_found = nb_non_empty - to_be_found
        allow_above_l0 = .true.
        ptr_procs(1:nprocs) = iptr_perm_working(1:nprocs)
      END DO
      DO i=1,size_perm
        IF(irhs_ptr(perm_rhs(i)+1)-irhs_ptr(perm_rhs(i)).EQ.0) THEN
          posintmprhs = posintmprhs+1
          tmp_rhs(posintmprhs) = perm_rhs(i)
          IF(posintmprhs.EQ.size_perm) EXIT
        END IF
      END DO
      IF(reorder) THEN
        posintmprhs = 0 
        ALLOCATE(perm_po(n),stat=ierr)
        IF(ierr.GT.0) THEN
          WRITE(*,*)'Allocation error in INTERLEAVE_RHS_AM1'
          CALL mumps_abort()
        END IF
        DO j=1,n
          perm_po(sym_perm(j))=j
        END DO
        nblocks = n/nbrhs
        DO i = 1, nblocks
          used = .false.
          DO j=1, nbrhs
            used(tmp_rhs(nbrhs*(i-1)+j))=.true.
          END DO
          DO j=1,n
            IF(used(perm_po(j))) THEN
              posintmprhs = posintmprhs + 1
              perm_rhs(posintmprhs) = perm_po(j)
            END IF
          END DO
        END DO
        IF(mod(n,nbrhs).NE.0) THEN
          used = .false.
          DO j=1, mod(n,nbrhs)
            used(tmp_rhs(nbrhs*nblocks+j))=.true.
          END DO
          DO j=1,n
            IF(used(perm_po(j))) THEN
              posintmprhs = posintmprhs + 1
              perm_rhs(posintmprhs) = perm_po(j)
            END IF
          END DO
        END IF
        DEALLOCATE(perm_po)
      ELSE
        perm_rhs = tmp_rhs
      END IF
      DEALLOCATE(tmp_rhs,
     &           ptr_procs,
     &           load_procs,
     &           used,
     &           iptr_perm_working,
     &           perm_working,
     &           mytypenode)
      RETURN
      END SUBROUTINE cmumps_interleave_rhs_am1