zsol_driver.F

Go to the documentation of this file.

C
C  This file is part of MUMPS 5.5.1, released
C  on Tue Jul 12 13:17:24 UTC 2022
C
C
C  Copyright 1991-2022 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
C  Mumps Technologies, University of Bordeaux.
C
C  This version of MUMPS is provided to you free of charge. It is
C  released under the CeCILL-C license 
C  (see doc/CeCILL-C_V1-en.txt, doc/CeCILL-C_V1-fr.txt, and
C  https://cecill.info/licences/Licence_CeCILL-C_V1-en.html)
C
      SUBROUTINE zmumps_solve_driver(id)
      USE zmumps_struc_def
      USE zmumps_sol_es
C     Lock Initialization (_LI) and Desruction (_LD)
      USE zmumps_sol_l0omp_m, ONLY: zmumps_sol_l0omp_li,
     &                              zmumps_sol_l0omp_ld
C
C  Purpose
C  =======
C
C  Performs solution phase (solve), Iterative Refinements
C  and Error analysis.
C
C
C
C
      USE zmumps_buf
      USE zmumps_ooc
      USE mumps_memory_mod
      USE zmumps_lr_data_m, only : zmumps_blr_struc_to_mod
     &                           , zmumps_blr_mod_to_struc
      USE mumps_front_data_mgt_m, only : mumps_fdm_struc_to_mod,
     &                                   mumps_fdm_mod_to_struc
      USE zmumps_save_restore
!$    USE OMP_LIB
      IMPLICIT NONE
C     -------------------
C     Explicit interfaces
C     -------------------
      INTERFACE
      SUBROUTINE zmumps_size_in_struct( id, NB_INT,NB_CMPLX,NB_CHAR )
      USE zmumps_struc_def
      TYPE (ZMUMPS_STRUC) :: id
      INTEGER(8)        :: NB_INT,NB_CMPLX,NB_CHAR
      END SUBROUTINE zmumps_size_in_struct
      SUBROUTINE zmumps_check_dense_rhs
     &(idrhs, idinfo, idn, idnrhs, idlrhs)
      COMPLEX(kind=8), DIMENSION(:), POINTER :: idRHS
      INTEGER, intent(in)    :: idN, idNRHS, idLRHS
      INTEGER, intent(inout) :: idINFO(:)
      END SUBROUTINE zmumps_check_dense_rhs
      END INTERFACE
C
      include 'mpif.h'
      include 'mumps_headers.h'
      include 'mumps_tags.h'
#if defined(V_T)
      include 'VT.inc'
#endif
      INTEGER :: STATUS(MPI_STATUS_SIZE)
      INTEGER :: IERR
      INTEGER, PARAMETER :: MASTER = 0
c
C  Parameters
C  ==========
C
      TYPE (ZMUMPS_STRUC), TARGET :: id
C
C  Local variables
C  ===============
C
      INTEGER MP,LP, MPG
      LOGICAL PROK, PROKG, LPOK
      INTEGER MTYPE, ICNTL21
      LOGICAL LSCAL, POSTPros, GIVSOL
      INTEGER ICNTL10, ICNTL11
      INTEGER I,IPERM,K,JPERM, J, II, IZ2
      INTEGER IZ, NZ_THIS_BLOCK, PJ
C     pointers in IS
      INTEGER LIW
C     pointers in id%S
      INTEGER(8) :: LA, LA_PASSED
      INTEGER LIW_PASSED
      INTEGER(8) :: LWCB8_MIN, LWCB8, LWCB8_SOL_C
C     buffer sizes
      INTEGER ZMUMPS_LBUF, ZMUMPS_LBUF_INT
      INTEGER(8) :: ZMUMPS_LBUF_8
      INTEGER :: LBUFR, LBUFR_BYTES
      INTEGER :: MSG_MAX_BYTES_SOLVE, MSG_MAX_BYTES_GTHRSOL
      INTEGER(8) :: MSG_MAX_BYTES_SOLVE8
C     reception buffer
      INTEGER, ALLOCATABLE, DIMENSION(:) :: BUFR
C     null space
      INTEGER IBEG_ROOT_DEF, IEND_ROOT_DEF,
     &        IBEG_GLOB_DEF, IEND_GLOB_DEF,
     &        IROOT_DEF_RHS_COL1
C
      INTEGER NITREF, NOITER, SOLVET, KASE
C     Meaningful only with tree pruning and sparse RHS
      LOGICAL INTERLEAVE_PAR, DO_PERMUTE_RHS
C     true if ZMUMPS_SOL_C called during postprocessing
      LOGICAL FROM_PP
C
C     TIMINGS
      DOUBLE PRECISION TIMEIT, TIMEEA, TIMEEA1, TIMELCOND
      DOUBLE PRECISION TIME3
      DOUBLE PRECISION TIMEC1,TIMEC2
      DOUBLE PRECISION TIMEGATHER1,TIMEGATHER2
      DOUBLE PRECISION TIMESCATTER1,TIMESCATTER2
      DOUBLE PRECISION TIMECOPYSCALE1,TIMECOPYSCALE2
C     ------------------------------------------
C     Declarations related to exploit sparsity
C     ------------------------------------------
      INTEGER     :: NRHS_NONEMPTY
      INTEGER     :: STRAT_PERMAM1
      LOGICAL     :: DO_NULL_PIV
      INTEGER, DIMENSION(:), POINTER :: IRHS_PTR_COPY
      INTEGER, DIMENSION(:), POINTER :: IRHS_SPARSE_COPY
      COMPLEX(kind=8), DIMENSION(:), POINTER :: RHS_SPARSE_COPY
      LOGICAL IRHS_SPARSE_COPY_ALLOCATED, IRHS_PTR_COPY_ALLOCATED,
     &        RHS_SPARSE_COPY_ALLOCATED
C
      INTEGER, DIMENSION(:), ALLOCATABLE :: MAP_RHS_loc
      INTEGER, DIMENSION(:), POINTER :: IRHS_loc_PTR
      LOGICAL :: IRHS_loc_PTR_allocated
      COMPLEX(kind=8), DIMENSION(:), POINTER :: idRHS_loc
      INTEGER(8) :: DIFF_SOL_loc_RHS_loc
      INTEGER(8) :: RHS_loc_size, RHS_loc_shift
      INTEGER(8)  :: NBT
      INTEGER     :: NBCOL, COLSIZE, JBEG_RHS, JEND_RHS, JBEG_NEW,
     &               NBCOL_INBLOC, IPOS, IPOSRHSCOMP
      INTEGER, DIMENSION(:), ALLOCATABLE :: PERM_RHS
      INTEGER, DIMENSION(:), POINTER :: PTR_POSINRHSCOMP_FWD,
     &                                  PTR_POSINRHSCOMP_BWD
      COMPLEX(kind=8), DIMENSION(:), POINTER :: PTR_RHS
      INTEGER :: SIZE_IPTR_WORKING, SIZE_WORKING
C     NRHS_NONEMPTY: holds
C         either the original number of RHS (id%NRHS defined on host)
C         or, when the RHS is sparse, it holds the
C         number of non empty columns.
C         it is computed on master and is
C         then broadcasted on all processes.
C     IRHS_PTR_COPY holds a compressed local copy of IRHS_PTR (or points
C              on the master to id%IRHS_PTR if no permutation requested)
C     IRHS_SPARSE_COPY might be allocated or might also point to
C         id%IRHS_SPARSE. To test if we can deallocate it we trace
C         with IRHS_SPARSE_COPY_ALLOCATED when it was effectively
C         allocated.
C     NBCOL_INBLOC  total nb columns to process in this block
C     JBEG_RHS global ptr for starting column requested for this block
C     JEND_RHS global ptr for end column_number requested for this block
C     PERM_RHS -- Permutation of RHS computed on master and broadcasted
C         on all procs (of size id%NRHS orginal)
C         PERM_RHS(k) = i means that i is the kth column to be processed
C         Note that PERM_RHS will be used also in case of interleaving
C     ------------------------------------
      INTEGER :: NOMP
      COMPLEX(kind=8) ONE
      COMPLEX(kind=8) ZERO
      parameter( one = (1.0d0,0.0d0) )
      parameter( zero = (0.0d0,0.0d0) )
      DOUBLE PRECISION RZERO, RONE
      parameter( rzero = 0.0d0, rone = 1.0d0 )
C
C     RHS_IR is internal to ZMUMPS and used for iterative refinement
C     or the error analysis section. It either points to the user's
C     RHS (on the host when the solution is centralized or the RHS
C     is dense), or is a workarray allocated inside this routine
C     of size N.
      COMPLEX(kind=8), DIMENSION(:), POINTER :: RHS_IR
      COMPLEX(kind=8), DIMENSION(:), POINTER :: WORK_WCB
      COMPLEX(kind=8), DIMENSION(:), POINTER :: PTR_RHS_ROOT
      INTEGER(8) :: LPTR_RHS_ROOT
C
C  Local workarrays that will be dynamically allocated
C
      COMPLEX(kind=8), ALLOCATABLE :: SAVERHS(:), C_RW1(:),
     &                                 C_RW2(:),
     &                                 SRW3(:), C_Y(:),
     &                                 C_W(:)
      INTEGER :: LCWORK
      COMPLEX(kind=8), ALLOCATABLE :: CWORK(:)
      INTEGER, ALLOCATABLE :: MAP_RHS(:)
      DOUBLE PRECISION, ALLOCATABLE :: R_Y(:), D(:)
      DOUBLE PRECISION, ALLOCATABLE :: R_W(:)
C     The 2 following workarrays are temporary local
C     arrays only used for distributed matrix input
C     (KEEP(54) .NE. 0).
      DOUBLE PRECISION,    ALLOCATABLE, DIMENSION(:) :: R_LOCWK54
      COMPLEX(kind=8), ALLOCATABLE, DIMENSION(:) :: C_LOCWK54
      INTEGER   :: NBENT_RHSCOMP, NB_FS_RHSCOMP_F,
     &             NB_FS_RHSCOMP_TOT
      INTEGER, DIMENSION(:), ALLOCATABLE :: UNS_PERM_INV
      LOGICAL :: UNS_PERM_INV_NEEDED_INMAINLOOP,
     &           UNS_PERM_INV_NEEDED_BEFMAINLOOP
      INTEGER LIWK_SOLVE, LIWCB
      INTEGER, ALLOCATABLE :: IW1(:), IWK_SOLVE(:), IWCB(:)
      INTEGER ::  LIWK_PTRACB
      INTEGER(8), ALLOCATABLE :: PTRACB(:)
C
C  Parameters arising from the structure
C
      INTEGER(8)       :: MAXS
      DOUBLE PRECISION, DIMENSION(:), POINTER :: CNTL
      INTEGER, DIMENSION (:), POINTER :: KEEP,ICNTL,INFO
      INTEGER(8), DIMENSION (:), POINTER :: KEEP8
      INTEGER, DIMENSION (:), POINTER :: IS
      DOUBLE PRECISION, DIMENSION(:),POINTER::   RINFOG
C     ===============================================================
C     SCALING issues:
C       When scaling was performed
C       RHS holds the solution of the scaled system
C       The unscaled second member (b0) was given
C       then we have to scale both rhs adn solution:
C        A(sca) = LU  = D1*A*D2 , with D2 = COLSCA
C                                      D1 = ROWSCA
C        --------------
C        CASE OF A X =B
C        --------------
C         (ICNTL(9)=1 or MTYPE=1)
C           A*x0 = b0
C           b(sca) = D1 * b0 = ROWSCA*S(ISTW3)
C           A(sca) [(D2) **(-1)] x0 = b(sca)
C           so the computed solution by Check y0 of LU *y0 = b(sca)
C           is : y0 =[(D2) **(-1)] x0 and so x0= D2*y0 is modified
C        --------------
C        CASE OF AT X =B
C        --------------
C         (ICNTL(9).NE.1 or MTYPE=0)
C           A(sca) = LU  = D1*A*D2
C           AT*x0 = b0 => D2ATD1 D1-1 x0 = D2b0
C           b(sca) = D2 * b0 = COLSCA*S(ISTW3)
C           A(sca)T [(D1) **(-1)] x0 = b(sca)
C           so the computed solution by Check y0 of LU *y0 = b(sca)
C           is : y0 =[(D1) **(-1)] x0 and so x0= D1*y0 is modified
C
C      In case of distributed RHS we need
C      scaling information on each processor
C
          type scaling_data_t
            sequence
            DOUBLE PRECISION, dimension(:), pointer :: SCALING
            DOUBLE PRECISION, dimension(:), pointer :: SCALING_LOC
          end type scaling_data_t
          type (scaling_data_t) :: scaling_data_sol, scaling_data_dr
C      To scale on the fly during GATHER SOLUTION
          DOUBLE PRECISION, DIMENSION(:), POINTER :: PT_SCALING
          DOUBLE PRECISION, TARGET                :: Dummy_SCAL(1)
C
C  ==================== END OF SCALING related data ================
C
C  Local variables
C
C     Interval associated to the subblocks of RHS a node has to process
      INTEGER, DIMENSION(:), ALLOCATABLE, TARGET :: RHS_BOUNDS
      INTEGER                            :: LPTR_RHS_BOUNDS
      INTEGER, DIMENSION(:), POINTER     :: PTR_RHS_BOUNDS
      LOGICAL  :: DO_NBSPARSE, NBSPARSE_LOC
      LOGICAL  :: PRINT_MAXAVG
      DOUBLE PRECISION ARRET
      COMPLEX(kind=8) C_DUMMY(1)
      DOUBLE PRECISION R_DUMMY(1)
      INTEGER IDUMMY(1), JDUMMY(1), KDUMMY(1), LDUMMY(1), MDUMMY(1)
      INTEGER, TARGET :: IDUMMY_TARGET(1)
      COMPLEX(kind=8), TARGET :: CDUMMY_TARGET(1)
      INTEGER JJ
      INTEGER allocok
      INTEGER NBRHS, NBRHS_EFF, BEG_RHS, NB_RHSSKIPPED,
     &        LD_RHS,
     &        MASTER_ROOT, MASTER_ROOT_IN_COMM
      INTEGER SIZE_ROOT, LD_REDRHS
      INTEGER(8) :: IPT_RHS_ROOT
      INTEGER(8) :: IBEG, IBEG_RHSCOMP, KDEC, IBEG_loc, IBEG_REDRHS
      INTEGER LD_RHSCOMP, NCOL_RHS_loc
      INTEGER LD_RHS_loc, JBEG_RHS_loc
      INTEGER NB_K133, IRANK, TSIZE
      INTEGER KMAX_246_247
      INTEGER IFLAG_IR, IRStep
      LOGICAL TESTConv
      LOGICAL WORKSPACE_MINIMAL_PREFERRED, WK_USER_PROVIDED
      INTEGER(8) NB_BYTES       !size of data allocated during solve
      INTEGER(8) NB_BYTES_MAX   !MAX size of data allocated during solve
      INTEGER(8) NB_BYTES_EXTRA !For Step2Node, which may be freed later
      INTEGER(8) NB_BYTES_LOC   !For temp. computations
      INTEGER(8) NB_INT, NB_CMPLX, NB_CHAR, K34_8, K35_8
      INTEGER(8) K16_8, ITMP8, NB_BYTES_ON_ENTRY
#if defined(V_T)
C  Vampir
      INTEGER soln_drive_class, glob_comm_ini, perm_scal_ini, soln_dist,
     &        soln_assem, perm_scal_post
#endif
      LOGICAL I_AM_SLAVE, BUILD_POSINRHSCOMP
      LOGICAL :: BUILD_RHSMAPINFO
      LOGICAL WORK_WCB_ALLOCATED, IS_INIT_OOC_DONE
      LOGICAL :: IS_LR_MOD_TO_STRUC_DONE
      INTEGER :: KEEP350_SAVE
      LOGICAL STOP_AT_NEXT_EMPTY_COL
      INTEGER  MTYPE_LOC
      INTEGER(4) :: I4
      INTEGER  MAT_ALLOC_LOC, MAT_ALLOC
      INTEGER MUMPS_PROCNODE
      EXTERNAL mumps_procnode
      INTEGER(8) :: FILE_SIZE,STRUC_SIZE
C
C  First executable statement
C
#if defined(V_T)
      CALL vtclassdef( 'Soln driver',soln_drive_class,ierr)
      CALL vtfuncdef( 'glob_comm_ini',soln_drive_class,
     &     glob_comm_ini,ierr)
      CALL vtfuncdef( 'perm_scal_ini',soln_drive_class,
     &     perm_scal_ini,ierr)
      CALL vtfuncdef( 'soln_dist',soln_drive_class,soln_dist,ierr)
      CALL vtfuncdef( 'soln_assem',soln_drive_class,soln_assem,ierr)
      CALL vtfuncdef( 'perm_scal_post',soln_drive_class,
     &     perm_scal_post,ierr)
#endif
C     -- The following pointers xxCOPY might be allocated but then
C     -- the associated xxCOPY_ALLOCATED will be set to
C     -- enable deallocation
      irhs_ptr_copy => idummy_target
      irhs_ptr_copy_allocated = .false.
      irhs_sparse_copy => idummy_target
      irhs_sparse_copy_allocated=.false.
      rhs_sparse_copy => cdummy_target
      rhs_sparse_copy_allocated=.false.
      NULLIFY(rhs_ir)
      NULLIFY(work_wcb)
      NULLIFY(scaling_data_dr%SCALING)
      NULLIFY(scaling_data_dr%SCALING_LOC)
      NULLIFY(scaling_data_sol%SCALING)
      NULLIFY(scaling_data_sol%SCALING_LOC)
      irhs_loc_ptr_allocated = .false.
      is_init_ooc_done   = .false.
      is_lr_mod_to_struc_done = .false.
      wk_user_provided   = .false.
      work_wcb_allocated = .false.
      cntl =>id%CNTL
      keep =>id%KEEP
      keep8=>id%KEEP8
      is   =>id%IS
      icntl=>id%ICNTL
      info =>id%INFO
C      ASPK =>id%A
C      COLSCA =>id%COLSCA
C      ROWSCA =>id%ROWSCA
      rinfog =>id%RINFOG
      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))
      print_maxavg = .NOT.(id%NSLAVES.EQ.1 .AND. keep(46).EQ.1)
      IF (.not.prok)  mp =0
      IF (.not.prokg) mpg=0
      IF ( prok  ) WRITE(mp,100)
      IF ( prokg ) WRITE(mpg,100)
      nb_bytes       = 0_8
      nb_bytes_max   = 0_8
      nb_bytes_extra = 0_8
      k34_8    = int(keep(34), 8)
      k35_8    = int(keep(35), 8)
      k16_8    = int(keep(16), 8)
      nbent_rhscomp = 0
C     Used by DISTRIBUTED_SOLUTION to skip empty columns
C     that are skipped (case of sparse RHS)
      nb_rhsskipped = 0
C     next 4 initialisations needed in case of error
C     to free space allocated
      lscal              = .false.
      work_wcb_allocated = .false.
      icntl21  = -99998  ! will be bcasted later to slaves
      ibeg_rhscomp =-152525_8  ! Should not be used
      build_posinrhscomp = .true.
      ibeg_glob_def = -9888  ! unitialized state
      iend_glob_def = -9888  ! unitialized state
      ibeg_root_def = -9777  ! unitialized state
      iend_root_def = -9777  ! unitialized state
      iroot_def_rhs_col1 = -9666 ! unitialized state
C     Not needed anymore (since new version of gather)
C     LD_RHSCOMP = max(KEEP(89),1)  ! at the nb of pivots eliminated on
                                    ! that proc
      ld_rhscomp = 1
      nb_fs_rhscomp_tot = keep(89)
!     number of FS var of the pruned tree
!     mapped on this proc
      nb_fs_rhscomp_f = nb_fs_rhscomp_tot
C     Save value of KEEP(350), in case of LR solve
C     KEEP(350) may be overwritten and restored
C     Old unoptimized version before 5.0.2 not available anymore
      IF (keep(350).LE.0) keep(350)=1
      IF (keep(350).GT.2) keep(350)=1
      keep350_save = keep(350)
C
C     Depending on the type of parallelism,
C     the master can have the role of a slave
      i_am_slave = ( id%MYID .ne. master  .OR.
     &             ( id%MYID .eq. master .AND.
     &               keep(46) .eq. 1 ) )
C     
C     Compute the number of integers and nb of reals in the structure
      CALL zmumps_size_in_struct (id, nb_int, nb_cmplx, nb_char)
      nb_bytes = nb_bytes + nb_int * k34_8 + nb_cmplx * k35_8 + nb_char
      nb_bytes_on_entry = nb_bytes !used to check alloc/dealloc count ok
      CALL zmumps_compute_memory_save(id,file_size,struc_size)
      nb_bytes_max = max(nb_bytes_max,nb_bytes)
C     ======================================
C     BEGIN CHECK KEEP ENTRIES AND INTERFACE
C     ======================================
C     The checks below used to be in ZMUMPS_DRIVER. It is much better
C     to have them here in ZMUMPS_SOL_DRIVER because this enables
C     more flexibility in the management of priorities between various
C     checks.
      IF (id%MYID .EQ. master) THEN
c     subroutine only because called at facto and solve
          CALL zmumps_set_k221(id)
          id%KEEP(111) = id%ICNTL(25)
C         For the case of ICNTL(20)=1 one could
C         switch off exploit sparsity when RHS is too dense.
          IF (id%ICNTL(20) .EQ. 1) id%KEEP(235) = -1 !automatic
          IF (id%ICNTL(20) .EQ. 2) id%KEEP(235) = 0  !off
          IF (id%ICNTL(20) .EQ. 3) id%KEEP(235) = 1  !on
          IF (id%ICNTL(20).EQ.1 .or. id%ICNTL(20).EQ.2 .or.
     &       id%ICNTL(20).EQ.3) THEN
            id%KEEP(248) = 1 !sparse RHS
          ELSE IF (id%ICNTL(20).EQ.10 .OR. id%ICNTL(20).EQ.11) THEN
            id%KEEP(248) = -1 ! dist. RHS
          ELSE
            id%KEEP(248) = 0 !dense RHS
          ENDIF
          icntl21      = id%ICNTL(21)
          IF (icntl21 .ne.0.and.icntl21.ne.1) icntl21=0
          IF ( id%ICNTL(30) .NE.0 ) THEN
C           A-1 is on
            id%KEEP(237) = 1
          ELSE
C           A-1 is off
            id%KEEP(237) = 0
          ENDIF
          IF (id%KEEP(248) .eq.0.and. id%KEEP(237).ne.0) THEN
C            For A-1 we have a sparse RHS in the API.
C            Force KEEP(248) accordingly.
             id%KEEP(248)=1
          ENDIF
          IF ((id%KEEP(221).EQ.2 ).AND.(id%KEEP(248).NE.0) ) THEN
C          -- input RHS is indeed stored in REDRHS and RHSCOMP
           id%KEEP(248) = 0
          ENDIF
          IF ((id%KEEP(221).EQ.2 ).AND.(id%KEEP(235).NE.0) ) THEN
C          -- input RHS is in fact effectively
C          -- stored in REDRHS and RHSCOMP
           id%KEEP(235) = 0
          ENDIF
          IF ( (id%KEEP(248).EQ.0).AND.(id%KEEP(111).EQ.0) ) THEN
C           RHS is not sparse and thus exploit sparsity is reset to 0
            id%KEEP(235) = 0
          ENDIF
          IF (keep(248) .EQ. -1) THEN
C           V0 distributed RHS: no ES
            id%KEEP(235) = 0
          ENDIF
C         Case of Automatic setting of exploit sparsity (KEEP(235)=-1)
C         (in MUMPS_DRIVER original value of KEEP(235) is reset)
          IF(id%KEEP(111).NE.0) id%KEEP(235)=0
C
          IF (id%KEEP(235).EQ.-1) THEN
            IF (id%KEEP(237).NE.0) THEN
C            for A-1
              id%KEEP(235)=1
            ELSE
              id%KEEP(235)=1
            ENDIF
          ELSE IF (id%KEEP(235).NE.0) THEN
            id%KEEP(235)=1
          ENDIF
C         Setting of KEEP(242) (permute RHS)
          IF ((keep(111).NE.0)) THEN
C          In the context of null space, the null pivots
C          are by default permuted to post-order
C          However for null space there is in this case no need to
C          permute null pivots since they are already in correct order.
C          Setting KEEP(242)=1 would just force to go through
C          part of the code permuting to identity.
C          Apart for validation purposes this is not interesting
C          costly (and more risky).
              keep(242)   = 0
          ENDIF
          IF (keep(248).EQ.0.AND.keep(111).EQ.0) THEN
C             Permutation possible if sparse RHS
C             (KEEP(248).NE.0: A-1 or General Sparse)
C             or null space (even if in current version 
C                            it is deactived)
              keep(242)   = 0
          ENDIF
          IF ((keep(242).NE.0).AND.keep(237).EQ.0) THEN
            IF ((keep(242).NE.-9).AND.keep(242).NE.1.AND.
     &           keep(242).NE.-1) THEN
C            Reset it to 0 
             keep(242) = 0
            ENDIF
          ENDIF
          IF (keep(242).EQ.-9) THEN
C           {
C           Automatic setting of permute RHS
            IF (id%KEEP(237).NE.0) THEN
              keep(242) = 1  ! postorder for A-1
            ELSE ! dense or general sparse or distributed RHS
              keep(242) = 0  ! no permutation in most general case
              IF (keep(248) .EQ. 1) THEN ! sparse RHS
                IF (id%KEEP(235) .EQ. 1) THEN ! Tree pruning
                  IF (id%NRHS .GT. 1) THEN
                    IF (keep(497).EQ.-1 .OR. keep(497).GE.1) THEN
                      keep(242)=1
                    ENDIF
                  ENDIF
                ENDIF
              ENDIF
            ENDIF
C           }
          ENDIF
          IF ( (id%KEEP(221).EQ.1 ).AND.(id%KEEP(235).NE.0) ) THEN
C          -- Do not permute RHS with REDRHS for the time being
            id%KEEP(242) = 0
          ENDIF
          IF (keep(242).EQ.0) keep(243)=0 ! interleave off
          IF ((keep(237).EQ.0).OR.(keep(242).EQ.0)) THEN
C             Interleave (243) possible only
C             when permute RHS (242) is on and with A-1
              keep(243) = 0
          ENDIF
          IF (id%KEEP(237).EQ.1) THEN  ! A-1 entries
C           Case of automatic setting of KEEP(243), KEEP(493-498)
C           (exploit sparsity parameters)
           IF (id%NSLAVES.EQ.1) THEN
            IF (id%KEEP(243).EQ.-1) id%KEEP(243)=0  
            IF (id%KEEP(495).EQ.-1) id%KEEP(495)=1
            IF (id%KEEP(497).EQ.-1) id%KEEP(497)=1
           ELSE
            IF (id%KEEP(243).EQ.-1) id%KEEP(243)=1
            IF (id%KEEP(495).EQ.-1) id%KEEP(495)=1
            IF (id%KEEP(497).EQ.-1) id%KEEP(497)=1
           ENDIF
          ELSE ! dense or general sparse or distributed RHS
            id%KEEP(243)=0
            id%KEEP(495)=0
            IF (keep(248) .EQ. 1) THEN ! sparse RHS
              IF (id%KEEP(235) .EQ. 1) THEN ! Tree pruning
                IF (id%NRHS .GT. 1) THEN
                  IF (id%KEEP(497).EQ.-1) id%KEEP(497)=1
                ENDIF
              ENDIF
            ELSE
C             nbsparse meaningless for distributed or dense RHS 
C             Force it to 0 whatever was the initial value
              id%KEEP(497)=0
            ENDIF
          ENDIF
          mtype = id%ICNTL(  9 )
          IF (mtype.NE.1) mtype=0  ! see interface
          IF ((mtype.EQ.0).AND.keep(50).NE.0) mtype =1
!         suppress option Atx=b for A-1
          IF (id%KEEP(237).NE.0) mtype = 1
C
C         ICNTL(35) was defined at analysis and 
C         consistently reset at factorization
C         It was stored in KEEP(486) after factorization 
C         Set KEEP(485) accordingly.
C
          IF (keep(486) .EQ. 2) THEN 
            keep(485) = 1 ! BLR solve
          ELSE
            keep(485) = 0 ! FR solve
          ENDIF
      ENDIF
C     Bcast id%KEEP(401) strategy (which
C     may be switched off or on during solve)
      CALL mpi_bcast( id%KEEP(401), 1, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast(mtype,1,mpi_integer,master,
     &               id%COMM,ierr)
      CALL mpi_bcast( id%KEEP(111), 1, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( id%KEEP(221), 1, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( id%KEEP(235), 1, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( id%KEEP(237), 1, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( id%KEEP(242), 2, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( id%KEEP(248), 1, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( id%KEEP(350), 1, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( id%KEEP(485), 1, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( id%KEEP(495), 3, mpi_integer, master, id%COMM,
     &                  ierr )
      CALL mpi_bcast( icntl21, 1, mpi_integer, master, id%COMM, ierr )
C     Broadcast original id%NRHS (used at least for checks on SOL_loc
C                       and to allocate PERM_RHS in case of exploit sparsity)
      CALL mpi_bcast( id%NRHS,1, mpi_integer, master, id%COMM,ierr)
C
C     TIMINGS: reset to 0
      timec2=0.0d0
      timecopyscale2=0.0d0
      timegather2=0.0d0
      timescatter2=0.0d0
      id%DKEEP(112)=0.0d0
      id%DKEEP(113)=0.0d0
C     id%DKEEP(114) time for the iterative refinement
C     id%DKEEP(120) time for the error analysis
C     id%DKEEP(121) time for condition number
C     id%DKEEP(122) time for matrix redistribution (copy+scale solution)
      id%DKEEP(114)=0.0d0
      id%DKEEP(120)=0.0d0
      id%DKEEP(121)=0.0d0
      id%DKEEP(115)=0.0d0
      id%DKEEP(116)=0.0d0
      id%DKEEP(122)=0.0d0
C     Time for fwd, bwd and scalapack is
C     accumulated in DKEEP(117-119) within SOL_C
C     If requested time for each call to FWD/BWD
C     might be print but on output to solve
C     phase DKEEP will hold on each proc the accumulated time
      id%DKEEP(117)=0.0d0
      id%DKEEP(118)=0.0d0
      id%DKEEP(119)=0.0d0
      id%DKEEP(123)=0.0d0 
      id%DKEEP(124)=0.0d0 
      id%DKEEP(125)=0.0d0 
      id%DKEEP(126)=0.0d0 
      id%DKEEP(127)=0.0d0 
      id%DKEEP(128:134)=0.0d0
      id%DKEEP(140:153)=0.0d0
C
      CALL mumps_secdeb(time3)
C     ------------------------------
C     Check parameters on the master
C     ------------------------------
      IF ( id%MYID .EQ. master ) THEN
         IF ((keep(23).NE.0).AND.keep(50).NE.0) THEN
C          Maximum transversal permutation
C          has not been saved (KEEP(23)>0 and UNS_PERM allocated)
C          when matrix is symmetric.
           IF (prokg) WRITE(mpg,'(A)')
     &       ' Internal Error 1 in solution driver '
           id%INFO(1)=-444
           id%INFO(2)=keep(23)
          ENDIF
C         ------------------------------------
C         Check that factors are available
C         either in-core or on disk, case
C         where factors were discarded during
C         factorization (e.g. useful to simulate
C         an OOC factorization or just get nb of
C         negative pivots or determinant)
C         ------------------------------------
          IF (keep(201) .EQ. -1) THEN
             IF (prokg) THEN
                WRITE(mpg,'(A)')
     &               ' ERROR: Solve impossible because factors not kept'
             ENDIF
             id%INFO(1)=-44
             id%INFO(2)=keep(251)
             GOTO 333
          ELSE IF (keep(221).EQ.0 .AND. keep(251) .EQ. 2
     &            .AND. keep(252).EQ.0) THEN
             IF (prokg) THEN
                WRITE(mpg,'(A)')
     &               ' ERROR: Solve impossible because factors not kept'
             ENDIF
             id%INFO(1)=-44
             id%INFO(2)=keep(251)
             GOTO 333
          ENDIF
C         ------------------
          IF (keep(252).NE.0 .AND. id%NRHS .NE. id%KEEP(253)) THEN
C            Fwd in facto
C            KEEP(252-253) available on all procs since analysis phase
C            Error: id%NRHS is not allowed to change since analysis
C            because fwd has been performed during facto with
C            KEEP(253) RHS
             IF (prokg) THEN
                WRITE(mpg,'(A)')
     &               ' ERROR: id%NRHS not allowed to change when',
     &               ' ICNTL(32)=1'
             ENDIF
             id%INFO(1)=-42
             id%INFO(2)=id%KEEP(253)
             GOTO 333
          ENDIF
C     Testing MTYPE instead of ICNTL(9)
          IF (keep(252).NE.0 .AND. mtype.NE.1) THEN
C            Fwd in facto is not compatible with transpose system
             info(1) = -43
             info(2) = 9
             IF (prokg) THEN
                WRITE(mpg,'(A)')
     &               .NE.' ERROR: Transpose system (ICNTL(9)0) not ',
     &               ' compatible with forward performed during',
     &               ' factorization (ICNTL(32)=1)'
             ENDIF
             GOTO 333
          ENDIF
          IF (keep(248) .NE. 0.AND.keep(252).NE.0) THEN
C            Fwd during facto incompatible with sparse RHS
C            Forbid sparse RHS when Fwd performed during facto
C            Sparse RHS may be due to A-1 (ICNTL(30)
             info(1) = -43
             IF (keep(237).NE.0) THEN
               info(2) = 30 ! ICNTL(30)
               IF (prokg) THEN
                  WRITE(mpg,'(A)')
     &                 ' ERROR: A-1 functionality incompatible with',
     &                 ' forward performed during factorization',
     &                 ' (ICNTL(32)=1)'
               ENDIF
             ELSE
               info(2) = 20      ! ICNTL(20)
               IF (prokg) THEN
                  WRITE(mpg,'(A)')
     &         ' ERROR: sparse or dist. RHS incompatible with forward',
     &         ' elimination during factorization (ICNTL(32)=1)'
               ENDIF
             ENDIF
             GOTO 333
          ENDIF
          IF (keep(237) .NE. 0 .AND. icntl21.NE.0) THEN
             IF (prokg) THEN
                WRITE(mpg,'(A)')
     &               ' ERROR: A-1 functionality is incompatible',
     &               ' with distributed solution.'
             ENDIF
             info(1)=-48
             info(2)=21
             GOTO 333
          ENDIF
          IF (keep(237) .NE. 0 .AND. keep(60) .NE.0) THEN
             IF (prokg)  THEN
                WRITE(mpg,'(A)')
     &               ' ERROR: A-1 functionality is incompatible',
     &               ' with Schur.'
             ENDIF
             info(1)=-48
             info(2)=19
             GOTO 333
          ENDIF
          IF (keep(237) .NE. 0 .AND. keep(111) .NE.0) THEN
             IF (prokg) THEN
                WRITE(mpg,'(A)')
     &               ' ERROR: A-1 functionality is incompatible',
     &               ' with null space.'
             ENDIF
             info(1)=-48
             info(2)=25
             GOTO 333
          ENDIF
          IF (id%NRHS .LE. 0) THEN
             id%INFO(1)=-45
             id%INFO(2)=id%NRHS
             IF ((id%KEEP(111).NE.0).AND.(id%INFOG(28).EQ.0)) THEN
               IF (prokg) THEN
                WRITE(mpg,'(A)')
     &              'ICNTL(25) NE 0 but INFOG(28)=0',
     &              '  the matrix is not deficient'
               ENDIF
             ENDIF
             GOTO 333
          ENDIF
C         Entries of A-1 are stored in place of the input sparse RHS
C         thus no need for RHS to be allocated.
          IF ( (id%KEEP(237).EQ.0) ) THEN
             IF ((id%KEEP(248) == 0 .AND.keep(221).NE.2)
     &            .OR. icntl21==0) THEN
C              RHS must be of size N on the master either to
C              store the dense centralized RHS, either to store
C              the dense centralized solution.
               CALL zmumps_check_dense_rhs
     &         (id%RHS,id%INFO,id%N,id%NRHS,id%LRHS)
               IF (id%INFO(1) .LT. 0) GOTO 333
             ENDIF
          ELSE
C           Check that the constraint NRHS=N is respected
C           Check for valid sparse RHS structure done
            IF (id%NRHS .NE. id%N) THEN
              id%INFO(1)=-47
              id%INFO(2)=id%NRHS
              GOTO 333
            ENDIF
          ENDIF
          IF (id%KEEP(248) == 1) THEN
C           ------------------------------------
C           RHS_SPARSE, IRHS_SPARSE and IRHS_PTR
C           must be allocated of adequate size
C           ------------------------------------
            IF (( id%NZ_RHS .LE.0 ).AND.(keep(237).NE.0)) THEN
C             At least one entry of A-1 must be requested
              id%INFO(1)=-46
              id%INFO(2)=id%NZ_RHS
              GOTO 333
            ENDIF
            IF (( id%NZ_RHS .LE.0 ).AND.(keep(221).EQ.1)) THEN
C             At least one entry of RHS must be nonzero with
c             Schur reduced RHS option
              id%INFO(1)=-46
              id%INFO(2)=id%NZ_RHS
              GOTO 333
            ENDIF
            IF ( id%NZ_RHS .GT. 0 ) THEN
              IF ( .not. associated(id%RHS_SPARSE) )THEN
                id%INFO(1)=-22
                id%INFO(2)=10
                GOTO 333
              ENDIF
            ENDIF
            IF (id%NZ_RHS .GT. 0) THEN
              IF ( .not. associated(id%IRHS_SPARSE) )THEN
                id%INFO(1)=-22
                id%INFO(2)=11
                GOTO 333
              ENDIF
            ENDIF
            IF ( .not. associated(id%IRHS_PTR) )THEN
              id%INFO(1)=-22
              id%INFO(2)=12
              GOTO 333
            ENDIF
C
            IF (size(id%IRHS_PTR) < id%NRHS + 1) THEN
              id%INFO(1)=-22
              id%INFO(2)=12
              GOTO 333
            END IF
            IF (id%IRHS_PTR(id%NRHS + 1).ne.id%NZ_RHS+1) THEN
              id%INFO(1)=-27
              id%INFO(2)=id%IRHS_PTR(id%NRHS+1)
              GOTO 333
            END IF
C           compare with dble to prevent overflow
            IF (dble(id%N)*dble(id%NRHS).LT.dble(id%NZ_RHS)) THEN
C             Possible in case of dupplicate entries in Sparse RHS
              IF (prokg) THEN
                 write(mpg,*)
     &              " WARNING: many dupplicate entries in ", 
     &              " sparse RHS  provided by the user ",
     &              " id%NZ_RHS,id%N,id%NRHS =",
     &              id%NZ_RHS,id%N,id%NRHS
              ENDIF
            END IF
            IF (id%IRHS_PTR(1).ne.1) THEN
              id%INFO(1)=-28
              id%INFO(2)=id%IRHS_PTR(1)
              GOTO 333
            END IF
            IF (size(id%IRHS_SPARSE) < id%NZ_RHS) THEN
              id%INFO(1)=-22
              id%INFO(2)=11
              GOTO 333
            END IF
            IF (size(id%RHS_SPARSE) < id%NZ_RHS) THEN
              id%INFO(1)=-22
              id%INFO(2)=10
              GOTO 333
            END IF
          ENDIF
C         --------------------------------
C         Set null space options for solve
C         --------------------------------
          CALL zmumps_get_ns_options_solve(icntl(1),keep(1),
     &                                     id%NRHS,
     &                                     mpg,info(1))
          IF (info(1) .LT. 0) GOTO 333
C
      END IF                   ! MASTER
C     --------------------------------------
C     Check distributed solution vectors
C     --------------------------------------
      IF (icntl21==1) THEN
          IF ( i_am_slave ) THEN
C           (I)SOL_loc should be allocated to hold the
C           distributed solution on exit
            IF ( id%LSOL_loc < id%KEEP(89) ) THEN
              id%INFO(1)= -29
              id%INFO(2)= id%LSOL_loc
              GOTO 333
            ENDIF
            IF (id%KEEP(89) .NE. 0) THEN
              IF ( .not. associated(id%ISOL_loc) )THEN
                id%INFO(1)=-22
                id%INFO(2)=13
                GOTO 333
              ENDIF
              IF ( .not. associated(id%SOL_loc) )THEN
                id%INFO(1)=-22
                id%INFO(2)=14
                GOTO 333
              ENDIF
              IF (size(id%ISOL_loc) < id%KEEP(89) ) THEN
                id%INFO(1)=-22
                id%INFO(2)=13
                GOTO 333
              END IF
#             if defined(MUMPS_F2003)
              IF (size(id%SOL_loc,kind=8) <
     &              int(id%NRHS-1,8)*int(id%LSOL_loc,8)+
     &              int(id%KEEP(89),8)) THEN
                id%INFO(1)=-22
                id%INFO(2)=14
                GOTO 333
              END IF
#             else
C             Warning: size returns a standard INTEGER and could
C             overflow if id%SOL_loc was allocated of size > 2^31-1;
C             still we prefer to perform this test since only (1) very
C             large problems with large NRHS and small numbers of MPI
C             can result in such a situation; (2) the test could be
C             suppressed if needed but might be still be ok in case
C             the right-hand side overflows too.
              IF (size(id%SOL_loc) <
     &              (id%NRHS-1)*id%LSOL_loc+id%KEEP(89)) THEN
                id%INFO(1)=-22
                id%INFO(2)=14
                GOTO 333
              END IF
#             endif
            ENDIF
          ENDIF
      ENDIF
      IF (id%MYID .NE. master) THEN
          IF (id%KEEP(248) == 1) THEN
C          RHS should NOT be associated
C          if I am not master since it is
C          not even used to store the solution
           IF ( associated( id%RHS ) ) THEN
             id%INFO( 1 ) = -22
             id%INFO( 2 ) = 7
             GOTO 333
           END IF
           IF ( associated( id%RHS_SPARSE ) ) THEN
             id%INFO( 1 ) = -22
             id%INFO( 2 ) = 10
             GOTO 333
           END IF
           IF ( associated( id%IRHS_SPARSE ) ) THEN
             id%INFO( 1 ) = -22
             id%INFO( 2 ) = 11
             GOTO 333
           END IF
           IF ( associated( id%IRHS_PTR ) ) THEN
             id%INFO( 1 ) = -22
             id%INFO( 2 ) = 12
             GOTO 333
           END IF
          END IF
      ENDIF
      IF (i_am_slave .AND. id%KEEP(248).EQ.-1) THEN
        CALL zmumps_check_distrhs(
     &       id%Nloc_RHS,
     &       id%LRHS_loc,
     &       id%NRHS,
     &       id%IRHS_loc,
     &       id%RHS_loc,
     &       id%INFO)
        IF (id%INFO(1) .LT. 0) GOTO 333
      ENDIF
C     Prepare pointers to pass POINTERS(1) to
C     routines with implicit interfaces which
C     will then assume contiguous information
C     without needing to copy pointer arrays
C     in and out. Do this even if KEEP(248)
C     is different from -1 because of the
C     call to ZMUMPS_DISTSOL_INDICES
      IF (associated(id%IRHS_loc)) THEN
        IF (size(id%IRHS_loc) .NE. 0) THEN
          irhs_loc_ptr=>id%IRHS_loc
        ELSE
C         so that IRHS_loc_PTR(1) is ok
          irhs_loc_ptr=>idummy_target
        ENDIF
      ELSE
        irhs_loc_ptr=>idummy_target
      ENDIF
      IF (associated(id%RHS_loc)) THEN
        IF (size(id%RHS_loc) .NE. 0) THEN
          idrhs_loc=>id%RHS_loc
        ELSE
          idrhs_loc=>cdummy_target
        ENDIF
      ELSE
        idrhs_loc=>cdummy_target
      ENDIF
      IF (i_am_slave .AND. icntl21.EQ.1 .AND.
     &    keep(248) .EQ. -1) THEN  ! Dist RHS and dist solution
        IF (associated(id%RHS_loc) .AND.
     &      associated(id%SOL_loc)) THEN
          IF (id%KEEP(89).GT.0) THEN
C           ----------------------------------------------------
C           Check if RHS_loc and SOL_loc point to same object...
C           id%SOL_loc(1) ok otherwise an error -22/14
C           would have been raised earlier.
C           idRHS_loc(1) may point to CDUMMY but is ok
C           ----------------------------------------------------
            CALL mumps_size_c(idrhs_loc(1),id%SOL_loc(1),
     &           diff_sol_loc_rhs_loc)
C           ----------------------------------------
C           Check for compatible dimensions in case
C           SOL_loc and RHS_loc point to same memory
C           ----------------------------------------
            IF (diff_sol_loc_rhs_loc .EQ. 0_8 .AND.
     &          id%LSOL_loc .GT. id%LRHS_loc) THEN
C              Note that, depending on the block size,
C              if all columns are processed in one
C              shot, this could still work. However,
C              and since this was forbidden in the UG,
C              we raise the error systematically
               id%INFO(1)=-56
               id%INFO(2)=id%LRHS_loc
               IF (lpok) THEN
                WRITE(lp,'(A,I9,A,I9)')
     &" ** Error RHS_loc and SOL_loc pointers match but LRHS_loc="
     &,id%LRHS_loc, " and LSOL_loc=", id%LSOL_loc
               ENDIF
            ENDIF
          ENDIF
        ENDIF
      ENDIF
      IF (id%MYID.EQ.master) THEN
C         Do some checks (REDRHS), depending on KEEP(221)
          CALL zmumps_check_redrhs(id)
      END IF ! MYID.EQ.MASTER
      IF (id%INFO(1) .LT. 0) GOTO 333
C     -------------------------
C     Propagate possible errors
C     -------------------------
 333  CONTINUE
      CALL mumps_propinfo( id%ICNTL(1),
     &                      id%INFO(1),
     &                      id%COMM, id%MYID )
      IF ( id%INFO(1) .LT. 0 ) GO TO 90
C     ====================================
C     END CHECK INTERFACE AND KEEP ENTRIES
C     ====================================
C     ====================================
C     Process case of NZ_RHS = 0  with
C     sparse RHS and General Sparse (NOT A-1)
C     -----------------------------------
      IF ((id%KEEP(248).EQ.1).AND.(id%KEEP(237).EQ.0)) THEN
C
       CALL mpi_bcast(id%NZ_RHS,1,mpi_integer,master,
     &               id%COMM,ierr)
C
       IF (id%NZ_RHS.EQ.0) THEN
C       We reset solution to zero and we return
C       (first freeing working space at label 90)
        IF ((icntl21.EQ.1).AND.(i_am_slave)) THEN
C       ----------------------
C       SOL_loc reset to zero
C       ----------------------
C         ----------------------
C         Prepare ISOL_loc array
C         ----------------------
          liw_passed=max(1,keep(32))
C         Only called if more than 1 pivot
C         was eliminated by the processor.
C         Note that LSOL_loc >= KEEP(89)
          IF (keep(89) .GT. 0) THEN
            CALL zmumps_distsol_indices( mtype, id%ISOL_loc(1),
     &               id%PTLUST_S(1),
     &               id%KEEP(1),id%KEEP8(1),
     &               id%IS(1), liw_passed,id%MYID_NODES,
     &               id%N, id%STEP(1), id%PROCNODE_STEPS(1),
     &               id%NSLAVES, scaling_data_sol, lscal
C                    For checking only
     &               , .false., idummy(1), 1
     &               )
            DO j=1, id%NRHS
              DO i=1, keep(89)
                id%SOL_loc((j-1)*id%LSOL_loc + i) =zero
              ENDDO
            ENDDO
          ENDIF
        ENDIF
        IF (icntl21.NE.1) THEN  ! centralized solution
C       ----------------------------
C       RHS reset to zero on master
C       ----------------------------
         IF (id%MYID.EQ.master) THEN
            DO j=1, id%NRHS
              DO i=1, id%N
                id%RHS(int(j-1,8)*int(id%LRHS,8) + int(i,8)) =zero
              ENDDO
            ENDDO
         ENDIF
        ENDIF
C
C       print solve phase stats if requested
        IF ( prokg )  THEN
C          write(6,*) " NZ_RHS is zero "
           WRITE( mpg, 150 )
C             ICNTL(35) should not been accessed during SOLVE thus
C             print KEEP(486) value set during factorization
     &        id%NRHS, icntl(27), icntl(9), icntl(10), icntl(11),
     &        icntl(20), icntl(21), icntl(30), keep(486)
           IF (keep(221).NE.0) THEN
            WRITE (mpg, 152) keep(221)
           ENDIF
           IF (keep(252).GT.0) THEN   ! Fwd during facto
            WRITE (mpg, 153) keep(252)
           ENDIF
        ENDIF
C
C       --------
        GOTO 90 ! end of solve deallocate what is needed
C       ====================================
C       END CHECK INTERFACE AND KEEP ENTRIES
C       ====================================
       ENDIF  ! test NZ_RHS.EQ.0
C      --------
      ENDIF ! (id%KEEP(248).EQ.1).AND.(id%KEEP(237).EQ.0)
      interleave_par   =.false.
      do_permute_rhs   =.false.
C
      IF ((id%KEEP(235).NE.0).or.(id%KEEP(237).NE.0)) THEN
C        Case of pruned elimination tree or selected entries in A-1
         IF (id%KEEP(237).NE.0.AND.
     &        id%KEEP(248).EQ.0) THEN
C         When A-1 is requested (keep(237).ne.0)
C         sparse RHS has been forced to be on.
          IF (lpok) THEN
           WRITE(lp,'(A,I4,I4)')
     &     ' Internal Error 2 in solution driver (A-1) ',
     &       id%KEEP(237), id%KEEP(248)
          ENDIF
          CALL mumps_abort()
         ENDIF
C        NBT is inout in MUMPS_REALLOC and should be initialized.
         nbt = 0
C        -- Allocate Step2node on each proc
         CALL mumps_realloc(id%Step2node, id%KEEP(28), id%INFO, lp,
     &        force=.true.,
     &        string='id%Step2node (Solve)', memcnt=nbt, errcode=-13)
         CALL mumps_propinfo( icntl(1), info(1),
     &        id%COMM, id%MYID )
         IF ( info(1).LT.0 ) RETURN
C        -- build Step2node on each proc;
C        -- this is usefull to have at each step a unique
C        -- representative node (associated with principal variable of
C        -- that node.
         IF (nbt.NE.0) THEN
          ! Step2node was reallocated and needs be recomputed
          DO i=1, id%N
           IF (id%STEP(i).LE.0) cycle  ! nonprincipal variables
           id%Step2node(id%STEP(i)) = i
          ENDDO
C        ELSE
C          we reuse Step2node computed in a previous solve phase
C          Step2node is deallocated each time a new analysis is
C          performed or when job=-2 is called
         ENDIF
         nb_bytes = nb_bytes + nbt*k34_8
         nb_bytes_max = max(nb_bytes_max,nb_bytes)
         nb_bytes_extra = nb_bytes_extra + nbt * k34_8
C        Mapping information used during solve. In case of several
C        facto+solve it has to be recomputed.
C        In case of several solves with the same
C        facto, it is not recomputed.
C        It used to compute the interleaving
C        for A-1, and, in dev_version, passed to sol_c to compute
C        some stats
         IF((keep(235).NE.0).OR.(keep(237).NE.0)) THEN
           IF(.NOT.associated(id%IPTR_WORKING)) THEN
             CALL zmumps_build_mapping_info(id)
           END IF
         END IF
      ENDIF
C
C     Initialize SIZE_OF_BLOCK from MUMPS_SOL_ES module
      IF ( i_am_slave )
     &  CALL zmumps_sol_es_init(id%OOC_SIZE_OF_BLOCK, id%KEEP(201))
      do_null_piv = .true.
      nbcol_inbloc = -9998
      nz_this_block= -9998
      jbeg_rhs  = -9998
c
      IF (id%MYID.EQ.master) THEN ! Compute NRHS_NONEMPTY
C
C       -- Sparse RHS does
        IF ( keep(111)==0 .AND. keep(248)==1
     &   ) THEN
C       -- Note that KEEP(111).NE.0 (null space on)
C       -- and KEEP(248).NE.0 will be made incompatible
C       -- When computing entries of A-1 (or SparseRHS only)
           nrhs_nonempty = 0
           DO i=1, id%NRHS
              IF (id%IRHS_PTR(i).LT.id%IRHS_PTR(i+1))
     &             nrhs_nonempty = nrhs_nonempty+1 !ith col in non empty
           ENDDO
           IF (nrhs_nonempty.LE.0) THEN
C           Internal error: tested before in mumps_driver
            IF (lpok)
     &        WRITE(lp,*) " Internal Error 3 in solution driver ",
     &                    " NRHS_NONEMPTY= ",
     &        nrhs_nonempty
            CALL mumps_abort()
           ENDIF
        ELSE
           nrhs_nonempty = id%NRHS
        ENDIF
      ENDIF
C     ------------------------------------
C     If there is a special root node,
C     precompute mapping of root's master
C     ------------------------------------
      size_root   = -33333
      IF ( keep( 38 ) .ne. 0 ) THEN
            master_root = mumps_procnode(
     &                    id%PROCNODE_STEPS(id%STEP( keep(38))),
     &                    keep(199) )
            IF (id%MYID_NODES .eq. master_root) THEN
              size_root = id%root%TOT_ROOT_SIZE
            ELSE IF ((id%MYID.EQ.master).AND.keep(60).NE.0) THEN
C             SIZE_ROOT also used for KEEP(221).NE.0
              size_root=id%KEEP(116)
            ENDIF
      ELSE IF (keep( 20 ) .ne. 0 ) THEN
            master_root = mumps_procnode(
     &                    id%PROCNODE_STEPS(id%STEP(keep(20))),
     &                    keep(199) )
            IF (id%MYID_NODES .eq. master_root) THEN
              size_root = id%IS(
     &               id%PTLUST_S(id%STEP(keep(20)))+keep(ixsz) + 3)
            ELSE IF ((id%MYID.EQ.master).AND.keep(60).NE.0) THEN
C             SIZE_ROOT also used for KEEP(221).NE.0
              size_root=id%KEEP(116)
            ENDIF
      ELSE
            master_root = -44444
      END IF
C     --------------
C     Get block size
C     --------------
C     We work on a maximum of NBRHS at a time.
C     The leading dimension of RHS is id%LRHS on the host process
C     and it is set to N on slave processes.
      IF (id%MYID .eq. master) THEN
        keep(84) = icntl(27)
C     Treating ICNTL(27)=0 as if ICNTL(27)=1
        IF(icntl(27).EQ.0) keep(84)=1
        IF (keep(252).NE.0) THEN
!       Fwd in facto: all rhs (KEEP(253) need be processed in one pass
          nbrhs = keep(253)
        ELSE
          IF (keep(201) .EQ. 0 .OR. keep(84) .GT. 0) THEN
            nbrhs = abs(keep(84))
          ELSE
            nbrhs = -2*keep(84)
          END IF
          IF (nbrhs .GT. nrhs_nonempty ) nbrhs = nrhs_nonempty
C
        ENDIF
      ENDIF
#if defined(V_T)
      CALL vtbegin(glob_comm_ini,ierr)
#endif
C     NRHS_NONEMPTY needed on all procs to allocate RHSCOMP on slaves
      CALL mpi_bcast(nrhs_nonempty,1,mpi_integer,master,
     &               id%COMM,ierr)
      CALL mpi_bcast(nbrhs,1,mpi_integer,master,
     &               id%COMM,ierr)
C
      IF (keep(201).GT.0) THEN
C         --- id%KEEP(201) indicates if OOC is on (=1) of not (=0)
C         -- 107: number of buffers
C         Define number of types of files (L, possibly U)
          workspace_minimal_preferred = .false.
          IF (id%MYID .eq. master) THEN
             keep(107) = max(0,keep(107))
             IF ((keep(107).EQ.0).AND.
     &            (keep(204).EQ.0).AND.(keep(211).NE.1) ) THEN
C             -- default setting for release 4.8
              ! Case of
              !  -Emmergency buffer only and
              !  -Synchronous mode
              !  -NO_O_DIRECT (because of synchronous choice)
              ! THEN
              !   "Basic system-based version"
              !   We can force to allocate S to a minimal
              !   value.
              workspace_minimal_preferred=.true.
             ENDIF
          ENDIF
          CALL mpi_bcast( keep(107), 1, mpi_integer,
     &                  master, id%COMM, ierr )
          CALL mpi_bcast( keep(204), 1, mpi_integer,
     &                  master, id%COMM, ierr )
          CALL mpi_bcast( keep(208), 2, mpi_integer,
     &                  master, id%COMM, ierr )
          CALL mpi_bcast( workspace_minimal_preferred, 1,
     &                  mpi_logical,
     &                  master, id%COMM, ierr )
C       --- end of OOC case
      ENDIF
      IF ( i_am_slave ) THEN
C
C       NB_K133:  Max number of simultaneously processed
C           active fronts.
C         Why more than one active node ?
C         1/ In parallel when we start a level 2 node
C           then we do not know exactly when we will
C           have received all contributions from the
C           slaves.
C           This is very critical in OOC since the
C           size provided to the solve phase is
C           much smaller and since we need
C           to determine the size fo the buffers for IO.
C            We pospone the allocation of the block NFRONT*NB_NRHS
C            and solve the problem.
C
C
C         2/ While processing a node and sending information
C            if we have not enough memory in send buffer
C            then we must receive.
C            We feel that this is not so critical.
C
        nb_k133     = 3
C
C       To this we must add one time KEEP(133) to store
C       the RHS of the root node if the root is local.
C       Furthermore this quantity has to be multiplied by the
C       blocking size in case of multiple RHS.
C
        IF ( keep( 38 ) .NE. 0 .OR. keep( 20 ) .NE. 0 ) THEN
          IF ( master_root .eq. id%MYID_NODES ) THEN
            IF (
     &          .NOT. associated(id%root%RHS_CNTR_MASTER_ROOT)
     &         ) THEN
                nb_k133 = nb_k133 + 1
            ENDIF
          END IF
        ENDIF
        lwcb8_min = int(nb_k133,8)*int(keep(133),8)*int(nbrhs,8)
C
C       ---------------------------------------------------------------
C       Set WK_USER_PROVIDED to true when workspace WK_USER is provided
C       by user
C       We can accept WK_USER to be provided on only one proc and
C       different values of WK_USER per processor. Note that we are
C       inside a block "IF (I_AM_SLAVE)"
        wk_user_provided = (id%LWK_USER.NE.0)
        IF (id%LWK_USER.EQ.0) THEN
          itmp8 = 0_8
        ELSE IF (id%LWK_USER.GT.0) THEN
          itmp8= int(id%LWK_USER,8)
        ELSE
          itmp8 = -int(id%LWK_USER,8)* 1000000_8
        ENDIF
C       Incore: Check if the provided size is equal to that used during
C               facto (case of ITMP8/=0 and KEEP8(24)/=ITMP8)
C               But also check case of space not provided during solve
C               but was provided during facto
C                (case of ITMP8=0 and KEEP8(24)/=0)
        IF (keep(201).EQ.0) THEN  ! incore
C         Compare provided size with previous size
          IF (itmp8.NE.keep8(24)) THEN
C           -- error when reusing space allocated
            info(1) = -41
            info(2) = id%LWK_USER
            GOTO 99    ! jump to propinfo
                       ! (S is used in between and not allocated)
                       ! NO COMM must occur then before next propinfo
                       ! it happens in Mila's code but only with
                       ! KEEP(209) > 0
          ENDIF
        ELSE
          keep8(24)=itmp8
        ENDIF
C       KEEP8(24) holds the size of WK_USER provided by user.
C
        maxs = 0_8
        IF (wk_user_provided) THEN
           maxs = keep8(24)
           IF (maxs.LT. keep8(20)) THEN
                  info(1)= -11
                  ! MAXS should be increased by at least ITMP8
                  itmp8  = keep8(20)+1_8-maxs
                  CALL  mumps_set_ierror(itmp8, info(2))
           ENDIF
           IF (info(1) .GE. 0 ) id%S => id%WK_USER(1:keep8(24))
        ELSE IF (associated(id%S)) THEN
C          Avoid the use of "size(id%S)" because it returns
C          a default integer that may overflow. Also "size(id%S,kind=8)"
C          will only be available with Fortran 2003 compilers.
           maxs = keep8(23)
        ELSE
          ! S not allocated and WK_USER not provided ==> must be in OOC
          IF (keep(201).EQ.0) THEN ! incore
            WRITE(*,*) ' Working array S not allocated ',
     &                ' on entry to solve phase (in core) '
            CALL mumps_abort()
          ELSE
C         -- OOC and WK_USER not provided:
C            define size (S) and allocate it
C           ---- modify size of MAXS: in a simple
C           ---- system-based version, we want to
C           ---- use a small size for MAXS, to
C           ---- avoid the system pagecache to be
C           ---- polluted by 'our memory'
C
            IF ( keep(209).EQ.-1 .AND. workspace_minimal_preferred)
     &        THEN
C             We need space to load at least the largest factor
              maxs = keep8(20) + 1_8
            ELSE IF ( keep(209) .GE.0 ) THEN
C             Use suggested value of MAXS provided in KEEP(209)
              maxs = max(int(keep(209),8), keep8(20) + 1_8)
            ELSE
              maxs  = id%KEEP8(14) ! initial value: do not use more than
                               ! minimum (non relaxed) size of OOC facto
            ENDIF
C
            maxs = max(maxs, id%KEEP8(20)+1_8)
            ALLOCATE (id%S(maxs), stat = allocok)
            keep8(23)=maxs
            IF ( allocok .GT. 0 ) THEN
              IF (lpok) THEN
                 WRITE(lp,*) id%MYID,': problem allocation of S ',
     &                'at solve'
              ENDIF
              info(1) = -13
              CALL mumps_set_ierror(maxs, info(2))
              NULLIFY(id%S)
              keep8(23)=0_8
            ENDIF
            nb_bytes = nb_bytes + keep8(23) * k35_8
            nb_bytes_max = max(nb_bytes_max,nb_bytes)
C           --- end of OOC case
          ENDIF
C         -- end of id%S already associated
        ENDIF
C
C       On the slaves, S is divided as follows:
C       S(1..LA) holds the factors,
C       S(LA+1..MAXS) is free workspace
        IF(keep(201).EQ.0)THEN
           la  = keep8(31)
        ELSE
C          MAXS has normally be dimensionned to store only factors.
           la = maxs
           IF(maxs.GT.keep8(31)+keep8(20)*int(keep(107)+1,8))THEN
C            If we have a very large MAXS, the size reserved for
C            loading the factors into memory does not need to exceed the
C            total size of factors. The (KEEP8(20)*(KEEP(107)+1)) term
C            is here in order to ensure that even with round-off
C            problems (linked to the number of solve zones) factors can
C            all be stored in-core
             la=keep8(31)+keep8(20)*int(keep(107)+1,8)
           ENDIF
        ENDIF
C
C       We need to allocate a workspace of size LWCB8 for the solve phase.
C       Either it is available at the end of MAXS, or we perform a
C       dynamic allocation.
        IF ( maxs-la .GT. lwcb8_min ) THEN
           lwcb8 = maxs - la
           work_wcb => id%S(la+1_8:la+lwcb8)
           work_wcb_allocated=.false.
        ELSE
           lwcb8 = lwcb8_min
           ALLOCATE(work_wcb(lwcb8), stat = allocok)
           IF (allocok < 0 ) THEN
              info(1)=-13
              CALL mumps_set_ierror(lwcb8,info(2))
           ENDIF
           work_wcb_allocated=.true.
           nb_bytes = nb_bytes + lwcb8*k35_8
           nb_bytes_max = max(nb_bytes_max,nb_bytes)
        ENDIF
      ENDIF ! I_AM_SLAVE
C -----------------------------------
  99  CONTINUE
      CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
      IF (info(1) < 0) GOTO 90
C -----------------------------------
      IF ( i_am_slave ) THEN
        IF (keep(201).GT.0) THEN
          CALL zmumps_init_fact_area_size_s(la)
C         -- This includes thread creation
C         -- for asynchronous strategies
          CALL zmumps_ooc_init_solve(id)
          is_init_ooc_done = .true.
        ENDIF ! KEEP(201).GT.0
      ENDIF
C
      CALL mumps_propinfo( icntl(1), info(1),
     &                     id%COMM,id%MYID)
      IF (info(1) < 0) GOTO 90
C
      IF (i_am_slave) THEN
        IF (keep(485).EQ.1) THEN
          IF (.NOT. (associated(id%FDM_F_ENCODING))) THEN
            WRITE(*,*) "Internal error 18 in ZMUMPS_SOL_DRIVER"
            CALL mumps_abort()
          ENDIF
          IF (.NOT. (associated(id%BLRARRAY_ENCODING))) THEN
            WRITE(*,*) "Internal error 19 in ZMUMPS_SOL_DRIVER"
            CALL mumps_abort()
          ENDIF
C         Access to OOC data in module during solve
          CALL mumps_fdm_struc_to_mod('F',id%FDM_F_ENCODING)
          CALL zmumps_blr_struc_to_mod(id%BLRARRAY_ENCODING)
          is_lr_mod_to_struc_done = .true.
        ENDIF
      ENDIF
      IF (id%MYID.EQ.master) THEN
        IF ( prokg )  THEN
          WRITE( mpg, 150 )
C             ICNTL(35) should not been accessed during SOLVE thus
C             print KEEP(486) value set during factorization
     &        id%NRHS, nbrhs, icntl(9), icntl(10), icntl(11),
     &        icntl(20), icntl(21), icntl(30), keep(486)
          IF (keep(111).NE.0) THEN
            WRITE (mpg, 151) keep(111)
          ENDIF
          IF (keep(221).NE.0) THEN
            WRITE (mpg, 152) keep(221)
          ENDIF
          IF (keep(252).GT.0) THEN   ! Fwd during facto
            WRITE (mpg, 153) keep(252)
          ENDIF
        ENDIF
C
C     ====================================
C       Define LSCAL, ICNTL10 and ICNTL11
C     ====================================
C
        lscal = (((keep(52) .GT. 0) .AND. (keep(52) .LE. 8)) .OR. (
     &    keep(52) .EQ. -1) .OR. keep(52) .EQ. -2)
        icntl10 = icntl(10)
        icntl11 = icntl(11)
C       Values of ICNTL(11) out of range
        IF ((icntl11 .LT. 0).OR.(icntl11 .GE. 3)) THEN
           icntl11 = 0
           IF (prokg) WRITE(mpg,'(A)')
     &    ' WARNING: ICNTL(11) out of range'
        ENDIF
        postpros = .false.
        IF (icntl11.NE.0 .OR. icntl10.NE.0) THEN
          postpros = .true.
C       FORBID ERROR ANALYSIS AND ITERATIVE REFINEMENT
C       if there are options that are not compatible
          IF (keep(111).NE.0) THEN
C       IF WE RETURN A NULL SPACE BASIS or compute entries in A-1
C        of Fwd in facto
C       -When only one columns of A-1 is requested then
C        we could try to reactivate IR even if
C          -code need be updated
C          -accuracy could be # when one or more columns are requested
              IF (prokg) WRITE(mpg,'(A,A)')
     &       ' WARNING: Incompatible features: null space basis ',
     &                ' and Iter. Ref and/or Err. Anal.'
              postpros = .false.
           ELSE IF ( keep(237) .NE.0 ) THEN
              IF (prokg) WRITE(mpg,'(A,A)')
     &       ' WARNING: Incompatible features: AM1',
     &                ' and Iter. Ref and/or Err. Anal.'
              postpros = .false.
           ELSE IF ( keep(252) .NE.0 ) THEN
              IF (prokg) WRITE(mpg,'(A,A)')
     &       ' WARNING: Incompatible features: Fwd in facto ',
     &                ' and Iter. Ref and/or Err. Anal.'
              postpros = .false.
           ELSE IF (keep(221).NE.0) THEN
C       Forbid error analysis and iterative refinement
C       in case of reduced rhs/solution
              IF (prokg) WRITE(mpg,'(A,A)')
     &       ' WARNING: Incompatible features: reduced RHS ',
     &       '          and Iter. Ref and/or Err. Anal.'
              postpros = .false.
            ELSE IF  (nbrhs.GT. 1 .OR. icntl(21) .GT. 0) THEN
C          Forbid error analysis and iterative refinement if
C       the solution is distributed or
C       in the case where nrhs > 1
              IF (prokg) WRITE(mpg,'(A,A)')
     &       ' WARNING:  Incompatible features: nrhs>1 or distrib sol',
     &       '           and Iter. Ref and/or Err. Anal.'
              postpros = .false.
            ELSE IF ( keep(248) .EQ. -1 ) THEN
C          Forbid error analysis and iterative refinement
C          in case of distributed RHS
              IF (prokg) WRITE(mpg,'(A,A)')
     &       ' WARNING:  Incompatible features: distrib rhs',
     &       '           and Iter. Ref and/or Err. Anal.'
              postpros = .false.
            ENDIF  
            IF (.NOT.postpros) THEN
              icntl11 = 0
              icntl10 = 0
            ENDIF
        ENDIF
C   Write a warning.
        IF ((icntl(10) .NE. 0) .AND. (icntl10 .EQ. 0)) THEN
            IF (prokg) WRITE(mpg,'(A)')
     &    ' WARNING: ICNTL(10) treated as if set to 0 '
        ENDIF
        IF ((icntl(11) .NE. 0)
     &        .AND.(icntl11 .EQ. 0)) THEN
            IF (prokg) WRITE(mpg,'(A)')
     &    ' WARNING: ICNTL(11) treated as if set to 0 '
        ENDIF
C     -- end of test master
      END IF
       CALL mpi_bcast(postpros,1,mpi_logical,master,
     &               id%COMM,ierr)
C  We need the original matrix only in the case of 
C  we want to perform IR or Error Analysis, i.e. if
C  POSTPros = TRUE
      mat_alloc_loc = 0
      IF ( postpros ) THEN
       mat_alloc_loc = 1
C      Check if the original matrix has been allocated.
        IF ( keep(54) .EQ. 0 ) THEN
C        The original  matrix is centralized
         IF ( id%MYID .eq. master ) THEN
            IF (keep(55).eq.0) THEN 
C            Case of matrix assembled centralized
             IF (.NOT.associated(id%A) .OR. 
     &          (.NOT.associated(id%IRN)) .OR. 
     &                  ( .NOT.associated(id%JCN))) THEN
               IF (prokg) WRITE(mpg,'(A)')
     &           ' WARNING: original centralized assembled', 
     &           ' matrix is not allocated '
                mat_alloc_loc = 0
             ENDIF
            ELSE 
C            Case of matrix in elemental format
             IF (.NOT.associated(id%A_ELT).OR. 
     &           .NOT.associated(id%ELTPTR).OR.
     &           .NOT.associated(id%ELTVAR)) THEN
              IF (prokg) WRITE(mpg,'(A)')
     &        ' WARNING: original elemental matrix is not allocated '
              mat_alloc_loc = 0
             ENDIF
            ENDIF
         ENDIF  !end master, centralized matrix
        ELSE                                     
C        The original  matrix is assembled distributed
         IF ( i_am_slave .AND. (id%KEEP8(29) .GT. 0_8) ) THEN
C        If MAT_ALLOC_LOC = 1 the local distributed matrix is
C        allocated, otherwise MAT_ALLOC_LOC = 0
           IF ((.NOT.associated(id%A_loc)) .OR. 
     &        (.NOT.associated(id%IRN_loc)) .OR. 
     &        (.NOT.associated(id%JCN_loc))) THEN
             IF (prokg) WRITE(mpg,'(A)')
     &           ' WARNING: original distributed assembled', 
     &           '  matrix is not allocated '
             mat_alloc_loc = 0
           ENDIF
         ENDIF
        ENDIF ! end test allocation matrix (keep(54))
      ENDIF ! POSTPros
      CALL mpi_reduce( mat_alloc_loc, mat_alloc, 1,
     &                 mpi_integer,
     &                 mpi_min, master, id%COMM, ierr)
      IF ( id%MYID .eq. master ) THEN
        IF (mat_alloc.EQ.0) THEN
              postpros = .false.
              icntl11 = 0
              icntl10 = 0
C   Write a warning.
           IF ((icntl(10) .NE. 0) .AND. (icntl10 .EQ. 0)) THEN
            IF (prokg) WRITE(mpg,'(A)')
     &    ' WARNING: ICNTL(10) treated as if set to 0 '
           ENDIF
           IF ((icntl(11) .EQ. 1).OR.(icntl(11) .EQ. 2) 
     &        .AND.(icntl11 .EQ. 0)) THEN
            IF (prokg) WRITE(mpg,'(A)')
     &    ' WARNING: ICNTL(11) treated as if set to 0 '
          ENDIF
        ENDIF
        IF (postpros) THEN
            ALLOCATE(saverhs(id%N*nbrhs),stat = allocok)
            IF ( allocok .GT. 0 ) THEN
              IF (lpok) THEN
                 WRITE(lp,*) id%MYID,
     &                ':Problem in solve: error allocating SAVERHS'
              ENDIF
              info(1) = -13
              info(2) = id%N*nbrhs
            END IF
            nb_bytes = nb_bytes + int(size(saverhs),8)*k35_8
            nb_bytes_max = max(nb_bytes_max,nb_bytes)
        ENDIF
C
C     Forbid entries in a-1, in case of null space computations
c
        IF (keep(237).NE.0 .AND.keep(111).NE.0) THEN
C         Ignore ENTRIES IN A-1 in case we compute
C         vectors of the null space (KEEP(111)).NE.0.)
C         We should still allocate IRHS_SPARSE
          IF (prokg) WRITE(mpg,'(A)')
     &    ' WARNING: KEEP(237) treated as if set to 0 (null space)'
          keep(237)=0
        ENDIF
C     -- end of test master
      END IF
      CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
      IF (info(1) .LT.0 ) GOTO 90
C     --------------------------------------------------
C     Broadcast information to have all processes do the
C     same thing (error analysis/iterative refinements/
C     scaling/distribution of solution)
C     --------------------------------------------------
      CALL mpi_bcast(icntl10,1,mpi_integer,master,
     &               id%COMM,ierr)
      CALL mpi_bcast(icntl11,1,mpi_integer,master,
     &               id%COMM,ierr)
      CALL mpi_bcast(icntl21,1,mpi_integer,master,
     &               id%COMM,ierr)
      CALL mpi_bcast(postpros,1,mpi_logical,master,
     &               id%COMM,ierr)
      CALL mpi_bcast(lscal,1,mpi_logical,master,
     &               id%COMM,ierr)
      CALL mpi_bcast(keep(237),1,mpi_integer,master,
     &               id%COMM,ierr)
C     KEEP(248)==1 if not_NullSpace (KEEP(111)=0)
C           and sparse RHS on input (id%ICNTL(20)/KEEP(248)==1)
C     (KEEP(248)==1 implies KEEP(111) = 0, otherwise error was raised)
C     We cant thus isolate the case of
C     sparse RHS associated to Null space computation because
C     in this case preparation is different since
C        -we skip the forward step and
C        -the pattern of the RHS
C         of the bwd is related to null pivot indices found and not
C         to information contained in the sparse rhs input format.
      do_permute_rhs = (keep(242).NE.0)
C      apply interleaving in parallel (FOR A-1 or Null space only)
      IF ( (id%NSLAVES.GT.1) .AND. (keep(243).NE.0)
     &        )  THEN
C       -- Option to interleave RHS only makes sense when
C       -- A-1 option is on or Null space compution are on
C          (note also that KEEP(243).NE.0 only when PERMUTE_RHS is on)
        IF ((keep(237).NE.0).or.(keep(111).GT.0)) THEN
           interleave_par= .true.
        ELSE
          IF (prokg) THEN
            write(mpg,*) ' Warning incompatible options ',
     &      ' interleave RHS reset to false '
          ENDIF
        ENDIF
      ENDIF
C       --------------------------------------
C       Compute an upperbound of message size
C       for forward and backward solutions:
C       --------------------------------------
        msg_max_bytes_solve8 = int(( 4 + keep(133) ) * keep(34),8) +
     &                         int(keep(133)*keep(35),8) * int(nbrhs,8)
     &  + int(16*keep(34),8) ! for request id, pointer to next + safety
C       Note that
        IF ( msg_max_bytes_solve8 .GT.
     &                int(huge(i4),8)) THEN
          info(1) = -18
          info(2) = ( huge(i4) -
     &    ( 16 + 4 + keep(133) ) ) / 
     &    ( keep(133) * keep(35) )
        ENDIF
        IF (info(1) .LT.0 ) GOTO 111
        msg_max_bytes_solve = int(msg_max_bytes_solve8)
C       ------------------------------------------
C       Compute an upperbound of message size
C       for ZMUMPS_GATHER_SOLUTION. Except
C       possibly on the non working host, it
C       should be smaller than MSG_MAX_BYTES_SOLVE
C       ------------------------------------------
        IF (keep(237).EQ.0) THEN
C         Note that for ZMUMPS_GATHER_SOLUTION LBUFR buffer should
C         be larger that MAX_inode(NPIV))*NBRHS + NPIV
C         which is covered by next formula since KMAX_246_247  is larger
C         than  MAX_inode(NPIV))
C                   2 integers packed (npiv and termination)
C         Note that MSG_MAX_BYTES_GTHRSOL < MSG_MAX_BYTES_SOLVE
C         so that it should not overflow
          kmax_246_247 = max(keep(246),keep(247))
          msg_max_bytes_gthrsol =  ( 2 + kmax_246_247 ) * keep(34) +
     &                             kmax_246_247 * nbrhs * keep(35)
        ELSE IF (icntl21.EQ.0) THEN
C         Each message from a slave is of size max 4:
C            2 integers  : I,J
C            1 complex   : (Aij)-1
C            1 terminaison
          msg_max_bytes_gthrsol =  (  3  * keep(34) + keep(35) )
        ELSE
C         Not needed in case of distributed solution and A-1
C         because the entries of A −1 are 
C         returned in RHS SPARSE on the host.
          msg_max_bytes_gthrsol =  0
        ENDIF
C       The buffer is used both for solve and for ZMUMPS_GATHER_SOLUTION
        lbufr_bytes = max(msg_max_bytes_solve, msg_max_bytes_gthrsol)
        tsize = int(min(100_8*int(msg_max_bytes_gthrsol,8),
     &              10000000_8))
        lbufr_bytes = max(lbufr_bytes,tsize)
        lbufr = ( lbufr_bytes + keep(34) - 1 ) / keep(34)
        ALLOCATE (bufr(lbufr),stat=allocok)
        IF ( allocok .GT. 0 ) THEN
            IF (lpok) THEN
               WRITE(lp,*) id%MYID,
     &              ' Problem in solve: error allocating BUFR'
            ENDIF
            info(1) = -13
            info(2) = lbufr
            GOTO 111
        ENDIF
        nb_bytes = nb_bytes + int(size(bufr),8)*k34_8
        nb_bytes_max = max(nb_bytes_max,nb_bytes)
      IF ( i_am_slave .AND. id%NSLAVES .GT. 1 ) THEN
C       ------------------------------------------------------
C       Dimension send buffer for small integers, e.g. TRACINE
C       ------------------------------------------------------
        zmumps_lbuf_int = ( 20 + id%NSLAVES * id%NSLAVES  * 4 )
     &                 * keep(34)
        CALL zmumps_buf_alloc_small_buf( zmumps_lbuf_int, ierr )
        IF ( ierr .NE. 0 ) THEN
          info(1) = -13
          info(2) = zmumps_lbuf_int
          IF ( lpok) THEN
            WRITE(lp,*) id%MYID,
     &      ':Error allocating small Send buffer:IERR=',ierr
          END IF
          GOTO 111
        END IF
C
C       ---------------------------------------
C       Dimension cyclic send buffer for normal
C       messages, based on largest message
C       size during forward and backward solves
C       ---------------------------------------
C       Compute buffer size in BYTES (ZMUMPS_LBUF)
C       using integer8 in ZMUMPS_LBUF_8 
C       then convert it in integer4 and bound it to largest integer value
C
        zmumps_lbuf_8 = 
     &       (int(msg_max_bytes_solve,8)+2_8*int(keep(34),8))*
     &        int(id%NSLAVES,8)
C       Avoid buffers larger than 100 Mbytes ...
        zmumps_lbuf_8 = min(zmumps_lbuf_8, 100000000_8)
C       ... as long as we can send messages to at least 3
C       destinations simultaneously
        zmumps_lbuf_8 = max(zmumps_lbuf_8,
     &      int((msg_max_bytes_solve+2*keep(34)),8) * 
     &      int(min(id%NSLAVES,3),8) )
        zmumps_lbuf_8 = zmumps_lbuf_8 + 2_8*int(keep(34),8)
C       Convert to integer and bound it to largest 32-bit integer
C       and suppress 10 integers (one should be enough!) 
C       to enable computation of integer size.
        zmumps_lbuf_8 = min(zmumps_lbuf_8, 
     &                      int(huge(i4),8)
     &                      - 10_8*int(keep(34),8)
     &                     )
        zmumps_lbuf   = int(zmumps_lbuf_8, kind(zmumps_lbuf))
        CALL zmumps_buf_alloc_cb( zmumps_lbuf, ierr )
        IF ( ierr .NE. 0 ) THEN
          info(1) = -13
          info(2) = zmumps_lbuf/keep(34) + 1
          IF ( lpok) THEN
            WRITE(lp,*) id%MYID,
     &      ':Error allocating Send buffer:IERR=', ierr
          END IF
          GOTO 111
        END IF
C
C
C     -- end of I am slave
      ENDIF
C
      IF ( postpros )  THEN
C       When Iterative refinement of error analysis requested
C       Allocate RHS_IR on slave processors
C       (note that on MASTER RHS_IR points to RHS)
        IF ( id%MYID .NE. master ) THEN
C
          ALLOCATE(rhs_ir(id%N),stat=ierr)
          nb_bytes = nb_bytes + int(size(rhs_ir),8)*k35_8
          nb_bytes_max = max(nb_bytes_max,nb_bytes)
          IF ( ierr .GT. 0 ) THEN
            info(1)=-13
            info(2)=id%N
            IF (lpok) THEN
               WRITE(lp,*) 'ERROR while allocating RHS on a slave'
            ENDIF
            GOTO 111
          END IF
        ELSE
          rhs_ir=>id%RHS
        ENDIF
      ENDIF
C
C     Parallel A-1 or General sparse and 
C     exploit sparsity between columns
      do_nbsparse = ( ( (keep(237).NE.0).OR.(keep(235).NE.0) )
     &                  .AND.
     &                ( keep(497).NE.0 ) 
     &              )
      IF ( i_am_slave ) THEN
        IF(do_nbsparse) THEN
c         --- ALLOCATE outside loop RHS_BOUNDS is needed
          lptr_rhs_bounds = 2*keep(28)
          ALLOCATE(rhs_bounds(lptr_rhs_bounds), stat=ierr)
          IF (ierr.GT.0) THEN
            info(1)=-13
            info(2)=lptr_rhs_bounds
            IF (lpok) THEN
               WRITE(lp,*) 'ERROR while allocating RHS_BOUNDS on',
     &              ' a slave'
            ENDIF
            GOTO 111
          END IF
          nb_bytes = nb_bytes +
     &        int(size(rhs_bounds),8)*k34_8
          nb_bytes_max = max(nb_bytes_max,nb_bytes)
          ptr_rhs_bounds => rhs_bounds
        ELSE
          lptr_rhs_bounds = 1
          ptr_rhs_bounds => idummy_target
        ENDIF
      ENDIF
C     --------------------------------------------------
      IF ( i_am_slave ) THEN
        IF ((keep(221).EQ.2 .AND. keep(252).EQ.0)) THEN
C          -- RHSCOMP must have been allocated in
C          -- previous solve step (with option KEEP(221)=1)
           IF (.NOT.associated(id%RHSCOMP)) THEN
             info(1) = -35
             info(2) = 1
             GOTO 111
           ENDIF
C          IF ((KEEP(248).EQ.0) .OR. (id%NRHS.EQ.1)) THEN
C          POSINRHSCOMP_ROW/COL are meaningful and could even be reused
           IF (.NOT.associated(id%POSINRHSCOMP_ROW) ) ! .OR.
!    &        .NOT.(id%POSINRHSCOMP_COL_ALLOC))
     &     THEN
             info(1) = -35
             info(2) = 2
             GOTO 111
           ENDIF
           IF (.not.id%POSINRHSCOMP_COL_ALLOC) THEN
C             POSINRHSCOMP_COL that is kept from
C             previous call to solve must then (already)
C             point to id%POSINRHSCOMP_ROW
              id%POSINRHSCOMP_COL => id%POSINRHSCOMP_ROW
           ENDIF
        ELSE
C         ----------------------
C         Allocate POSINRHSCOMP_ROW/COL
C         ----------------------
C         The size of POSINRHSCOMP arrays
C         does not depend on the block of RHS
C         POSINRHSCOMP_ROW/COL are initialized in the loop of RHS
          IF (associated(id%POSINRHSCOMP_ROW)) THEN
            nb_bytes = nb_bytes -
     &          int(size(id%POSINRHSCOMP_ROW),8)*k34_8
            DEALLOCATE(id%POSINRHSCOMP_ROW)
          ENDIF
          ALLOCATE (id%POSINRHSCOMP_ROW(id%N), stat = allocok)
          IF ( allocok .GT. 0 ) THEN
             info(1)=-13
             info(2)=id%N
             GOTO 111
          END IF
          nb_bytes = nb_bytes +
     &          int(size(id%POSINRHSCOMP_ROW),8)*k34_8
          nb_bytes_max = max(nb_bytes_max,nb_bytes)
          IF (id%POSINRHSCOMP_COL_ALLOC) THEN
            nb_bytes = nb_bytes -
     &          int(size(id%POSINRHSCOMP_COL),8)*k34_8
            DEALLOCATE(id%POSINRHSCOMP_COL)
            NULLIFY(id%POSINRHSCOMP_COL)
            id%POSINRHSCOMP_COL_ALLOC = .false.
          ENDIF
C
          IF ((keep(50).EQ.0).OR.keep(237).NE.0) THEN
           ALLOCATE (id%POSINRHSCOMP_COL(id%N), stat = allocok)
           IF ( allocok .GT. 0 ) THEN
             info(1)=-13
             info(2)=id%N
             GOTO 111
           END IF
           id%POSINRHSCOMP_COL_ALLOC = .true.
           nb_bytes = nb_bytes +
     &          int(size(id%POSINRHSCOMP_COL),8)*k34_8
           nb_bytes_max = max(nb_bytes_max,nb_bytes)
          ELSE
C          Do no allocate POSINRHSCOMP_COL
           id%POSINRHSCOMP_COL => id%POSINRHSCOMP_ROW
           id%POSINRHSCOMP_COL_ALLOC = .false.
          ENDIF
          IF (keep(221).NE.2) THEN
C         -- only in the case of bwd after reduced RHS
C         -- we have to keep "old" RHSCOMP
            IF (associated(id%RHSCOMP)) THEN
             nb_bytes = nb_bytes - id%KEEP8(25)*k35_8
             DEALLOCATE(id%RHSCOMP)
             NULLIFY(id%RHSCOMP)
             id%KEEP8(25)=0_8
            ENDIF
          ENDIF
        ENDIF
C       ---------------------------
C       Allocate local workspace
C       for the solve (ZMUMPS_SOL_C)
C       ---------------------------
        liwk_solve = 2 * keep(28) + id%NA(1)+1
        liwk_ptracb= keep(28) 
C       KEEP(228)+1 temporary integer positions
C       will be needed in ZMUMPS_SOL_S
        IF (keep(201).EQ.1) THEN
          liwk_solve = liwk_solve + keep(228) + 1
        ELSE
C         Reserve 1 position to pass array of size 1 in routines
          liwk_solve = liwk_solve + 1
        ENDIF
        ALLOCATE ( iwk_solve(liwk_solve), 
     &             ptracb(liwk_ptracb), stat = allocok )
        IF (allocok .GT. 0 ) THEN
         info(1)=-13
         info(2)=liwk_solve + liwk_ptracb*keep(10)
         GOTO 111
        END IF
        nb_bytes = nb_bytes + int(liwk_solve,8)*k34_8 +
     &             int(liwk_ptracb,8)*k34_8 *int(keep(10),8)
        nb_bytes_max = max(nb_bytes_max,nb_bytes)
C       array IWCB used temporarily to hold
C       indices of a front unpacked from a message
C       and to stack (potentially in a recursive call)
C       headers of size 2 positions of CB blocks.
        liwcb =  20*nb_k133*2 + keep(133)
        ALLOCATE ( iwcb( liwcb), stat = allocok )
        IF (allocok .GT. 0 ) THEN
         info(1)=-13
         info(2)=liwcb
         GOTO 111
        END IF
        nb_bytes = nb_bytes + int(liwcb,8)*k34_8
        nb_bytes_max = max(nb_bytes_max,nb_bytes)
C
C       -- Code for a slave
C       -----------
C       Subdivision
C       of array IS
C       -----------
        liw = keep(32)
C       Define a work array of size maximum global frontal
C       size (KEEP(133)) for the call to ZMUMPS_SOL_C
C       This used to be of size id%N.
        ALLOCATE(srw3(keep(133)), stat = allocok )
        IF ( allocok .GT. 0 ) THEN
          info(1)=-13
          info(2)=keep(133)
          GOTO 111
        END IF
        nb_bytes = nb_bytes + int(size(srw3),8)*k35_8
        nb_bytes_max = max(nb_bytes_max,nb_bytes)
C     -----------------
C     End of slave code
C     -----------------
      ELSE
C       I am the master with host not working
C
C       LIW is used on master when calling
C       the routine ZMUMPS_GATHER_SOLUTION.
        liw=0
      END IF
C
C     Precompute inverse of UNS_PERM outside loop
      IF (allocated(uns_perm_inv)) DEALLOCATE(uns_perm_inv)
      uns_perm_inv_needed_inmainloop = .false.
      IF ( ( id%MYID .eq. master.AND.(keep(23).GT.0) .AND.
     &         (mtype .NE. 1).AND.(keep(248).NE.0)
     &       )
C          Permute UNS_PERM on master only with
C          sparse RHS (KEEP(248).NE.0 ) when AT x = b is solved
     &      .OR. ( keep(237).NE.0 .AND. keep(23).NE.0  )
C          When A-1 is active and when the matrix is unsymmetric
C          and a column permutation has been applied (Max transversal)
C          then  we have performed a
C          factorization of a column permuted matrix AQ = LU.
C          In this case,
C          the permuted entry must be used to select the target
C          entries for the BWD (note that a diagonal entry of A-1
C          is not anymore a diagonal of AQ. Thus a diagonal
C          of A-1 does not correspond to the same path
C          in the tree during FWD and BWD steps when MAXTRANS is on
C          and permutation is not identity.)
C          Note that the inverse permutation
C          UNS_PERM_INV needs to be allocated on each proc
C          since it is used in ZMUMPS_SOL_C routine for pruning.
C          It is allocated only once and its allocation has been
C          migrated outside the blocking on the right hand sides.
     & ) THEN
            uns_perm_inv_needed_inmainloop = .true.
      ENDIF
      uns_perm_inv_needed_befmainloop = .false.
      IF ( keep(23) .GT.0 .AND.
     &            mtype .NE. 1 .AND. keep(248).EQ.-1 ) THEN
C          Similar to sparse RHS case, we need to modify IRHS_loc
C          indices in the distributed RHS case. However, we need
C          UNS_PERM_INV on all processors. But only before theC
C          main loop on the RHS blocks.
           uns_perm_inv_needed_befmainloop = .true.
      ENDIF
      IF ( uns_perm_inv_needed_inmainloop .OR.
     &     uns_perm_inv_needed_befmainloop ) THEN
           ALLOCATE(uns_perm_inv(id%N),stat=allocok)
           if (allocok .GT.0 ) THEN
             info(1)=-13
             info(2)=id%N
             GOTO 111
           endif
           nb_bytes = nb_bytes + int(id%N,8)*k34_8
           nb_bytes_max = max(nb_bytes_max,nb_bytes)
           IF (id%MYID.EQ.master) THEN
C           Build inverse permutation
            DO i = 1, id%N
              uns_perm_inv(id%UNS_PERM(i))=i
            ENDDO
           ENDIF
C
      ELSE
           ALLOCATE(uns_perm_inv(1), stat=allocok)
           if (allocok .GT.0 ) THEN
             info(1)=-13
             info(2)=1
             GOTO 111
           endif
           nb_bytes = nb_bytes + 1_8*k34_8
           nb_bytes_max = max(nb_bytes_max,nb_bytes)
      ENDIF
C
 111  CONTINUE
#if defined(V_T)
      CALL vtend(glob_comm_ini,ierr)
#endif
C
C     Synchro point + Broadcast of errors
C
C     Propagate errors
      CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
      IF (info(1) .LT.0 ) GOTO 90
C
C     UNS_PERM_INV needed on slaves:
      IF ( keep(23).NE.0 .AND.
     &     ( keep(237).NE.0 .OR.
     &     ( mtype.NE.1 .AND. keep(248).EQ.-1 ) ) ) THEN
C       Broadcast UNS_PERM_INV
        CALL mpi_bcast( uns_perm_inv,id%N,mpi_integer,master,
     &                  id%COMM,ierr )
      ENDIF
C     -------------------------------
C     BEGIN
C     Preparation for distributed RHS
C     -------------------------------
      IF (i_am_slave .AND. keep(248).EQ.-1) THEN
C       Distributed RHS case
        ALLOCATE(map_rhs_loc(max(id%Nloc_RHS,1)), stat=allocok)
        IF (allocok .GT. 0) THEN
          id%INFO(1)=-13
          id%INFO(2)=max(id%Nloc_RHS,1)
          GOTO 20
        ENDIF
        nb_bytes = nb_bytes + max(int(id%Nloc_RHS,8),1_8)*k34_8
      ENDIF
C     MAP_RHS_loc will be built in the main
C     loop, when processing the first block.
C     It requires POSINRHSCOMP to be built.
      build_rhsmapinfo = .true.
 20   CONTINUE
      CALL mumps_propinfo( icntl(1), info(1),
     &                         id%COMM,id%MYID)
      IF ( info(1) .LT.0 ) GOTO 90
C     In case of Unsymmetric column permutation and
C     transpose system, use MUMPS internal indices
C     for IRHS_loc_PTR. Done before scaling since
C     scaling is on permuted matrix
      IF ( i_am_slave .AND. keep(23).GT.0 .AND. keep(248).EQ.-1
     &     .AND. mtype.NE.1 ) THEN
          IF (id%Nloc_RHS .GT. 0) THEN
            ALLOCATE(irhs_loc_ptr(id%Nloc_RHS),stat=allocok)
            IF (allocok.GT.0) THEN
              info(1)=-13
              info(2)=id%Nloc_RHS
              GOTO 25
            ENDIF
            irhs_loc_ptr_allocated = .true.
            nb_bytes = nb_bytes + max(int(id%Nloc_RHS,8),1_8)*k34_8
            nb_bytes_max = max(nb_bytes_max,nb_bytes)
            DO i=1, id%Nloc_RHS
              IF (id%IRHS_loc(i).GE.1 .AND. id%IRHS_loc(i).LE.id%N)
     &          THEN
                irhs_loc_ptr(i)=uns_perm_inv(id%IRHS_loc(i))
              ELSE
C               Keep track of out-of range entries
                irhs_loc_ptr(i)=id%IRHS_loc(i)
              ENDIF
            ENDDO
          ENDIF
      ENDIF
C     Check if UNS_PERM_INV still needed
C     to free memory
      IF (uns_perm_inv_needed_befmainloop .AND.
     & .NOT. uns_perm_inv_needed_inmainloop) THEN
        nb_bytes = nb_bytes - int(size(uns_perm_inv),8)*k34_8
        DEALLOCATE(uns_perm_inv)
        ALLOCATE(uns_perm_inv(1))  ! to posibly pass it as an argument
        nb_bytes = nb_bytes + k34_8
      ENDIF
      IF (lscal .AND. id%KEEP(248).EQ.-1) THEN
C       Scaling done based on original indices
C       provided by user
        IF (mtype == 1) THEN
C         No transpose
          scaling_data_dr%SCALING=>id%ROWSCA
        ELSE
C         Transpose
          scaling_data_dr%SCALING=>id%COLSCA
        ENDIF
        CALL zmumps_set_scaling_loc( scaling_data_dr, id%N,
     &       irhs_loc_ptr(1), id%Nloc_RHS,
     &       id%COMM, id%MYID, i_am_slave, master,
     &       nb_bytes, nb_bytes_max, k16_8, lp, lpok,
     &       icntl(1), info(1) )
      ENDIF
C     -------------------------------
C     END
C     Preparation for distributed RHS
C     -------------------------------
 25   CONTINUE
      CALL mumps_propinfo( icntl(1), info(1),
     &                         id%COMM,id%MYID)
      IF ( info(1) .LT.0 ) GOTO 90
C     -------------------------------------
C     BEGIN
C     Preparation for distributed solution
C     -------------------------------------
      IF ( icntl21==1 ) THEN
        IF (lscal) THEN
C         In case of scaling we will need to scale
C         back the sol. Put the values of the scaling
C         arrays needed to do that on each processor.
          IF (id%MYID.NE.master) THEN
            IF (mtype == 1) THEN
              ALLOCATE(id%COLSCA(id%N),stat=allocok)
            ELSE
              ALLOCATE(id%ROWSCA(id%N),stat=allocok)
            ENDIF
            IF (allocok > 0) THEN
              IF (lpok) THEN
               WRITE(lp,*) 'Error allocating temporary scaling array'
              ENDIF
              info(1)=-13
              info(2)=id%N
              GOTO 37
            ENDIF
            nb_bytes = nb_bytes + int(id%N,8)*k16_8
            nb_bytes_max = max(nb_bytes_max,nb_bytes)
          ENDIF ! MYID .NE. MASTER
 37       CALL mumps_propinfo( icntl(1), info(1),
     &                         id%COMM,id%MYID)
          IF (info(1) .LT.0 ) GOTO 90
          IF (i_am_slave) THEN
            ALLOCATE(scaling_data_sol%SCALING_LOC(id%KEEP(89)),
     &               stat=allocok)
            IF (allocok > 0) THEN
              IF (lpok) THEN
                WRITE(lp,*) 'Error allocating local scaling array'
              ENDIF
              info(1)=-13
              info(2)=id%KEEP(89)
              GOTO 38
            ENDIF
            nb_bytes = nb_bytes + int(id%KEEP(89),8)*k16_8
            nb_bytes_max = max(nb_bytes_max,nb_bytes)
          ENDIF ! I_AM_SLAVE
 38       CONTINUE
          CALL mumps_propinfo( icntl(1), info(1),
     &                         id%COMM,id%MYID)
          IF (info(1) .LT.0 ) THEN
            GOTO 90
          ENDIF
          IF (mtype == 1) THEN
              CALL mpi_bcast(id%COLSCA(1),id%N,
     &                       mpi_double_precision,master,
     &                       id%COMM,ierr)
              scaling_data_sol%SCALING=>id%COLSCA
          ELSE
              CALL mpi_bcast(id%ROWSCA(1),id%N,
     &                       mpi_double_precision,master,
     &                       id%COMM,ierr)
              scaling_data_sol%SCALING=>id%ROWSCA
          ENDIF
        ENDIF ! LSCAL
        IF ( i_am_slave ) THEN
C         ----------------------
C         Prepare ISOL_loc array
C         ----------------------
          liw_passed=max(1,liw)
C         Only called if more than 1 pivot
C         was eliminated by the processor.
C         Note that LSOL_loc >= KEEP(89)
          IF (keep(89) .GT. 0) THEN
            CALL zmumps_distsol_indices( mtype, id%ISOL_loc(1),
     &               id%PTLUST_S(1),
     &               id%KEEP(1),id%KEEP8(1),
     &               id%IS(1), liw_passed,id%MYID_NODES,
     &               id%N, id%STEP(1), id%PROCNODE_STEPS(1),
     &               id%NSLAVES, scaling_data_sol, lscal
C                    For checking only
     &               , (keep(248).EQ.-1), irhs_loc_ptr(1), id%Nloc_RHS
     &               )
          ENDIF 
          IF (id%MYID.NE.master .AND. lscal) THEN
C           ---------------------------------
C           Local (small) scaling arrays have
C           been built, free temporary copies
C           ---------------------------------
            IF (mtype == 1) THEN
              DEALLOCATE(id%COLSCA)
              NULLIFY(id%COLSCA)
            ELSE
              DEALLOCATE(id%ROWSCA)
              NULLIFY(id%ROWSCA)
            ENDIF
            nb_bytes = nb_bytes - int(id%N,8)*k16_8
          ENDIF
        ENDIF ! I_AM_SLAVE
        IF (keep(23) .NE. 0 .AND. mtype==1) THEN
C         Broadcast the unsymmetric permutation and
C         permute the indices in ISOL_loc
          IF (id%MYID.NE.master) THEN
            ALLOCATE(id%UNS_PERM(id%N),stat=allocok)
            IF (allocok > 0) THEN
              info(1)=-13
              info(2)=id%N
              GOTO 40
            ENDIF
          ENDIF
        ENDIF
C
C =====================  ERROR handling and propagation ================
 40     CONTINUE
        CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
        IF (info(1) .LT.0 ) GOTO 90
C ======================================================================
C
        IF (keep(23) .NE. 0 .AND. mtype==1) THEN
          CALL mpi_bcast(id%UNS_PERM(1),id%N,mpi_integer,master,
     &               id%COMM,ierr)
          IF (i_am_slave) THEN
            DO i=1, keep(89)
              id%ISOL_loc(i) = id%UNS_PERM(id%ISOL_loc(i))
            ENDDO
          ENDIF
          IF (id%MYID.NE.master) THEN
            DEALLOCATE(id%UNS_PERM)
            NULLIFY(id%UNS_PERM)
          ENDIF
        ENDIF
      ENDIF ! ICNTL(21)=1
C     --------------------------------------
C     Preparation for distributed solution
C     END
C     --------------------------------------
C     ----------------------------
C     Preparation for reduced RHS
C     ----------------------------
      IF ( ( keep(221) .EQ. 1 ) .OR.
     &     ( keep(221) .EQ. 2 )
     &   ) THEN
C       -- First compute MASTER_ROOT_IN_COMM proc number in
C          COMM_NODES on which is mapped the master of the root.
         IF (keep(46).EQ.1) THEN
             master_root_in_comm=master_root
         ELSE
             master_root_in_comm =master_root+1
         ENDIF
         IF ( id%MYID .EQ. master ) THEN
C            --------------------------------
C            Avoid using LREDRHS when id%NRHS is
C            equal to 1, as was done for RHS
C            --------------------------------
             IF (id%NRHS.EQ.1) THEN
               ld_redrhs = id%KEEP(116)
             ELSE
               ld_redrhs = id%LREDRHS
             ENDIF
         ENDIF
         IF (master.NE.master_root_in_comm) THEN
C        -- Make available LD_REDRHS on MASTER_ROOT_IN_COMM
C           This will then be used to test if a single
C           message can be sent
C           (this is possible if LD_REDRHS=SIZE_SCHUR)
            IF ( id%MYID .EQ. master ) THEN
C            -- send LD_REDRHS to MASTER_ROOT_IN_COMM
C               using COMM communicator
             CALL mpi_send(ld_redrhs,1,mpi_integer,
     &       master_root_in_comm, 0, id%COMM,ierr)
            ELSEIF ( id%MYID.EQ.master_root_in_comm) THEN
C            -- recv LD_REDRHS
             CALL mpi_recv(ld_redrhs,1,mpi_integer,
     &       master, 0, id%COMM,status,ierr)
            ENDIF
C          -- other procs not concerned
         ENDIF
      ENDIF
C
      IF ( keep(248)==1 ) THEN  ! Sparse RHS (A-1 or general sparse)
!        JBEG_RHS - current starting column within A-1 or sparse rhs
!                      set in the loop below and used to obtain the
!                      global index of the column of the sparse RHS
!                      Also used to get index in global permutation.
!                      It also allows to skip empty columns;
        jend_rhs = 0 ! last column in current blockin A-1
C
C       Compute and apply permutations
        IF (do_permute_rhs) THEN
C          Allocate PERM_RHS
           ALLOCATE(perm_rhs(id%NRHS),stat=allocok)
           IF (allocok > 0) THEN
                info(1) = -13
                info(2) = id%NRHS
                GOTO 109
           ENDIF
           nb_bytes = nb_bytes +  int(id%NRHS,8)*k34_8
           nb_bytes_max = max(nb_bytes_max,nb_bytes)
           IF (id%MYID.EQ.master) THEN
C           PERM_RHS is computed on MASTER, it might be modified 
C           in case of interleaving and will thus be distributed 
C           (BCAST) to all slaves only later.
C           Compute PERM_RHS
C           on output: PERM_RHS(k) = i means that i is the kth column
C                                    to be processed
            IF (keep(237).EQ.0) THEN
C               Permute RHS : case of GS (General Sparse) RHS
C               IRHS_SPARSE is of size at least NZ_RHS > 0
C               since all this is skipped when NZ_RHS=0. So
C               accessing IRHS_SPARSE(1) is ok.
                CALL zmumps_permute_rhs_gs(
     &            lp, lpok, prokg, mpg, keep(242),
     &            id%SYM_PERM(1), id%N, id%NRHS,
     &            id%IRHS_PTR(1),  id%NRHS+1, 
     &            id%IRHS_SPARSE(1), id%NZ_RHS,
     &            perm_rhs, ierr)
                 IF (ierr.LT.0) THEN
                   info(1) = -9999
                   info(2) = ierr
                   GOTO 109  ! propagate error
                 ENDIF
            ELSE
C            Case of A-1  : 
C            We compute the permutation of the RHS (sparse matrix)
C                     (to compute all inverse entries)
C            We apply permutation to IRHS_SPARSE ONLY.
C            Note NRHS_NONEMPTY holds the nb of non empty columns 
C            in A-1.
              strat_permam1 = keep(242)
                 CALL zmumps_permute_rhs_am1
     &             (strat_permam1, id%SYM_PERM(1),
     &             id%IRHS_PTR(1), id%NRHS+1,
     &             perm_rhs, id%NRHS,
     &             ierr
     &           )
            ENDIF
           ENDIF
        ENDIF
      ENDIF
C
C     Note that within ZMUMPS_SOL_C, PERM_RHS could be used
C     for A-1 case (with DO_PERMUTE_RHS OR INTERLEAVE_RHS
C     being tested) to get the column index for the
C     original matrix of RHS (column index in A-1)
C     of the permuted columns that have been selected.
C     PERM_RHS is also used in ZMUMPS_GATHER_SOLUTION
C     in case of sparse RHS awith DO_PERMUTE_RHS.
C
C     Allocate PERM_RHS of size 1 if not allocated
      IF (.NOT. allocated(perm_rhs)) THEN
          ALLOCATE(perm_rhs(1),stat=allocok)
          IF (allocok > 0) THEN
            info(1) = -13
            info(2) = 1
            GOTO 109
          ENDIF
          nb_bytes = nb_bytes +  int(size(perm_rhs),8)*k34_8
          nb_bytes_max = max(nb_bytes_max,nb_bytes)
      ENDIF
C     Propagate errors
109   CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
      IF (info(1) .LT.0 ) GOTO 90
c --------------------------
c --------------------------
      IF (id%NSLAVES .EQ. 1) THEN
c     - In case of NS/A-1 we may want to permute RHS
c     - for NS thus is to apply permutation to PIVNUL_LIST
*     - before starting loop of NBRHS
       IF (do_permute_rhs .AND. keep(111).NE.0 ) THEN
C     NOTE:
C        when host not working both master and slaves have
C        in this case the complete list
            WRITE(*,*) id%MYID, ':INTERNAL ERROR 1 : ',
     &                 ' PERMUTE RHS during null space computation ',
     &                 ' not available yet '
            CALL mumps_abort()
       ENDIF                     ! End Permute_RHS
      ELSE
         IF (do_permute_rhs .AND. keep(111).NE.0 ) THEN
            WRITE(*,*) id%MYID, ':INTERNAL ERROR 2 : ',
     &                 ' PERMUTE RHS during null space computation ',
     &                 ' not available yet '
            CALL mumps_abort()
C
C
         ENDIF                  !  End DO_PERMUTE_RHS
         IF (interleave_par.AND. (keep(111).NE.0)) THEN
            WRITE(*,*) id%MYID, ':INTERNAL ERROR 3 : ',
     &        ' INTERLEAVE RHS during null space computation ',
     &        ' not available yet '
            CALL mumps_abort()
         ENDIF
         IF (interleave_par.AND.keep(111).EQ.0) THEN
C     -   A-1 + Interleave:
C     permute RHS on master
           IF (id%MYID.EQ.master) THEN
C            -- PERM_RHS must have been already set or initialized
C            -- it is then modified in next routine
             size_working = id%IPTR_WORKING(id%NPROCS+1)-1
             size_iptr_working = id%NPROCS+1
             CALL zmumps_interleave_rhs_am1(
     &        perm_rhs, id%NRHS,
     &        id%IPTR_WORKING(1), size_iptr_working,
     &        id%WORKING(1), size_working,
     &        id%IRHS_PTR(1),
     &        id%STEP(1), id%SYM_PERM(1), id%N, nbrhs,
     &        id%PROCNODE_STEPS(1), keep(28), id%NSLAVES,
     &        keep(199),
     &        keep(493).NE.0, 
     &        keep(495).NE.0, keep(496), prokg, mpg 
     &        )
           ENDIF               !      End Master
         ENDIF                 !  End A-1 and INTERLEAVE_PAR
C -------------
      ENDIF                  ! End Parallel Case
c --------------------------
c
      IF (do_permute_rhs.AND.(keep(111).EQ.0)) THEN
C     --- Distribute PERM_RHS before loop of RHS
C     --- (with null space option PERM_RHS is not allocated / needed
C          to permute the null column pivot list)
        CALL mpi_bcast(perm_rhs(1),
     &            id%NRHS,
     &            mpi_integer,
     &            master, id%COMM,ierr)
      ENDIF
C   L0-threads to be activated iff KEEP(401)>0 and KEEP(400)>0
          IF (keep(401) .GT. 0) THEN
C          L0-threads was requested for solve phase
C          and will be effective only if KEEP(400) >0 
C          which indicates that L0-threads was 
C          performed during analysis+factorization
           IF ( keep(400) .GT. 0 ) THEN
C           Check if number of threads is consistent with
C           the one used during factorization
C{
            nomp = 1
!$          NOMP=omp_get_max_threads()
              IF (keep(400).NE.nomp) THEN
C               NOMP should be the one from analysis
                id%INFO(1) = -58
                id%INFO(2) = keep(400)
            IF (lpok) WRITE(lp,'(A,A,I5,A,I5)')
     &" FAILURE DETECTED IN SOLVE: #threads for KEEP(401)",
     &" changed from",keep(400)," at analysis to", nomp
              ENDIF
           ENDIF
C}
          ENDIF
      IF (keep(400) .GT. 0) THEN
        CALL zmumps_sol_l0omp_li(keep(400))
      ENDIF
C     ==============================
C     BLOCKING ON the number of RHS
C      We work on  a maximum of NBRHS at a time.
C      the leading dimension of RHS is id%LRHS on master
C      and is set to N on slaves
C     ==============================
C  We may want to allow to have NBRHS that varies
C  this is typically the case when a partitionning of
C  the right hand side is performed and leads to
C  irregular partitions.
C  We only have to be sure that the size of each partition
C  is smaller than NBRHS.
      beg_rhs=1
      DO WHILE (beg_rhs.LE.nrhs_nonempty)
C     {
C       ==========================
C       -- NBRHS     : Original block size
C       -- BEG_RHS   : Column index of the first RHS in the list of
C                      non empty RHS (RHS_LOC) to
C                      be processed during this iteration
C       -- NBRHS_EFF : Effective block size at current iteration
C          In case of sparse RHS (KEEP(248)==1) NBRHS_EFF only refers to
C                  non-empty columns and is used to compute NBCOL_INBLOC
C          -- NBCOL_INBLOC  : the number of columns of sparse RHS needed
C                        to get NBRHS_EFF non empty columns columns of
C                        sparse RHS processed at each step
C
        nbrhs_eff    = min(nrhs_nonempty-beg_rhs+1, nbrhs)
C
C       Sparse RHS
C       Free space and reset pointers if needed
        IF (irhs_sparse_copy_allocated) THEN
            nb_bytes =  nb_bytes -
     &       int(size(irhs_sparse_copy),8)*k34_8
            DEALLOCATE(irhs_sparse_copy)
            irhs_sparse_copy_allocated=.false.
            NULLIFY(irhs_sparse_copy)
        ENDIF
        IF (irhs_ptr_copy_allocated) THEN
            nb_bytes =  nb_bytes -
     &       int(size(irhs_ptr_copy),8)*k34_8
            DEALLOCATE(irhs_ptr_copy)
            irhs_ptr_copy_allocated=.false.
            NULLIFY(irhs_ptr_copy)
        ENDIF
        IF (rhs_sparse_copy_allocated) THEN
            nb_bytes =  nb_bytes -
     &       int(size(rhs_sparse_copy),8)*k35_8
            DEALLOCATE(rhs_sparse_copy)
            rhs_sparse_copy_allocated=.false.
            NULLIFY(rhs_sparse_copy)
        ENDIF
C
C       ===========================================================
C       Set LD_RHS and IBEG for the accesses to id%RHS (in cases
C       id%RHS is accessed). Remark that IBEG might still be
C       overwritten later, in case of general sparse right-hand side
C       and centralized solution to skip empty columns 
C       ===========================================================
        IF (
C           slave procs
     &      ( id%MYID .NE. master )
C       even on master when RHS not allocated
     &     .or.
C         Case of Master working but with distributed sol and
C            ( sparse RHS or null space )
C         -- Allocate not needed on host not working
     &    ( i_am_slave .AND. id%MYID .EQ. master .AND.
     &      icntl21 .NE.0 .AND.
     &      ( keep(248).ne.0 .OR. keep(221).EQ.2
     &          .OR. keep(111).NE.0 )
     &    )
     &     .or.
C         Case of Master and
C          (compute entries of INV(A))
C         Even when I am a master with host not working I
C         am in charge of gathering solution to scale it
C         and to copy it back in the sparse RHS format
     &    ( id%MYID .EQ. master .AND. (keep(237).NE.0) )
C
     &    ) THEN
          ld_rhs = id%N
          ibeg   = 1
        ELSE
          ! (id%MYID .eq. MASTER)
          IF ( associated(id%RHS) ) THEN
C             Leading dimension of RHS on master is id%LRHS
              ld_rhs    = max(id%LRHS, id%N)
          ELSE
C             --- LRHS might not be defined (dont use it)
              ld_rhs    = id%N
          ENDIF
          ibeg      = int(beg_rhs-1,8) * int(ld_rhs,8) + 1_8
        ENDIF
C       JBEG_RHS might also be used in DISTRIBUTED_SOLUTION
C       even when RHS is not sparse on input. In this case,
C       there are no empty columns. (If RHS is sparse JBEG_RHS
C       is overwritten).
        jbeg_rhs = beg_rhs
C       ==========================================
C       Shift empty columns in case of sparse RHS
C       ==========================================
        IF ( (id%MYID.EQ.master) .AND.
     &        keep(248)==1  ) THEN
C         update position of JBEG_RHS on first non-empty
C         column of this block
          jbeg_rhs = jend_rhs + 1
          IF (do_permute_rhs.OR.interleave_par) THEN
            DO WHILE ( id%IRHS_PTR(perm_rhs(jbeg_rhs)) .EQ.
     &          id%IRHS_PTR(perm_rhs(jbeg_rhs)+1) )
C             Empty column
              IF ((keep(237).EQ.0).AND.(icntl21.EQ.0).AND.
     &              (keep(221).NE.1) ) THEN
C                General sparse RHS (NOT A-1) and centralized solution
C                Set to zero part of the
C                solution corresponding to empty columns
                 DO i=1, id%N
                   id%RHS(int(perm_rhs(jbeg_rhs) -1,8)*int(ld_rhs,8)+
     &                    int(i,8)) = zero
                 ENDDO
              ENDIF
              jbeg_rhs = jbeg_rhs +1
            ENDDO
          ELSE
            DO WHILE( id%IRHS_PTR(jbeg_rhs) .EQ.
     &        id%IRHS_PTR(jbeg_rhs+1) )
              IF ((keep(237).EQ.0).AND.(icntl21.EQ.0).AND.
     &          (keep(221).NE.1) ) THEN
C               Case of general sparse RHS (NOT A-1) and
C               centralized solution: set to zero part of
C               the solution corresponding to empty columns
                DO i=1, id%N
                  id%RHS(int(jbeg_rhs -1,8)*int(ld_rhs,8) +
     &                   int(i,8)) = zero
                ENDDO
              ENDIF
              IF (keep(221).EQ.1) THEN
C               Reduced RHS set to ZERO
                DO i = 1, id%SIZE_SCHUR
                  id%REDRHS(int(jbeg_rhs-1,8)*int(ld_redrhs,8) +
     &            int(i,8)) =  zero
                ENDDO
              ENDIF
              jbeg_rhs = jbeg_rhs +1
            ENDDO
          ENDIF                 ! End DO_PERMUTE_RHS.OR.INTERLEAVE_PAR
C         Count nb of RHS columns skipped: useful for
C         * ZMUMPS_DISTRIBUTED_SOLUTION to reset those
C           columns to zero.
C         * in case of reduced right-hand side, to set
C           corresponding entries of RHSCOMP to 0 after
C           forward phase.
          nb_rhsskipped = jbeg_rhs - (jend_rhs + 1)
          IF ((keep(248).EQ.1).AND.(keep(237).EQ.0)
     &         .AND. (icntl21.EQ.0))
     &         THEN
             ! case of general sparse rhs with centralized solution,
             !set IBEG to shifted columns
             ! (after empty columns have been skipped)
             ibeg      = int(jbeg_rhs-1,8) * int(ld_rhs,8) + 1_8
          ENDIF
        ENDIF ! of if (id%MYID.EQ.MASTER) .AND.  KEEP(248)==1
        CALL mpi_bcast( jbeg_rhs, 1, mpi_integer,
     &            master, id%COMM, ierr )
C
C       Shift on REDRHS in reduced RHS functionality
C
        IF (id%MYID.EQ.master .AND. keep(221).NE.0) THEN
C         Initialize IBEG_REDRHS
C         Note that REDRHS always has id%NRHS Colmuns
          ibeg_redrhs= int(jbeg_rhs-1,8)*int(ld_redrhs,8) + 1_8
        ELSE
          ibeg_redrhs=-142424_8  ! Should not be used
        ENDIF
C
C       =====================
C       BEGIN
C       Prepare RHS on master
C
#if defined(V_T)
        CALL vtbegin(perm_scal_ini,ierr)
#endif
        IF (id%MYID .eq. master) THEN
C         ======================
          IF (keep(248)==1) THEN
C         ======================
C
C         Sparse RHS format ( A-1 or sparse input format)
C         is provided as input by the user (IRHS_SPARSE ...)
C         --------------------------------------------------
C         Compute NZ_THIS_BLOCK and NBCOL_INBLOC
C         where
C         NZ_THIS_BLOCK is defined
C         as the number of entries in the next NBRHS_EFF
C         non empty columns (note that since they might be permuted
C         then the following formula is not always valid:
C            NZ_THIS_BLOCK=id%IRHS_PTR(BEG_RHS+NBRHS_EFF)-
C     &                    id%IRHS_PTR(BEG_RHS)
C         anyway NBCOL_INBLOC also need be computed so going through
C         columns one at a time is needed.
C
          nbcol        = 0
          nbcol_inbloc = 0
          nz_this_block = 0
C         With exploit sparsity we skip empty rows up to reaching
C         the first non empty column; then we process a block of
C         maximum size NBRHS_EFF except if we reach another empty
C         column. (We are not sure to have a copy allocated
C         and thus cannot compress on the fly, as done naturally
C         for A-1).
          stop_at_next_empty_col = .false.
          DO i=jbeg_rhs, id%NRHS
            nbcol_inbloc = nbcol_inbloc +1
            IF (do_permute_rhs.OR.interleave_par) THEN
C              PERM_RHS(k) = i means that i is the kth
C                            column to be processed
C              PERM_RHS should also be defined for
C                       empty columns i in A-1 (PERM_RHS(K) = i)
               colsize = id%IRHS_PTR(perm_rhs(i)+1)
     &                    - id%IRHS_PTR(perm_rhs(i))
            ELSE
               colsize = id%IRHS_PTR(i+1) - id%IRHS_PTR(i)
            ENDIF
            IF ((.NOT.stop_at_next_empty_col).AND.(colsize.GT.0).AND.
     &          (keep(237).EQ.0)) THEN
C              -- set STOP_NEXT_EMPTY_COL only for general
C              --  sparse case (not AM-1)
               stop_at_next_empty_col =.true.
            ENDIF
            IF (colsize.GT.0
     &      ) THEN
              nbcol = nbcol+1
              nz_this_block = nz_this_block + colsize
            ELSE IF (stop_at_next_empty_col) THEN
C We have reached an empty column with already selected non empty 
C columns: reduce block size to non empty columns reached so far.
              nbcol_inbloc = nbcol_inbloc -1
              nbrhs_eff = nbcol
              EXIT
            ENDIF
            IF (nbcol.EQ.nbrhs_eff) EXIT
          ENDDO
          IF (nz_this_block.EQ.0) THEN
           WRITE(*,*) " Internal Error 16 in sol driver NZ_THIS_BLOCK=",
     &               nz_this_block
           CALL mumps_abort()
          ENDIF
C
          IF (nbcol.NE.nbrhs_eff.AND. (keep(237).NE.0)
     &         .AND.keep(221).NE.1) THEN
C           With exploit sparsity for general sparse RHS (Not A-1)
C           we skip empty rows up to reaching
C           the first non empty column; then we process a block of
C           maximum size NBRHS_EFF except if we reach another empty
C           column. (We are not sure to have a copy allocated
C           and thus cannot compress on the fly, as done naturally
C           for A-1). Thus NBCOL might be smaller than NBRHS_EFF
            WRITE(6,*) ' Internal Error 8 in solution driver ',
     &            nbcol, nbrhs_eff
            call mumps_abort()
          ENDIF
C         -------------------------------------------------------------
C
          IF (nz_this_block .NE. 0) THEN
C           -----------------------------------------------------------
C           We recall that
C           NBCOL_INBLOC is the number of columns of sparse RHS needed 
C           to get NBRHS_EFF non empty columns:
            ALLOCATE(irhs_ptr_copy(nbcol_inbloc+1),stat=allocok)
            if (allocok .GT.0 ) then
              info(1)=-13
              info(2)=nbcol_inbloc+1
              GOTO 30
            endif
            irhs_ptr_copy_allocated = .true.
            nb_bytes =  nb_bytes +  int(nbcol_inbloc+1,8)*k34_8
            nb_bytes_max = max(nb_bytes_max,nb_bytes)
C
            jend_rhs =jbeg_rhs + nbcol_inbloc - 1
C           -----------------------------------------------------------
C           Initialize IRHS_PTR_COPY
C           compute local copy (compressed) of id%IRHS_PTR on Master
            IF (do_permute_rhs.OR.interleave_par) THEN
              ipos = 1
              j = 0
              DO i=jbeg_rhs, jbeg_rhs + nbcol_inbloc -1
                j = j+1
                irhs_ptr_copy(j) = ipos
                colsize = id%IRHS_PTR(perm_rhs(i)+1)
     &                    - id%IRHS_PTR(perm_rhs(i))
                ipos = ipos + colsize
              ENDDO
            ELSE
              ipos = 1
              j = 0
              DO i=jbeg_rhs, jbeg_rhs + nbcol_inbloc -1
                j = j+1
                irhs_ptr_copy(j) = ipos
                colsize = id%IRHS_PTR(i+1)
     &                     - id%IRHS_PTR(i)
                ipos = ipos + colsize
              ENDDO
            ENDIF                 ! End DO_PERMUTE_RHS.OR.INTERLEAVE_PAR
            irhs_ptr_copy(nbcol_inbloc+1)= ipos
            IF ( ipos-1 .NE. nz_this_block ) THEN
                WRITE(*,*) "Error in compressed copy of IRHS_PTR"
                ierr = 99
                call mumps_abort()
            ENDIF
C           -----------------------------------------------------------
C           IRHS_SPARSE : do a copy or point to the original indices
C
C           Check whether IRHS_SPARSE_COPY need be allocated
            IF (keep(23) .NE. 0 .and. mtype .NE. 1) THEN
C             AP = LU and At x = b ==> b need be permuted
              ALLOCATE(irhs_sparse_copy(nz_this_block)
     &               ,stat=allocok)
              if (allocok .GT.0 ) then
                info(1)=-13
                info(2)=nz_this_block
                GOTO 30
              endif
              irhs_sparse_copy_allocated=.true.
              nb_bytes = nb_bytes + int(nz_this_block,8)*k34_8
              nb_bytes_max = max(nb_bytes_max,nb_bytes)
            ELSE IF (do_permute_rhs.OR.interleave_par.OR.
     &           (keep(237).NE.0)) THEN
C               Columns are not contiguous and need be copied one by one
C               IRHS_SPARSE_COPY will hold a copy of contiguous permuted
C               columns so an explicit copy is needed.
C               IRHS_SPARSE_COPY is also allways allocated with A-1,
C               to enable receiving during mumps_gather_solution
C     .         on the master in any order.
                ALLOCATE(irhs_sparse_copy(nz_this_block),
     &                   stat=allocok)
                IF (allocok .GT.0 ) THEN
                    ierr = 99
                    GOTO 30
                ENDIF
                irhs_sparse_copy_allocated=.true.
                nb_bytes = nb_bytes + int(nz_this_block,8)*k34_8
                nb_bytes_max = max(nb_bytes_max,nb_bytes)
C
            ENDIF
C          
C     Initialize IRHS_SPARSE_COPY
            IF (irhs_sparse_copy_allocated) THEN
                IF ( do_permute_rhs.OR.interleave_par ) THEN
                  ipos = 1
                  DO i=jbeg_rhs, jbeg_rhs + nbcol_inbloc -1
                    colsize = id%IRHS_PTR(perm_rhs(i)+1)
     &                       - id%IRHS_PTR(perm_rhs(i))
                    irhs_sparse_copy(ipos:ipos+colsize-1) =
     &              id%IRHS_SPARSE(id%IRHS_PTR(perm_rhs(i)):
     &              id%IRHS_PTR(perm_rhs(i)+1) -1)
                    ipos = ipos + colsize
                  ENDDO
                ELSE
                  irhs_sparse_copy = id%IRHS_SPARSE(
     &              id%IRHS_PTR(jbeg_rhs):
     &              id%IRHS_PTR(jbeg_rhs)+nz_this_block-1)
                ENDIF
            ELSE
                irhs_sparse_copy
c    *                   (1:NZ_THIS_BLOCK)
     &           =>
     &          id%IRHS_SPARSE(id%IRHS_PTR(jbeg_rhs):
     &             id%IRHS_PTR(jbeg_rhs)+nz_this_block-1)
            ENDIF
            IF (lscal.OR.do_permute_rhs.OR.interleave_par.OR.
     &          (keep(237).NE.0)) THEN
C             if scaling is on or if columns of the RHS are
C             permuted then a copy of RHS_SPARSE is needed.
C             Also always allocated with A-1,
c             to enable receiving during mumps_gather_solution
C             on the master in any order.
C
              ALLOCATE(rhs_sparse_copy(nz_this_block),
     &               stat=allocok)
              IF (allocok .GT.0 ) THEN
                info(1)=-13
                info(2)=nz_this_block
                GOTO 30
              ENDIF
              rhs_sparse_copy_allocated = .true.
              nb_bytes = nb_bytes + int(nz_this_block,8)*k35_8
              nb_bytes_max = max(nb_bytes_max,nb_bytes)
            ELSE
              IF  ( keep(248)==1 ) THEN
                rhs_sparse_copy
c    *            (1:NZ_THIS_BLOCK)
     &          => id%RHS_SPARSE(id%IRHS_PTR(jbeg_rhs):
     &             id%IRHS_PTR(jbeg_rhs)+nz_this_block-1)
              ELSE
                rhs_sparse_copy
c    *            (1:NZ_THIS_BLOCK)
     &          => id%RHS_SPARSE(id%IRHS_PTR(beg_rhs):
     &             id%IRHS_PTR(beg_rhs)+nz_this_block-1)
              ENDIF
            ENDIF
            IF (do_permute_rhs.OR.interleave_par.OR.
     &          (id%KEEP(237).NE.0)) THEN
              IF (id%KEEP(237).NE.0) THEN
C               --initialized to one; it might be
C               modified if scaling is on (one first entry in each col is scaled)
                rhs_sparse_copy = one
              ELSE IF (.NOT. lscal) THEN
C               -- Columns are not contiguous and need be copied one by one
C               -- This need not be done if scaling is on because it
C               -- will done and scaled later.
                ipos = 1
                DO i=jbeg_rhs, jbeg_rhs + nbcol_inbloc -1
                  colsize = id%IRHS_PTR(perm_rhs(i)+1)
     &                    - id%IRHS_PTR(perm_rhs(i))
                  IF (colsize .EQ. 0) cycle
                  rhs_sparse_copy(ipos:ipos+colsize-1) =
     &            id%RHS_SPARSE(id%IRHS_PTR(perm_rhs(i)):
     &            id%IRHS_PTR(perm_rhs(i)+1) -1)
                  ipos = ipos + colsize
                ENDDO
              ENDIF
            ENDIF
C           =========================
            IF (keep(23) .NE. 0) THEN
C           =========================
*             maximum transversal was performed
              IF (mtype .NE. 1) THEN
*               At x = b is asked while
*               we have AP = LU   where P is the column permutation
*               due to max trans.
*               Therefore we need to modify rhs:
*                 b' = P-1 b   (P-1=Pt)
*               Apply column permutation to the right hand side RHS
*               Column J of the permuted matrix corresponds to
*               column PERMW(J) of the original matrix.
*
C               ==========
C               SPARSE RHS : permute indices rather than values
C               ==========
C               Solve with At X = B should never occur for A-1 
                ipos = 1
                DO i=1, nbcol_inbloc
C                 Note that: (i) IRHS_PTR_COPY is compressed;
C                 (ii) columns might have been permuted
                  colsize = irhs_ptr_copy(i+1) - irhs_ptr_copy(i)
                  DO k = 1, colsize
                   jperm = uns_perm_inv(irhs_sparse_copy(ipos+k-1))
                   irhs_sparse_copy(ipos+k-1) = jperm
                  ENDDO
                  ipos = ipos + colsize
                ENDDO
              ENDIF ! MTYPE.NE.1
            ENDIF ! KEEP(23).NE.0
          ENDIF ! NZ_THIS_BLOCK .NE. 0
C         -----
          ENDIF  !  ============ KEEP(248)==1
C         -----
        ENDIF   !  (id%MYID .eq. MASTER)
C
C =====================  ERROR handling and propagation ================
 30     CONTINUE
        CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
        IF (info(1) .LT.0 ) GOTO 90
C ======================================================================
C
C       NBCOL_INBLOC depends on loop
        IF (keep(248)==1) THEN
            CALL mpi_bcast( nbcol_inbloc,1, mpi_integer,
     &           master, id%COMM,ierr)
        ELSE
          nbcol_inbloc = nbrhs_eff
        ENDIF
        jend_rhs =jbeg_rhs + nbcol_inbloc - 1
        IF ((keep(111).eq.0).AND.(keep(252).EQ.0)
     &      .AND.(keep(221).NE.2 ).AND.(keep(248).EQ.1) ) THEN
C         ----------------------------
C         -- SPARSE RIGHT-HAND-SIDE
C         ----------------------------
          CALL mpi_bcast( nz_this_block,1, mpi_integer,
     &                    master, id%COMM,ierr)
          IF (id%MYID.NE.master .and. nz_this_block.NE.0) THEN
            ALLOCATE(irhs_sparse_copy(nz_this_block),
     &               stat=allocok)
            if (allocok .GT.0 ) then
               info(1)=-13
               info(2)=nz_this_block
               GOTO 45
            endif
            irhs_sparse_copy_allocated=.true.
C           RHS_SPARSE_COPY is broadcasted
C           for A-1 even if on the slaves the initialisation of the RHS
C           could be only based on the pattern. Doing so we
C           broadcast the scaled version of the RHS (scaling arrays
C           that are not available on slaves).
            ALLOCATE(rhs_sparse_copy(nz_this_block),
     &               stat=allocok)
            if (allocok .GT.0 ) then
               info(1)=-13
               info(2)=nz_this_block
               GOTO 45
            endif
            rhs_sparse_copy_allocated=.true.
            nb_bytes = nb_bytes + int(nz_this_block,8)*(k34_8+k35_8)
            nb_bytes_max = max(nb_bytes_max,nb_bytes)
C
            ALLOCATE(irhs_ptr_copy(nbcol_inbloc+1),stat=allocok)
            if (allocok .GT.0 ) then
               info(1)=-13
               info(2)=nbcol_inbloc+1
               GOTO 45
            endif
            irhs_ptr_copy_allocated = .true.
            nb_bytes = nb_bytes + int(nbcol_inbloc+1,8)*k34_8
            nb_bytes_max = max(nb_bytes_max,nb_bytes)
          ENDIF
C
C =====================  ERROR handling and propagation ================
 45       CONTINUE
          CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
          IF (info(1) .LT.0 ) GOTO 90
C ======================================================================
          IF (nz_this_block > 0) THEN
            CALL mpi_bcast(irhs_sparse_copy(1),
     &                 nz_this_block,
     &                 mpi_integer,
     &                 master, id%COMM,ierr)
            CALL mpi_bcast(irhs_ptr_copy(1),
     &             nbcol_inbloc+1,
     &             mpi_integer,
     &             master, id%COMM,ierr)
            IF (ierr.GT.0) THEN
               WRITE (*,*)'NOT OK FOR ALLOC PTR ON SLAVES'
               call mumps_abort()
            ENDIF
          ENDIF
        ENDIF
C
C       =========================================================
C       INITIALIZE POSINRHSCOMP_ROW/COL, RHSCOMP and related data
C       For distributed RHS, initialize RHSMAPINFO (at 1st block)
C       =========================================================
        IF ( i_am_slave ) THEN
C         --------------------------------------------------
C         If I am involved in the solve and if
C         either
C           no null space comput (keep(111)=0) and sparse rhs
C         or
C             null space computation
C         then
C           compute POSINRHSCOMP
C         endif
C
C         Fwd in facto: in this case only POSINRHSCOMP need be computed
C
C         (POSINRHSCOMP_ROW/COL indirection arrays should
C         have been allocated once outside loop)
C         Compute size of RHSCOMP since it might depend
C            on the process index and of the sparsity of the RHS
C            if it is exploited.
C         Initialize POSINRHSCOMP_ROW/COL
C
C         Note that LD_RHSCOMP and id%KEEP8(25)
C         are not set on the host in this routine in
C         the case of a non-working host.
C         Note that POSINRHSCOMP is now always computed in SOL_DRIVER
C         at least during the first block of RHS when sparsity of RHS
C         is not exploited.
C         -------------------------------
C         INITTIALZE POSINRHSCOMP_ROW/COL
C         -------------------------------
C
          IF ( keep(221).EQ.2 .AND. keep(252).EQ.0
     &      .AND.  (keep(248).NE.1 .OR. (id%NRHS.EQ.1))
     &      ) THEN
C           Reduced RHS was already computed during
C           a previous forward step AND is valid.
C           By valid we mean:
C            -no forward in facto  (KEEP(252)==0) during which
C             POSINRHSCOMP was not computed
C           AND
C            -no exploit sparsity with multiple RHS
C            because in this case POSINRHSCOMP would
C            be valid only for the last block processed during fwd.
C           In those cases since we only perform the backward step, we do not
C           need to compute POSINRHSCOMP
            build_posinrhscomp = .false.
          ENDIF
C         ------------------------
C         INITIALIZE POSINRHSCOMP
C         ------------------------
          IF (build_posinrhscomp) THEN
C           -- we first set MTYPE_LOC and
C           -- reset BUILD_POSINRHSCOMP for next iteration in loop
C
C           general case only POSINRHSCOMP is computed
            build_posinrhscomp = .false.
!           POSINRHSCOMP does not change between blocks
            mtype_loc = mtype
C
            IF ( (keep(111).NE.0) .OR. (keep(237).NE.0) .OR.
     &         (keep(252).NE.0) ) THEN
C
              IF (keep(111).NE.0) THEN
C               -- in the context of null space, we need to
C               -- build RHSCOMP to skip SOL_R. Therefore
C               -- we need to know for each concerned
C               -- row index its position in
C               -- RHSCOMP
C               We use row indices, as these are the ones that
C               were used to detect zero pivots during factorization.
C               POSINRHSCOMP_ROW will allow to find the (row) index of a
C               zero in RHSCOMP before calling ZMUMPS_SOL_S. Then
C               ZMUMPS_SOL_S uses column indices to build the solution
C               (corresponding to null space vectors)
                mtype_loc = 1
              ELSE IF  (keep(252).NE.0) THEN
C               -- Fwd in facto: since fwd is skipped we need to build POSINRHSCOMP
                mtype_loc = 1  ! (no transpose)
C     BUILD_POSINRHSCOMP = .FALSE.  ! POSINRHSCOMP does not change between blocks
              ELSE
C               -- A-1 only
                mtype_loc = mtype
                build_posinrhscomp = .true.
              ENDIF
            ENDIF
C           -- compute POSINRHSCOMP
            liw_passed=max(1,liw)
            IF (keep(237).EQ.0) THEN
              CALL zmumps_build_posinrhscomp(
     &           id%NSLAVES,id%N,
     &           id%MYID_NODES, id%PTLUST_S(1),
     &           id%KEEP(1),id%KEEP8(1),
     &           id%PROCNODE_STEPS(1), id%IS(1), liw_passed,
     &           id%STEP(1),
     &           id%POSINRHSCOMP_ROW(1), id%POSINRHSCOMP_COL(1),
     &           id%POSINRHSCOMP_COL_ALLOC,
     &           mtype_loc,
     &           nbent_rhscomp, nb_fs_rhscomp_tot )
                nb_fs_rhscomp_f = nb_fs_rhscomp_tot
            ELSE
              CALL zmumps_build_posinrhscomp_am1(
     &           id%NSLAVES,id%N,
     &           id%MYID_NODES, id%PTLUST_S(1), id%DAD_STEPS(1),
     &           id%KEEP(1),id%KEEP8(1),
     &           id%PROCNODE_STEPS(1), id%IS(1), liw,
     &           id%STEP(1),
     &           id%POSINRHSCOMP_ROW(1), id%POSINRHSCOMP_COL(1),
     &           id%POSINRHSCOMP_COL_ALLOC,
     &           mtype_loc,
     &           irhs_ptr_copy(1), nbcol_inbloc, irhs_sparse_copy(1),
     &           nz_this_block,perm_rhs, size(perm_rhs) , jbeg_rhs,
     &           nbent_rhscomp,
     &           nb_fs_rhscomp_f, nb_fs_rhscomp_tot,
     &            uns_perm_inv, size(uns_perm_inv) ! size 1 if not used
     &            )
            ENDIF
          ENDIF   ! BUILD_POSINRHSCOMP=.TRUE.
          IF (build_rhsmapinfo .AND. keep(248).EQ.-1) THEN
C
C           Prepare symbolic data for sends.
C           For the moment: MAP_RHS_loc
C
            CALL mumps_sol_rhsmapinfo( id%N, id%Nloc_RHS, id%KEEP(89),
     &            irhs_loc_ptr(1), map_rhs_loc, id%POSINRHSCOMP_ROW(1),
     &            id%NSLAVES, id%MYID_NODES,
     &            id%COMM_NODES, id%ICNTL(1), id%INFO(1) )
            build_rhsmapinfo = .false.
C           MUMPS_SOL_RHSMAPINFO does not propagate errors
          ENDIF
        ENDIF
        CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
        IF (info(1) .LT.0 ) GOTO 90
        IF (i_am_slave) THEN
          IF (keep(221).EQ.1) THEN
C           we need to save the reduced RHS for all RHS to perform
C           later the backward phase with an updated reduced RHS
C           thus we allocate NRHS_NONEMPTY columns in one shot.
C           Note that RHSCOMP might have been allocated in previous block
C           and RHSCOMP has been deallocated previous to entering loop on RHS
            IF (.not. associated(id%RHSCOMP)) THEN
C             So far we cannot combine this to exploit sparsity
C             so that NBENT_RHSCOMP will not change in the loop
C             and can be used to dimension RHSCOMP
C             Furthermore, during bwd phase the REDRHS provided
C             by the user might also have a different non empty
C             column pattern than the sparse RHS provided on input to
C             this phase: thus we need to allocate id%NRHS columns too.
              ld_rhscomp = max(nbent_rhscomp,1)
              id%KEEP8(25) = int(ld_rhscomp,8)*int(id%NRHS,8)
              ALLOCATE (id%RHSCOMP(id%KEEP8(25)),  stat = allocok)
              IF ( allocok .GT. 0 ) THEN
                 info(1)=-13
                 CALL mumps_set_ierror(id%KEEP8(25),info(2))
                 id%KEEP8(25)=0_8
                 GOTO 41
              END IF
              nb_bytes = nb_bytes + id%KEEP8(25)*k35_8
              nb_bytes_max = max(nb_bytes_max,nb_bytes)
            ENDIF
          ENDIF
          IF ((keep(221).NE.1).AND.
     &        ((keep(221).NE.2).OR.(keep(252).NE.0))
     &       ) THEN
C           ------------------
C           Allocate RHSCOMP (case of RHSCOMP allocated at each block of RHS)
C           ------------------
C           RHSCOMP allocated per block of maximum size NBRHS
            ld_rhscomp = max(nbent_rhscomp, ld_rhscomp)
C           NBRHS_EFF could be used instead on NBRHS
            IF (associated(id%RHSCOMP)) THEN
              IF ( (id%KEEP8(25).LT.int(ld_rhscomp,8)*int(nbrhs,8))
     &          .OR. (keep(235).NE.0).OR.(keep(237).NE.0) ) THEN
                ! deallocate and reallocate if:
                ! _larger array needed
                ! OR 
                ! _exploit sparsity/A-1: since size of RHSCOMP 
                ! is expected to vary much in these cases
                ! this should improve locality
                 nb_bytes = nb_bytes - id%KEEP8(25)*k35_8
                 DEALLOCATE(id%RHSCOMP)
                 NULLIFY(id%RHSCOMP)
                 id%KEEP8(25)=0_8
              ENDIF
            ENDIF
            IF (.not. associated(id%RHSCOMP)) THEN
              ld_rhscomp = max(nbent_rhscomp, 1)
              id%KEEP8(25) = int(ld_rhscomp,8)*int(nbrhs,8)
              ALLOCATE (id%RHSCOMP(id%KEEP8(25)),  stat = allocok )
              IF ( allocok .GT. 0 ) THEN
               info(1)=-13
               CALL mumps_set_ierror(id%KEEP8(25),info(2))
               GOTO 41
              END IF
              nb_bytes = nb_bytes + id%KEEP8(25)*k35_8
              nb_bytes_max = max(nb_bytes_max,nb_bytes)
            ENDIF
          ENDIF
          IF (keep(221).EQ.2) THEN
C           RHSCOMP has been allocated (call with KEEP(221).EQ.1)
C           even in the case fwd in facto
            ! Not correct: LD_RHSCOMP = LENRHSCOMP/id%NRHS_NONEMPTY
            ld_rhscomp = int(id%KEEP8(25)/int(id%NRHS,8))
          ENDIF
C
C         Shift on RHSCOMP
C
          IF ( keep(221).EQ.0 ) THEN
C            -- RHSCOMP reused in the loop
             ibeg_rhscomp= 1_8
          ELSE
C            Initialize IBEG_RHSCOMP
C
             ibeg_rhscomp= int(jbeg_rhs-1,8)*int(ld_rhscomp,8) + 1_8
          ENDIF
        ENDIF   ! I_AM_SLAVE
C =====================  ERROR handling and propagation ================
 41     CONTINUE
        CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
        IF (info(1) .LT.0 ) GOTO 90
C ======================================================================
C
C       ---------------------------
C       Prepare RHS on master (case
C       of dense and sparse RHS)
C       ---------------------------
        IF (id%MYID .eq. master) THEN
C         =========================
          IF (keep(23) .NE. 0) THEN
C         =========================
*         maximum transversal was performed
            IF (mtype .NE. 1) THEN
*             At x = b is asked while
*             we have AP = LU   where P is the column permutation
*             due to max trans.
*             Therefore we need to modify rhs:
*               b' = P-1 b   (P-1=Pt)
*             Apply column permutation to the right hand side RHS
*             Column J of the permuted matrix corresponds to
*             column PERMW(J) of the original matrix.
*
              IF (keep(248)==0) THEN
C               =========
C               DENSE RHS : permute values in RHS
C               =========
                ALLOCATE( c_rw2( id%N ),stat =allocok )
                IF ( allocok .GT. 0 ) THEN
                  info(1)=-13
                  info(2)=id%N
                  IF (lpok) THEN
                    WRITE(lp,*) id%MYID,
     &              ':Error allocating C_RW2 in ZMUMPS_SOLVE_DRIVE'
                  END IF
                  GOTO 30
                END IF
C               We directly permute in id%RHS.
                DO k = 1, nbrhs_eff
                  kdec = ibeg+int(k-1,8)*int(ld_rhs,8)
                  DO i = 1, id%N
                    c_rw2(i)=id%RHS(i-1+kdec)
                  END DO
                  DO i = 1, id%N
                    jperm = id%UNS_PERM(i)
                    id%RHS(i-1+kdec) = c_rw2(jperm)
                  END DO
                END DO
                DEALLOCATE(c_rw2)
              ENDIF
            ENDIF
          ENDIF
C
          IF (postpros) THEN
            IF ( keep(248) == 0 ) THEN
              DO k = 1, nbrhs_eff
                kdec = ibeg+int(k-1,8)*int(ld_rhs,8)
                DO i = 1, id%N
                  saverhs(i+(k-1)*id%N) = id%RHS(kdec+i-1)
                END DO
              ENDDO
            ELSE IF (keep(248)==1) THEN
              saverhs(:) = zero
              DO k = 1, nbrhs
                DO j = id%IRHS_PTR(k), id%IRHS_PTR(k+1)-1
                  i = id%IRHS_SPARSE(j)
                  saverhs(i+(k-1)*id%N) = id%RHS_SPARSE(j)
                ENDDO
              ENDDO
            ENDIF
          ENDIF
C
C         RHS is set to scaled right hand side
C
          IF (lscal) THEN
C           scaling was performed
            IF (keep(248)==0) THEN
C             dense RHS
              IF (mtype .EQ. 1) THEN
C               we solve Ax=b, use ROWSCA to scale the RHS
                DO k =1, nbrhs_eff
                  kdec = int(k-1,8) * int(ld_rhs,8) + int(ibeg-1,8)
                  DO i = 1, id%N
                    id%RHS(kdec+i) = id%RHS(kdec+i) *
     &                               id%ROWSCA(i)
                  ENDDO
                ENDDO
              ELSE
C               we solve Atx=b, use COLSCA to scale the RHS
                DO k =1, nbrhs_eff
                  kdec = int(k-1,8) * int(ld_rhs,8) + int(ibeg-1,8)
                  DO i = 1, id%N
                   id%RHS(kdec+i) = id%RHS(kdec+i) *
     &                              id%COLSCA(i)
                  ENDDO
                ENDDO
              ENDIF
            ELSE IF (keep(248)==1) THEN
C             -------------------------
C             KEEP(248)==1 (and MASTER)
C             -------------------------
              kdec=int(id%IRHS_PTR(jbeg_rhs),8)
C             Compute
              IF ((keep(248)==1) .AND.
     &          (do_permute_rhs.OR.interleave_par.OR.
     &           (id%KEEP(237).NE.0))
     &        ) THEN
C               -- copy from RHS_SPARSE need be done per
C                  column following PERM_RHS
C               Columns are not contiguous and need be copied one by one
                ipos = 1
                j    = 0
                DO i=jbeg_rhs, jbeg_rhs + nbcol_inbloc -1
                  IF (do_permute_rhs.OR.interleave_par) THEN
                      iperm = perm_rhs(i)
                  ENDIF 
                  j = j+1
C                 Note that we work here on compressed IRHS_PTR_COPY
                  colsize = irhs_ptr_copy(j+1) - irhs_ptr_copy(j)
C                 -- skip empty column
                  IF (colsize .EQ. 0) cycle
                  IF (id%KEEP(237).NE.0) THEN
                    IF (do_permute_rhs.OR.interleave_par) THEN
C                     if A-1 only, then, for each non empty target
C                     column PERM_RHS(I), scale in first position
C                     in column the diagonal entry
C                     build the scaled rhs ej on each slave.
                      rhs_sparse_copy(ipos) =   id%ROWSCA(iperm) *
     &                       one
                    ELSE
                      rhs_sparse_copy(ipos) =   id%ROWSCA(i) * one
                    ENDIF
                  ELSE
C                   Loop over nonzeros in column
                    DO k = 1, colsize
C                     Formula for II below is ok, except in case
C                     of maximum transversal (KEEP(23).NE.0) and
C                     transpose system (MTYPE .NE. 1):
C                     II = id%IRHS_SPARSE(id%IRHS_PTR(PERM_RHS(I))+K-1)
C                     In case of maximum transversal + transpose, one
C                     should then apply II=UNS_PERM_INV(II) after the
C                     above definition of II.
C
C                     Instead, we rely on IRHS_SPARSE_COPY, whose row
C                     indices have already been permuted in case of
C                     maximum transversal.
                      ii = irhs_sparse_copy(
     &                          irhs_ptr_copy(i-jbeg_rhs+1) 
     &                     +k-1)
C                     PERM_RHS(I) corresponds to column in original RHS.
C                     Original IRHS_PTR must be used to access id%RHS_SPARSE
                      IF (mtype.EQ.1) THEN
                        rhs_sparse_copy(ipos+k-1) =
     &                  id%RHS_SPARSE(id%IRHS_PTR(iperm)+k-1)*
     &                  id%ROWSCA(ii)
                      ELSE
                          rhs_sparse_copy(ipos+k-1) =
     &                    id%RHS_SPARSE(id%IRHS_PTR(iperm)+k-1)*
     &                    id%COLSCA(ii)
                      ENDIF
                    ENDDO
                  ENDIF
                  ipos = ipos + colsize
                ENDDO
              ELSE
                ! general sparse RHS
                ! without permutation
                IF (mtype .eq. 1) THEN
                  DO iz=1,nz_this_block
                    i=irhs_sparse_copy(iz)
                    rhs_sparse_copy(iz)=id%RHS_SPARSE(kdec+iz-1)*
     &                            id%ROWSCA(i)
                  ENDDO
                ELSE
                  DO iz=1,nz_this_block
                    i=irhs_sparse_copy(iz)
                    rhs_sparse_copy(iz)=id%RHS_SPARSE(kdec+iz-1)*
     &                                  id%COLSCA(i)
                  ENDDO
                ENDIF
              ENDIF
            ENDIF  ! KEEP(248)==1
          ENDIF  ! LSCAL
        ENDIF  ! id%MYID.EQ.MASTER
#if defined(v_t)
        CALL vtend(perm_scal_ini,ierr)
#endif
C
C       Prepare RHS on master
C       END
C       =====================
        IF ((keep(248).EQ.1).AND.(keep(237).EQ.0)) THEN
          ! case of general sparse: in case of empty columns
          ! modifed version of
          ! NBRHS_EFF need be broadcasted since it is used
          ! to update BEG_RHS at the end of the DO WHILE
            CALL mpi_bcast( nbrhs_eff,1, mpi_integer,
     &           master, id%COMM,ierr)
            CALL mpi_bcast(nb_rhsskipped,1,mpi_integer,master,
     &               id%COMM,ierr)
        ENDIF
C       -----------------------------------
C       Two main cases depending on option
C       for null space computation:
C
C       KEEP(111)=0 : use RHS from user
C                     (sparse or dense)
C       KEEP(111)!=0: build an RHS on each
C                     proc for null space
C                     computations
C       -----------------------------------
#if defined(V_T)
        CALL vtbegin(soln_dist,ierr)
#endif
        timescatter1=mpi_wtime()
        IF ((keep(111).eq.0).AND.(keep(252).EQ.0)
     &        .AND.(keep(221).NE.2 )) THEN
C           ------------------------
C           Use RHS provided by user
C           when not null space and not Fwd in facto
C           ------------------------
            IF (keep(248) == 0) THEN
C             ----------------------------
C             -- DENSE RIGHT-HAND-SIDE
C             ----------------------------
              IF ( .NOT.i_am_slave ) THEN
C               -- Master not working
                CALL zmumps_scatter_rhs(id%NSLAVES,id%N, id%MYID,
     &          id%COMM,
     &          mtype, id%RHS(ibeg), ld_rhs, nbrhs_eff,
     &          nbrhs_eff,
     &          c_dummy, 1, 1,
     &          idummy, 0,
     &          jdummy, id%KEEP(1), id%KEEP8(1), id%PROCNODE_STEPS(1),
     &          idummy, 1,
     &          id%STEP(1),
     &          id%ICNTL(1),id%INFO(1))
              ELSE
                IF (id%MYID .eq. master) THEN
                  ptr_rhs => id%RHS
                  ld_rhs_loc   = ld_rhs
                  ncol_rhs_loc = nbrhs_eff
                  ibeg_loc     = ibeg
                ELSE
                  ptr_rhs => cdummy_target
                  ld_rhs_loc     = 1
                  ncol_rhs_loc   = 1
                  ibeg_loc       = 1_8
                ENDIF
                liw_passed = max( liw, 1 )
                CALL zmumps_scatter_rhs(id%NSLAVES,id%N, id%MYID,
     &          id%COMM,
     &          mtype, ptr_rhs(ibeg_loc),ld_rhs_loc,ncol_rhs_loc,
     &          nbrhs_eff,
     &          id%RHSCOMP(ibeg_rhscomp), ld_rhscomp, nbrhs_eff,
     &          id%POSINRHSCOMP_ROW(1), nb_fs_rhscomp_f,
C
     &          id%PTLUST_S(1), id%KEEP(1), id%KEEP8(1),
     &          id%PROCNODE_STEPS(1),
     &          is(1), liw_passed,
     &          id%STEP(1),
     &          id%ICNTL(1),id%INFO(1))
              ENDIF
              IF (info(1).LT.0) GOTO 90
            ELSE IF (keep(248) .EQ. -1) THEN
              IF (i_am_slave) THEN
                IF (id%Nloc_RHS .NE. 0) THEN
                  rhs_loc_size=int(id%LRHS_loc,8)*int(nbrhs_eff-1,8)+
     &                         int(id%Nloc_RHS,8)
                  rhs_loc_shift=1_8+int(beg_rhs-1,8)*id%LRHS_loc
                ELSE
                  rhs_loc_size=1_8
                  rhs_loc_shift=1_8
                ENDIF
                CALL zmumps_scatter_dist_rhs(id%NSLAVES, id%N,
     &          id%MYID_NODES, id%COMM_NODES,
     &          nbrhs_eff, id%Nloc_RHS, id%LRHS_loc,
     &          map_rhs_loc,
     &          irhs_loc_ptr(1),
     &          idrhs_loc(rhs_loc_shift),
     &          rhs_loc_size,
     &          id%RHSCOMP(ibeg_rhscomp), ld_rhscomp,
     &          id%POSINRHSCOMP_ROW(1), nb_fs_rhscomp_f,
     &          lscal, scaling_data_dr,
     &          lp, lpok, keep(1), nb_bytes_loc, info(1))
C               NB_BYTES_LOC were allocated and freed above
                nb_bytes_max = max(nb_bytes_max,
     &                             nb_bytes_max+nb_bytes_loc)
              ENDIF
              CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
              IF (info(1).LT.0) GOTO 90
            ELSE
C             ===  KEEP(248)==1 =========
C             -- SPARSE RIGHT-HAND-SIDE
C             ----------------------------
              IF (nz_this_block > 0) THEN
                CALL mpi_bcast(rhs_sparse_copy(1),
     &                     nz_this_block,
     &                     mpi_double_complex,
     &                     master, id%COMM, ierr)
              ENDIF
C             -- At this point each process has a copy of the
C             -- sparse RHS. We need to store it into RHSCOMP.
C
              IF (keep(237).NE.0)  THEN
                IF ( i_am_slave ) THEN
C                 -----
C                 case of A-1
C                 -----
C                 - Take columns with non-zero entry, say j,
C                 - to build Ej and store it in RHSCOMP
                  k=1              ! Column index in RHSCOMP
                  id%RHSCOMP(1_8:int(nbrhs_eff,8)*int(ld_rhscomp,8))
     &            = zero
                  ipos = 1
                  DO i = 1, nbcol_inbloc
                    colsize = irhs_ptr_copy(i+1) - irhs_ptr_copy(i)
                    IF (colsize.GT.0) THEN
                      ! Find global column index J and set
                      ! column K of RHSCOMP to ej (here IBEG is one)
                      j = i - 1 + jbeg_rhs
                      IF (do_permute_rhs.OR.interleave_par) THEN
                        j = perm_rhs(j)
                      ENDIF
                      iposrhscomp = id%POSINRHSCOMP_ROW(j)
C                     IF ( (IPOSRHSCOMP.LE.NB_FS_RHSCOMP_F)
C     &               .AND.(IPOSRHSCOMP.GT.0) ) THEN
                      IF (iposrhscomp.GT.0)  THEN
C                       Columns J corresponds to ej and thus to variable j
C                       that is on my proc
C                       Note that :
C                       In first entry in column
C                       we have and MUST have already scaled value of diagonal.
C                       This need have been done on master because we do not
C                       have scaling arrays available on slaves.
C                       Furthermore we know that only one entry is
C                       needed the diagonal entry (for the forward with A-1).
C
                        id%RHSCOMP(int(k-1,8)*int(ld_rhscomp,8)+
     &                       int(iposrhscomp,8)) =
     &                       rhs_sparse_copy(ipos)
                      ENDIF               ! End of J on my proc
                      k = k + 1
                      ipos = ipos + colsize  ! go to next column
                    ENDIF
                  ENDDO
                  IF (k.NE.nbrhs_eff+1) THEN
                    WRITE(6,*) 'Internal Error 9 in solution driver ',
     &              k,nbrhs_eff
                    call mumps_abort()
                  ENDIF
                ENDIF ! I_AM_SLAVE
C               -------
c               END A-1
C               -------
              ELSE
C               --------------
C               General sparse
C               --------------
C               -- reset to zero RHSCOMP for skipped columns (if any)
                IF ((keep(221).EQ.1).AND.(nb_rhsskipped.GT.0)
     &            .AND.i_am_slave) THEN
                  DO k = jbeg_rhs-nb_rhsskipped, jbeg_rhs-1
                    DO i = 1,  ld_rhscomp
                      id%RHSCOMP(int(k-1,8)*int(ld_rhscomp,8)
     &                           + int(i,8)) =  zero
                    ENDDO
                  ENDDO
                ENDIF
                IF (i_am_slave) THEN
                  DO k = 1, nbcol_inbloc
!                   it is equal to NBRHS_EFF in this case
                    kdec = int(k-1,8) * int(ld_rhscomp,8) +
     &                     ibeg_rhscomp - 1_8
                    id%RHSCOMP(kdec+1_8:kdec+nbent_rhscomp) = zero
                    DO iz=irhs_ptr_copy(k), irhs_ptr_copy(k+1)-1
                      i=irhs_sparse_copy(iz)
                      iposrhscomp = id%POSINRHSCOMP_ROW(i)
C                     Since all fully summed variables mapped
C                     on each proc are stored at the beginning
C                     of RHSCOMP, we can compare to KEEP(89)
C                     to know if RHSCOMP should be initialized
C                     So far the tree has not been pruned to exploit
C                     sparsity to compress RHSCOMP so we compare to
C                     NB_FS_RHSCOMP_TOT
                      IF ( (iposrhscomp.LE.nb_fs_rhscomp_tot)
     &                .AND.(iposrhscomp.GT.0) ) THEN
C                        ! I is fully summed var mapped on my proc
                        id%RHSCOMP(kdec+iposrhscomp)=
     &                  id%RHSCOMP(kdec+iposrhscomp) +
     &                  rhs_sparse_copy(iz)
                      ENDIF
                    ENDDO
                  ENDDO
                END IF ! I_AM_SLAVE
              ENDIF ! KEEP(237)
            ENDIF  ! ==== KEEP(248)==1 =====
C
        ELSE IF (i_am_slave) THEN
            ! I_AM_SLAVE AND (null space or Fwd in facto)
            IF (keep(111).NE.0) THEN
C             -----------------------
C             Null space computations
C             -----------------------
C
C             We are working on columns BEG_RHS:BEG_RHS+NBRHS_EFF-1
C             of RHS.
C             Columns in 1..KEEP(112):
C                    Put a one in corresponding
C                    position of the right-hand-side,
C                    and zeros in other places.
C             Columns in KEEP(112)+1: KEEP(112)+KEEP(17):
C                    root node => set
C                    0 everywhere and compute the local range
C                    corresponding to IBEG/IEND in root
C                    that will be passed to ZMUMPS_SEQ_SOLVE_ROOT_RR
C                    Also keep track of which part of
C                    ZMUMPS_RHS must be passed to
C                    ZMUMPS_SEQ_SOLVE_ROOT_RR.
C
              IF (keep(111).GT.0) THEN
                ibeg_glob_def = keep(111)
                iend_glob_def = keep(111)
              ELSE
                ibeg_glob_def = beg_rhs
                iend_glob_def = beg_rhs+nbrhs_eff-1
              ENDIF
              IF ( id%KEEP(112) .GT. 0 .AND. do_null_piv) THEN
                IF (ibeg_glob_def .GT.id%KEEP(112)) THEN
                  id%KEEP(235) = 0
                  do_null_piv = .false.
                ENDIF
                IF (ibeg_glob_def .LT.id%KEEP(112)
     &          .AND. iend_glob_def .GT.id%KEEP(112)
     &          .AND. do_null_piv ) THEN
C                 IEND_GLOB_DEF = id%KEEP(112)
C                 forcing exploit sparsity
C                 - cannot be done at this point
C                 - and is not what the user would have expected the
C                   code to to do anyway !!!!
C                 suppress:  id%KEEP(235) = 1  ! End Block of sparsity ON
                  do_null_piv = .false.
                ENDIF
              ENDIF
              IF (id%KEEP(235).NE.0) THEN
C               Exploit Sparsity in null space computations
C               We build /allocate the sparse RHS on MASTER
C               based on pivnul_list. Then we broadcast it
C               on the slaves
C               In this case we have ONLY ONE ENTRY per RHS
C
                nz_this_block=iend_glob_def-ibeg_glob_def+1
                ALLOCATE(irhs_ptr_copy(nz_this_block+1),stat=allocok)
                IF (allocok .GT.0 ) THEN
                  info(1)=-13
                  info(2)=nz_this_block
                  GOTO 50
                ENDIF
                irhs_ptr_copy_allocated = .true.
                ALLOCATE(irhs_sparse_copy(nz_this_block),stat=allocok)
                IF (allocok .GT.0 ) THEN
                  info(1)=-13
                  info(2)=nz_this_block
                  GOTO 50
                ENDIF
                irhs_sparse_copy_allocated=.true.
                nb_bytes = nb_bytes +
     &                         int(nz_this_block,8)*(k34_8+k34_8)
     &                         + k34_8
                nb_bytes_max = max(nb_bytes_max,nb_bytes)
                IF (id%MYID.eq.master) THEN
                  !  compute IRHS_PTR and IRHS_SPARSE_COPY
                  ii = 1
                  DO i = ibeg_glob_def, iend_glob_def
                    irhs_ptr_copy(i-ibeg_glob_def+1)      = i
                      irhs_sparse_copy(ii) = id%PIVNUL_LIST(i)
                    ii = ii +1
                  ENDDO
                  irhs_ptr_copy(nz_this_block+1) = nz_this_block+1
                ENDIF
C
C =====================  ERROR handling and propagation ================
 50             CONTINUE
                CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
                IF (info(1) .LT.0 ) GOTO 90
C ======================================================================
                CALL mpi_bcast(irhs_sparse_copy(1),
     &                nz_this_block,
     &                mpi_integer,
     &                master, id%COMM,ierr)
                CALL mpi_bcast(irhs_ptr_copy(1),
     &                nz_this_block+1,
     &                mpi_integer,
     &                master, id%COMM,ierr)
C               End IF Exploit Sparsity
              ENDIF
c
C             Initialize RHSCOMP to 0  ! to be suppressed
              DO k=1, nbrhs_eff
                kdec = int(k-1,8) * int(ld_rhscomp,8)
                id%RHSCOMP(kdec+1_8:kdec+int(ld_rhscomp,8))=zero
              END DO
C             Loop over the columns.
C             Note that if ( KEEP(220)+KEEP(109)-1 < IBEG_GLOB_DEF
C             .OR. KEEP(220) > IEND_GLOB_DEF ) then we do not enter
C             the loop.
C             Note that local processor has indices
C             KEEP(220):KEEP(220)+KEEP(109)-1
C
C               Computation of null space and computation of backward
C               step incompatible, do one or the other.
                DO i=max(ibeg_glob_def,keep(220)),
     &               min(iend_glob_def,keep(220)+keep(109)-1)
C                 Local processor is concerned by I-th column of
C                 global right-hand side.
                  jj= id%POSINRHSCOMP_ROW(id%PIVNUL_LIST(i-keep(220)+1))
                  IF (jj.GT.0) THEN
                    IF (keep(50).EQ.0) THEN
                      ! unsymmetric : always set to fixation
                      id%RHSCOMP( ibeg_rhscomp+
     &                    int(i-ibeg_glob_def,8)*int(ld_rhscomp,8) +
     &                    int(jj-1,8) ) =
     &                    cmplx(id%DKEEP(2),kind=kind(id%RHSCOMP))
                    ELSE
                      ! Symmetric: always set to one
                      id%RHSCOMP( ibeg_rhscomp+
     &                int(i-ibeg_glob_def,8)*int(ld_rhscomp,8)+
     &                int(jj-1,8) )=
     &                one
                    ENDIF
                  ENDIF
                ENDDO
              IF ( keep(17).NE.0 .AND.
     &             id%MYID_NODES.EQ.master_root) THEN
C               ---------------------------
C               Deficiency of the root node
C               Find range relative to root
C               ---------------------------
C               Among IBEG_GLOB_DEF:IEND_GLOB_DEF, find
C               intersection with KEEP(112)+1:KEEP(112)+KEEP(17)
                ibeg_root_def  = max(ibeg_glob_def,keep(112)+1)
                iend_root_def  = min(iend_glob_def,keep(112)+keep(17))
C               First column of right-hand side that must
C               be passed to ZMUMPS_SEQ_SOLVE_ROOT_RR is:
                iroot_def_rhs_col1 = ibeg_root_def-ibeg_glob_def + 1
C               We look for indices relatively to the root node,
C               substract number of null pivots outside root node
                ibeg_root_def = ibeg_root_def-keep(112)
                iend_root_def = iend_root_def-keep(112)
C               Note that if IBEG_ROOT_DEF > IEND_ROOT_DEF, then this
C               means that nothing must be done on the root node
C               for this set of right-hand sides.
              ELSE
                ibeg_root_def = -90999
                iend_root_def = -95999
                iroot_def_rhs_col1= 1
              ENDIF
            ELSE  ! End of null space (test on KEEP(111))
C             case of Fwd in facto
C             id%RHSCOMP need not be initialized. It will be set on the fly
C             to zero for normal fully summed variables of the fronts and
C             to -1 on the roots for the id%N+KEEP(253) variables added
C             to the roots.
            ENDIF ! End of null space (test on KEEP(111))
        ENDIF  ! I am slave
        timescatter2=mpi_wtime()-timescatter1+timescatter2
C       -------------------------------------------
C       Reserve space at the end of WORK_WCB on the
C       master of the root node. It will be used to
C       store the reduced RHS.
C       -------------------------------------------
        IF ( i_am_slave ) THEN
          lwcb8_sol_c = lwcb8
          IF ( id%MYID_NODES .EQ. master_root ) THEN
C           This is a special root (otherwise MASTER_ROOT < 0)
            IF ( associated(id%root%RHS_CNTR_MASTER_ROOT) ) THEN
C             RHS_CNTR_MASTER_ROOT may have been allocated
C             during the factorization phase.
              ptr_rhs_root => id%root%RHS_CNTR_MASTER_ROOT
#             if defined(MUMPS_F2003)
              lptr_rhs_root = size(id%root%RHS_CNTR_MASTER_ROOT,kind=8)
#             else
              lptr_rhs_root = int(size(id%root%RHS_CNTR_MASTER_ROOT),8)
#             endif
            ELSE
C             Otherwise, we use workspace in WCB
              lptr_rhs_root = int(nbrhs_eff,8) * int(size_root,8)
              ipt_rhs_root  = lwcb8 - lptr_rhs_root + 1_8
              ptr_rhs_root => work_wcb(ipt_rhs_root:lwcb8)
              lwcb8_sol_c = lwcb8_sol_c - lptr_rhs_root
            ENDIF
          ELSE
            lptr_rhs_root = 1_8
            ipt_rhs_root = lwcb8 ! Will be passed, but not accessed
            ptr_rhs_root => work_wcb(ipt_rhs_root:lwcb8)
            lwcb8_sol_c = lwcb8_sol_c - lptr_rhs_root
          ENDIF
        ENDIF
        IF (keep(221) .EQ. 2 ) THEN
C         Copy/send REDRHS in PTR_RHS_ROOT
C         (column by column if leading dimension LD_REDRHS
C         of REDRHS is not equal to SIZE_ROOT).
C         REDRHS was provided on the host
          IF ( ( id%MYID .EQ. master_root_in_comm ) .AND.
     &       ( id%MYID .EQ. master ) ) THEN
C         -- Same proc : copy is possible:
            ii = 0
            DO k=1, nbrhs_eff
              kdec = ibeg_redrhs+int(k-1,8)*int(ld_redrhs,8)-1_8
              DO i = 1, size_root
                ptr_rhs_root(ii+i) = id%REDRHS(kdec+i)
              ENDDO
              ii = ii+size_root
            ENDDO
          ELSE
C         -- send REDRHS
            IF ( id%MYID .EQ. master) THEN
C           -- send to MASTER_ROOT_IN_COMM using COMM communicator
C              assert: id%KEEP(116).EQ.SIZE_ROOT
              IF (ld_redrhs.EQ.size_root) THEN
C              --  One send
                 kdec = ibeg_redrhs
                 CALL mpi_send(id%REDRHS(kdec),
     &              size_root*nbrhs_eff,
     &              mpi_double_complex,
     &              master_root_in_comm, 0, id%COMM,ierr)
              ELSE
C             --  NBRHS_EFF sends
                DO k=1, nbrhs_eff
                  kdec = ibeg_redrhs+int(k-1,8)*int(ld_redrhs,8)
                  CALL mpi_send(id%REDRHS(kdec),size_root,
     &               mpi_double_complex,
     &               master_root_in_comm, 0, id%COMM,ierr)
                ENDDO
              ENDIF
            ELSE IF ( id%MYID .EQ. master_root_in_comm ) THEN
C            -- receive from MASTER
              ii = 1
              IF (ld_redrhs.EQ.size_root) THEN
C                -- receive all in on shot
                 CALL mpi_recv(ptr_rhs_root(ii),
     &              size_root*nbrhs_eff,
     &              mpi_double_complex,
     &              master, 0, id%COMM,status,ierr)
              ELSE
                DO k=1, nbrhs_eff
                  CALL mpi_recv(ptr_rhs_root(ii),size_root,
     &            mpi_double_complex,
     &            master, 0, id%COMM,status,ierr)
                  ii = ii + size_root
                ENDDO
              ENDIF
            ENDIF
C         -- other procs are not concerned
          ENDIF
        ENDIF
        timec1=mpi_wtime()
        IF ( i_am_slave ) THEN
          liw_passed = max( liw, 1 )
          la_passed  = max( la, 1_8 )
C
          IF ((id%KEEP(235).EQ.0).and.(id%KEEP(237).EQ.0) ) THEN
C
C         --- Normal case : we do not exploit sparsity of the RHS
C
            from_pp = .false.
            nbsparse_loc = (do_nbsparse.AND.nbrhs_eff.GT.1)
            pruned_size_loaded = 0_8  ! From ZMUMPS_SOL_ES module
            CALL zmumps_sol_c(id%root, id%N, id%S(1), la_passed, is(1),
     & liw_passed, work_wcb(1), lwcb8_sol_c, iwcb, liwcb, nbrhs_eff,
     & id%NA(1),id%LNA,id%NE_STEPS(1), srw3, mtype, icntl(1), from_pp,
     & id%STEP(1), id%FRERE_STEPS(1), id%DAD_STEPS(1), id%FILS(1),
     & id%PTLUST_S(1), id%PTRFAC(1), iwk_solve, liwk_solve, ptracb,
     & liwk_ptracb, id%PROCNODE_STEPS(1), id%NSLAVES, info(1),keep(1),
     & keep8(1), id%DKEEP(1), id%COMM_NODES, id%MYID, id%MYID_NODES,
     & bufr(1), lbufr, lbufr_bytes, id%ISTEP_TO_INIV2(1),
     & id%TAB_POS_IN_PERE(1,1), ibeg_root_def, iend_root_def,
     & iroot_def_rhs_col1, ptr_rhs_root(1), lptr_rhs_root, size_root,
     & master_root, id%RHSCOMP(ibeg_rhscomp), ld_rhscomp,
     & id%POSINRHSCOMP_ROW(1), id%POSINRHSCOMP_COL(1)
     & , 1, 1, 1, 1, idummy, 1, jdummy, kdummy, 1, ldummy, 1, mdummy
     & , 1, 1, nbsparse_loc, ptr_rhs_bounds(1), lptr_rhs_bounds
     & , 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%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%L0_OMP_FACTORS(1), id%LL0_OMP_FACTORS
     &    )
          ELSE
C           Exploit sparsity of the RHS (all cases)
C           Remark that JBEG_RHS is already initialized
C
            from_pp = .false.
            nbsparse_loc = (do_nbsparse.AND.nbrhs_eff.GT.1)
            CALL zmumps_sol_c(id%root, id%N, id%S(1), la_passed,is(1),
     & liw_passed,work_wcb(1),lwcb8_sol_c,iwcb,liwcb,nbrhs_eff,id%NA(1),
     & id%LNA,id%NE_STEPS(1),srw3,mtype,icntl(1),from_pp,id%STEP(1),
     & id%FRERE_STEPS(1), id%DAD_STEPS(1), id%FILS(1), id%PTLUST_S(1),
     & id%PTRFAC(1), iwk_solve, liwk_solve, ptracb, liwk_ptracb,
     & id%PROCNODE_STEPS(1),id%NSLAVES,info(1),keep(1), keep8(1),
     & id%DKEEP(1),id%COMM_NODES,id%MYID,id%MYID_NODES,bufr(1),lbufr,
     & lbufr_bytes, id%ISTEP_TO_INIV2(1), id%TAB_POS_IN_PERE(1,1),
     & ibeg_root_def,iend_root_def,iroot_def_rhs_col1,ptr_rhs_root(1),
     & lptr_rhs_root, size_root, master_root, id%RHSCOMP(ibeg_rhscomp),
     & ld_rhscomp, id%POSINRHSCOMP_ROW(1), id%POSINRHSCOMP_COL(1),
     & nz_this_block, nbcol_inbloc, id%NRHS, jbeg_rhs, id%Step2node(1),
     & id%KEEP(28),irhs_sparse_copy(1),irhs_ptr_copy(1), size(perm_rhs),
     & perm_rhs, size(uns_perm_inv), uns_perm_inv, nb_fs_rhscomp_f,
     & nb_fs_rhscomp_tot,nbsparse_loc,ptr_rhs_bounds(1),lptr_rhs_bounds
     & ,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%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%L0_OMP_FACTORS(1), id%LL0_OMP_FACTORS )
          ENDIF   ! end of exploit sparsity (pruning nodes of the tree)
        END IF
C       -----------------
C       End of slave code
C       -----------------
C
C
C       Propagate errors
        CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
        timec2=mpi_wtime()-timec1+timec2
C
C       Change error code.
        IF (info(1).eq.-2) then
          info(1)=-11
          IF (lpok)
     &    write(lp,*)
     &   ' WARNING : -11 error code obtained in solve'
        END IF
        IF (info(1).eq.-3) then
          info(1)=-14
          IF (lpok)
     &    write(lp,*)
     &    ' WARNING : -14 error code obtained in solve'
        END IF
C
C       Return in case of error.
        IF (info(1).LT.0) GO TO 90
C
C       ======================================================
C       ONLY FORWARD was performed (case of reduced RHS with Schur
C                               option during factorisation)
C       ======================================================
        IF ( keep(221) .EQ. 1 ) THEN ! === Begin OF REDUCED RHS ======
C         --------------------------------------
C         Send (or copy) reduced RHS from PTR_RHS_ROOT located on
C         MASTER_ROOT_IN_COMM to REDRHS located on MASTER (host node).
C         (column by column if leading dimension LD_REDRHS
C         of REDRHS is not equal to SIZE_ROOT)
C         --------------------------------------
          IF ( ( id%MYID .EQ. master_root_in_comm ) .AND.
     &        ( id%MYID .EQ. master ) ) THEN
C           -- same proc  --> copy
            ii = 0
            DO k=1, nbrhs_eff
              kdec = ibeg_redrhs+int(k-1,8)*int(ld_redrhs,8) - 1_8
              DO i = 1, size_root
                id%REDRHS(kdec+i) = ptr_rhs_root(ii+i)
              ENDDO
              ii = ii+size_root
            ENDDO
          ELSE
C           -- recv in REDRHS
            IF ( id%MYID .EQ. master ) THEN
C             -- recv from MASTER_ROOT_IN_COMM
              IF (ld_redrhs.EQ.size_root) THEN
C             --  One message to receive
                kdec = ibeg_redrhs
                CALL mpi_recv(id%REDRHS(kdec),
     &              size_root*nbrhs_eff,
     &              mpi_double_complex,
     &              master_root_in_comm, 0, id%COMM,
     &              status,ierr)
              ELSE
C             --  NBRHS_EFF receives
                DO k=1, nbrhs_eff
                  kdec = ibeg_redrhs+int(k-1,8)*int(ld_redrhs,8)
                  CALL mpi_recv(id%REDRHS(kdec),size_root,
     &              mpi_double_complex,
     &              master_root_in_comm, 0, id%COMM,
     &              status,ierr)
                ENDDO
              ENDIF
            ELSE IF ( id%MYID .EQ. master_root_in_comm ) THEN
C           -- send to MASTER
              ii = 1
              IF (ld_redrhs.EQ.size_root) THEN
C               -- send all in on shot
                CALL mpi_send(ptr_rhs_root(ii),
     &              size_root*nbrhs_eff,
     &              mpi_double_complex,
     &              master, 0, id%COMM,ierr)
              ELSE
                DO k=1, nbrhs_eff
                  CALL mpi_send(ptr_rhs_root(ii),size_root,
     &            mpi_double_complex,
     &            master, 0, id%COMM,ierr)
                  ii = ii + size_root
                ENDDO
              ENDIF
            ENDIF
C         -- other procs are not concerned
          ENDIF
        ENDIF ! ====== END OF REDUCED RHS (Fwd only performed) ======
C       =======================================================
C       BACKWARD was PERFORMED
C       Postprocess solution that is distributed
        IF ( keep(221) .NE. 1 ) THEN  ! BACKWARD was PERFORMED
C       -- KEEP(221).NE.1 => we are sure that backward has been performed
          IF (icntl21 == 0) THEN ! CENTRALIZED SOLUTION
C           ========================================================
C           GATHER SOLUTION computed during bwd
C           Each proc holds the pieces of solution corresponding
C           to all fully summed variables mapped on that processor
C           (i.e. corresponding to master nodes mapped on that proc)
C           In case of A-1 we gather directly in RHS_SPARSE
C           the distributed solution.
C           Scaling is done in all case on the fly of the reception
C           Note that when only FORWARD has been performed
C           RSH_MUMPS holds the solution computed during forward step
C           (ZMUMPS_SOL_R)
C           there is no need to copy back in RSH_MUMPS the solution
C           ========================================================
C           centralized solution
            IF (keep(237).EQ.0) THEN
C             CWORK not needed for AM1
              lcwork = max(max(keep(247),keep(246)),1)
              ALLOCATE( cwork(lcwork), stat=allocok )
              IF (allocok > 0) THEN
                info(1)=-13
                info(2)=max(max(keep(247),keep(246)),1)
              ENDIF
            ENDIF
            IF ( (id%MYID.EQ.master).AND. (keep(237).NE.0)
     &        .AND. (id%NSLAVES.NE.1)) THEN
C             Precompute map of indices in current column
C             (no need to reset it between columns
              ALLOCATE (map_rhs(id%N), stat = allocok)
              IF ( allocok .GT. 0 ) THEN
                IF (lpok) THEN
                 WRITE(lp,*) ' Problem allocation of MAP_RHS at solve'
                ENDIF
                info(1) = -13
                info(2) = id%N
              ELSE
                nb_bytes = nb_bytes + int(id%N,8) * k34_8
                nb_bytes_max = max(nb_bytes_max,nb_bytes)
              ENDIF
            ENDIF
C           Propagate errors
            CALL mumps_propinfo( icntl(1), info(1),
     &                       id%COMM,id%MYID)
C           Return in case of error.
            IF (info(1).LT.0) GO TO 90
            IF ((id%MYID.NE.master).OR. .NOT.lscal) THEN
             pt_scaling => dummy_scal
            ELSE
              IF (mtype.EQ.1) THEN
                pt_scaling => id%COLSCA
              ELSE
                pt_scaling => id%ROWSCA
              ENDIF
            ENDIF
            liw_passed = max( liw, 1 )
            timegather1=mpi_wtime()
            IF ( .NOT.i_am_slave ) THEN
C             I did not participate to computing part of the solution
C             (id%RHSCOMP not set/allocate) : receive solution, store
C             it and scale it.
              IF (keep(237).EQ.0) THEN
C               We need a workspace of minimal size KEEP(247)
C               in order to unpack pieces of the solution.
                CALL zmumps_gather_solution(id%NSLAVES,id%N,
     &              id%MYID, id%COMM, nbrhs_eff,
     &              mtype, id%RHS(1), ld_rhs, id%NRHS, jbeg_rhs,
     &              jdummy, id%KEEP(1), id%KEEP8(1),
     &              id%PROCNODE_STEPS(1), idummy, 1,
     &              id%STEP(1), bufr(1), lbufr, lbufr_bytes,
     &              cwork(1), lcwork,
     &              lscal, pt_scaling(1), size(pt_scaling),
     &              c_dummy, 1 , 1, idummy, 1,
     &              perm_rhs, size(perm_rhs) ! for sparse permuted RHS
     &              )
              ELSE
C               only gather target entries of A-1
                CALL zmumps_gather_solution_am1(id%NSLAVES,id%N,
     &              id%MYID, id%COMM, nbrhs_eff,
     &              c_dummy, 1, 1,
     &              id%KEEP(1), bufr(1), lbufr, lbufr_bytes,
     &              lscal, pt_scaling(1), size(pt_scaling)
C                   --- A-1 related entries
     &             ,irhs_ptr_copy(1), size(irhs_ptr_copy),
     &              irhs_sparse_copy(1), size(irhs_sparse_copy),
     &              rhs_sparse_copy(1), size(rhs_sparse_copy),
     &              uns_perm_inv, size(uns_perm_inv),
     &              idummy, 1, 0
     &              )
              ENDIF
            ELSE
C             Avoid temporary copy (IS(1)) that some old
C             compilers would do otherwise
              IF (keep(237).EQ.0) THEN
                IF (id%MYID.EQ.master) THEN
                  ptr_rhs => id%RHS
                  ncol_rhs_loc = id%NRHS
                  ld_rhs_loc   = ld_rhs
                  jbeg_rhs_loc = jbeg_rhs
                ELSE
                  ptr_rhs => cdummy_target
                  ncol_rhs_loc = 1
                  ld_rhs_loc   = 1
                  jbeg_rhs_loc = 1
                ENDIF
                CALL zmumps_gather_solution(id%NSLAVES,id%N,
     &          id%MYID, id%COMM, nbrhs_eff, mtype,
     &          ptr_rhs(1), ld_rhs_loc, ncol_rhs_loc, jbeg_rhs_loc,
     &          id%PTLUST_S(1), id%KEEP(1), id%KEEP8(1),
     &          id%PROCNODE_STEPS(1), is(1), liw_passed,
     &          id%STEP(1), bufr(1), lbufr, lbufr_bytes,
     &          cwork(1), lcwork,
     &          lscal, pt_scaling(1), size(pt_scaling),
     &          id%RHSCOMP(ibeg_rhscomp), ld_rhscomp, nbrhs_eff,
     &          id%POSINRHSCOMP_COL(1), id%N,
     &          perm_rhs, size(perm_rhs) ! For sparse permuted RHS
     &          )
              ELSE ! only gather target entries of A-1
                CALL zmumps_gather_solution_am1(id%NSLAVES,id%N,
     &          id%MYID, id%COMM, nbrhs_eff,
     &          id%RHSCOMP(ibeg_rhscomp), ld_rhscomp, nbrhs_eff,
     &          id%KEEP(1), bufr(1), lbufr, lbufr_bytes,
     &          lscal, pt_scaling(1), size(pt_scaling)
C               --- A-1 related entries
     &          , irhs_ptr_copy(1), size(irhs_ptr_copy),
     &          irhs_sparse_copy(1), size(irhs_sparse_copy),
     &          rhs_sparse_copy(1), size(rhs_sparse_copy),
     &          uns_perm_inv, size(uns_perm_inv),
     &          id%POSINRHSCOMP_COL(1), id%N, nb_fs_rhscomp_tot
     &          )
              ENDIF
            ENDIF
            timegather2=mpi_wtime()-timegather1+timegather2
            IF (keep(237).EQ.0)  DEALLOCATE( cwork )
            IF ( (id%MYID.EQ.master).AND. (keep(237).NE.0)
     &        ) THEN
C             Copy back solution from RHS_SPARSE_COPY TO RHS_SPARSE
              DO j = jbeg_rhs, jbeg_rhs+nbcol_inbloc-1
                IF (do_permute_rhs.OR.interleave_par) THEN
                 pj = perm_rhs(j)
                ELSE
                 pj =j
                ENDIF
                colsize = id%IRHS_PTR(pj+1) -
     &                id%IRHS_PTR(pj)
                IF (colsize.EQ.0) cycle
                jj = j-jbeg_rhs+1
C               Precompute map of indices in current column
C               (no need to reset it between columns
                IF (id%NSLAVES.NE.1) THEN
                 DO ii=1, colsize
                   map_rhs(id%IRHS_SPARSE(
     &                id%IRHS_PTR(pj) + ii - 1)) = ii
                 ENDDO
                 DO iz2 = irhs_ptr_copy(jj),irhs_ptr_copy(jj+1)-1
                     ii = irhs_sparse_copy(iz2)
                     id%RHS_SPARSE(id%IRHS_PTR(pj)+map_rhs(ii)-1)= 
     &                      rhs_sparse_copy(iz2)
                 ENDDO
                ELSE
C                 Entries within a column are in order
C                 IZ - Column index in Sparse RHS
                  DO iz= id%IRHS_PTR(pj), id%IRHS_PTR(pj+1)-1
                    iz2 = irhs_ptr_copy(jj) + 
     &                    iz - id%IRHS_PTR(pj)
                    id%RHS_SPARSE(iz) = rhs_sparse_copy(iz2)
                  ENDDO
                ENDIF
              ENDDO
              IF (id%NSLAVES.NE.1) THEN
               nb_bytes = nb_bytes - int(size(map_rhs),8) * k34_8
               DEALLOCATE ( map_rhs )
              ENDIF
            ENDIF             ! end A-1 on master
C
C         -- END of backward was performed with centralized solution
          ELSE   ! (KEEP(221).NE.1) .AND.(ICNTL21.NE.0))
C
C           BEGIN of backward performed with distributed solution
C           time local copy + scaling
            timecopyscale1=mpi_wtime()
C           The non working host should not do this:
            IF ( i_am_slave ) THEN
              liw_passed = max( liw, 1 )
C             Only called if more than 1 pivot
C             was eliminated by the processor.
C             Note that LSOL_loc >= KEEP(89)
              IF ( keep(89) .GT. 0 ) THEN
                CALL zmumps_distributed_solution(id%NSLAVES,
     &          id%N,id%MYID_NODES,
     &          mtype, id%RHSCOMP(ibeg_rhscomp), ld_rhscomp,
     &          nbrhs_eff, id%POSINRHSCOMP_COL(1),
     &          id%ISOL_loc(1), id%SOL_loc(1), id%NRHS,
     &          jbeg_rhs-nb_rhsskipped, id%LSOL_loc,
     &          id%PTLUST_S(1), id%PROCNODE_STEPS(1),
     &          id%KEEP(1),id%KEEP8(1),
     &          is(1), liw_passed,
     &          id%STEP(1), scaling_data_sol, lscal, nb_rhsskipped,
     &          perm_rhs, size(perm_rhs) ) ! For permuted sparse RHS
              ENDIF
            ENDIF
            timecopyscale2=mpi_wtime()-timecopyscale1+timecopyscale2
          ENDIF
C         === BACKWARD was PERFORMED WITH DISTRIBUTED SOLUTION ===
C         ========================================================
        ENDIF ! ==== END of BACKWARD was PERFORMED (KEEP(221).NE.1)
C       note that the main DO-loop on blocks is not ended yet
C
C       ============================================
C       BEGIN
C
C       ITERATIVE REFINEMENT AND/OR ERROR ANALYSIS
C
C       ============================================
        IF ( icntl10 > 0 .AND. nbrhs_eff > 1 ) THEN
C
C         ----------------------------------
C         Multiple RHS: apply a fixed number
C         of iterative refinement steps
C         ----------------------------------
C         DO I = 1, ICNTL10
            write(6,*) ' Internal ERROR 15 in sol_driver '
C           Compute residual:  Y <- SAVERHS - A * RHS
C           Solve RHS <- A^-1 Y, Y modified
C           Assemble in RHS(REDUCE)
C           RHS <- RHS + Y
C         END DO
        END IF
        IF (postpros) THEN
C
C         SAVERHS holds the original right hand side
C         Sparse rhs are saved in SAVERHS as dense rhs
C
C         * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C         Start iterative refinements. The master is managing the
C         organisation of work, but slaves are used to solve systems of
C         equations and, in case of distributed matrix, perform
C         matrix-vector products. It is more complicated to do this with
C         the SPMD version than it was with the master/slave approach.
C
C         * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
c         IF ( PROK  .AND. ICNTL10 .NE. 0 ) WRITE( MP, 270 )
          IF ( prokg .AND. icntl10 .NE. 0 ) WRITE( mpg, 270 )
C         Initializations and allocations
          nitref = abs(icntl10)
          ALLOCATE(r_y(id%N), stat = allocok)
          IF ( allocok .GT. 0 ) THEN
            info(1)=-13
            info(2)=id%N
            GOTO 777
          ENDIF
          nb_bytes = nb_bytes + int(id%N,8)*k16_8
          ALLOCATE(c_y(id%N), stat = allocok)
          IF ( allocok .GT. 0 ) THEN
            info(1)=-13
            info(2)=id%N
            GOTO 777
          ENDIF
          nb_bytes = nb_bytes + int(id%N,8)*k35_8
          IF ( id%MYID .EQ. master ) THEN
            ALLOCATE( iw1( 2 * id%N ),stat = allocok )
            IF ( allocok .GT. 0 ) THEN
              info(1)=-13
              info(2)=2 * id%N
              GOTO 777
            ENDIF
            nb_bytes = nb_bytes + int(2*id%N,8)*k34_8
            ALLOCATE( c_w(id%N), stat = allocok )
            IF ( allocok .GT. 0 ) THEN
              info(1)=-13
              info(2)=id%N
              GOTO 777
            ENDIF
            nb_bytes = nb_bytes + int(id%N,8)*k35_8
            ALLOCATE( r_w(2*id%N), stat = allocok )
            IF ( allocok .GT. 0 ) THEN
              info(1)=-13
              info(2)=id%N
              GOTO 777
            ENDIF
            nb_bytes = nb_bytes + int(2*id%N,8)*k16_8
            IF ( prokg .AND. icntl10 .GT. 0 )
     &      WRITE( mpg, 240) 'MAXIMUM NUMBER OF STEPS =', nitref
C           end allocations on Master
          END IF
          ALLOCATE(c_locwk54(id%N),stat = allocok)
          IF ( allocok .GT. 0 ) THEN
            info(1)=-13
            info(2)=id%N
            GOTO 777
          ENDIF
          nb_bytes = nb_bytes + int(id%N,8)*k35_8
          ALLOCATE(r_locwk54(id%N),stat = allocok)
          IF ( allocok .GT. 0 ) THEN
            info(1)=-13
            info(2)=id%N
            GOTO 777
          ENDIF
          nb_bytes = nb_bytes + int(id%N,8)*k16_8
          kase = 0
C         Synchro point with broadcast of errors
 777      CONTINUE
          nb_bytes_max = max(nb_bytes_max,nb_bytes)
          CALL mumps_propinfo( icntl(1), info(1),
     &                       id%COMM,id%MYID)
          IF ( info(1) .LT. 0 ) GOTO 90
C         TIMEEA needed if EA and IR with stopping criterium
C         and IR with fixed n.of steps. 
          timeea = 0.0d0
C         TIMEEA1 needed if EA and IR with fixed n.of steps  
          timeea1 = 0.0d0
          CALL mumps_secdeb(timeit)
C         -------------------------
C
C         RHSOL holds the initial guess for the solution
C         We start the loop on the Iterative refinement procedure
C
C
C
C           |-   IRefin.   L O O P   -|
C           V                         V
C
C  =========================================================
C    Computation of the infinity norm of A
C  =========================================================
          IF ((icntl11.GT.0).OR.(icntl10.GT.0)) THEN
C           We don't get through these lines if ICNTL10<=0 AND ICNTL11<=0
            IF ( keep(54) .eq. 0 ) THEN
C             ------------------
C             Centralized matrix
C             ------------------
              IF ( id%MYID .eq. master ) THEN
C               -----------------------------------------
C               Call ZMUMPS_SOL_X outside, if needed,
C               in order to compute w(i,2)=sum|Aij|,j=1:n
C               in vector R_W(id%N+i)
C               -----------------------------------------
                IF (keep(55).NE.0) THEN
C                 unassembled matrix and norm of row required
                  CALL zmumps_sol_x_elt(mtype, id%N,
     &            id%NELT, id%ELTPTR(1),
     &            id%LELTVAR, id%ELTVAR(1),
     &            id%KEEP8(30), id%A_ELT(1),
     &            r_w(id%N+1), keep(1),keep8(1) )
                ELSE
C                 assembled matrix
                  IF ( mtype .eq. 1 ) THEN
                    CALL zmumps_sol_x
     &              ( id%A(1), id%KEEP8(28), id%N, id%IRN(1), id%JCN(1),
     &              r_w(id%N+1), keep(1),keep8(1),
     &              0, id%SYM_PERM(1) )
                  ELSE
                    CALL zmumps_sol_x
     &              ( id%A(1), id%KEEP8(28), id%N, id%JCN(1), id%IRN(1),
     &              r_w(id%N+1), keep(1),keep8(1),
     &              0, id%SYM_PERM(1) )
                  END IF
                ENDIF
              ENDIF
            ELSE
C             ---------------------
C             Matrix is distributed
C             ---------------------
              IF ( i_am_slave .and.
     &             id%KEEP8(29) .NE. 0_8 ) THEN
                IF ( mtype .eq. 1 ) THEN
                  CALL zmumps_sol_x(id%A_loc(1),
     &            id%KEEP8(29), id%N,
     &            id%IRN_loc(1), id%JCN_loc(1),
     &            r_locwk54, id%KEEP(1),id%KEEP8(1), 
     &            0, id%SYM_PERM(1) )
                ELSE
                  CALL zmumps_sol_x(id%A_loc(1),
     &            id%KEEP8(29), id%N,
     &            id%JCN_loc(1), id%IRN_loc(1),
     &            r_locwk54, id%KEEP(1),id%KEEP8(1),
     &            0, id%SYM_PERM(1) )
                END IF
              ELSE
                r_locwk54 = rzero
              END IF
C             -------------------------
C             Assemble result on master
C             -------------------------
              IF ( id%MYID .eq. master ) THEN
                CALL mpi_reduce( r_locwk54, r_w( id%N + 1 ),
     &            id%N, mpi_double_precision,
     &            mpi_sum,master,id%COMM, ierr)
              ELSE
                CALL mpi_reduce( r_locwk54, r_dummy,
     &            id%N, mpi_double_precision,
     &            mpi_sum,master,id%COMM, ierr)
              END IF
C           End if KEEP(54)
            END IF
C
            IF ( id%MYID .eq. master ) THEN
C             R_W is available on the master process only
              rinfog(4) = dble(zero)
              DO i = 1, id%N
                rinfog(4) = max(r_w( id%N +i), rinfog(4))
              ENDDO
            ENDIF
C           end ICNTL11 =/0  v ICNTL10>0
          ENDIF
C  =========================================================
C    END norm of A
C  =========================================================
C         Initializations for the IR
          noiter = 0
          iflag_ir = 0
          testconv = .false.
C         Test of convergence should be made
          IF (( id%MYID .eq. master ).AND.(icntl10.GT.0)) THEN
            testconv = .true.
            arret = cntl(2)
            IF (arret .LT. 0.0d0) THEN
              arret = sqrt(epsilon(0.0d0))
            END IF
          ENDIF
C  =========================================================
C         Starting IR
          DO  22 irstep = 1, nitref +1
C  =========================================================
C
C  =========================================================
C   Refine the solution starting from the second step of do loop
C  =========================================================
            IF (( id%MYID .eq. master ).AND.(irstep.GT.1)) THEN
              noiter = noiter + 1
              DO i = 1, id%N
                 id%RHS(ibeg+i-1) = id%RHS(ibeg+i-1) + c_y(i)
              ENDDO
            ENDIF
C  ===========================================
C   Computation of the RESIDUAL and of |A||x|
C  ===========================================
            IF ( keep(54) .eq. 0 ) THEN
              IF ( id%MYID .eq. master ) THEN
                IF (keep(55).NE.0) THEN
C                 input matrix by element
                  CALL zmumps_eltyd( mtype, id%N,
     &            id%NELT, id%ELTPTR(1), id%LELTVAR,
     &            id%ELTVAR(1), id%KEEP8(30), id%A_ELT(1),
     &            saverhs, id%RHS(ibeg),
     &            c_y, r_w, keep(50))
                ELSE
                  IF ( mtype .eq. 1 ) THEN
                   CALL zmumps_sol_y(id%A(1), id%KEEP8(28),
     &                   id%N, id%IRN(1),
     &                   id%JCN(1), saverhs,
     &                   id%RHS(ibeg), c_y, r_w, keep(1),keep8(1))
                  ELSE
                   CALL zmumps_sol_y(id%A(1), id%KEEP8(28),
     &                        id%N, id%JCN(1),
     &                        id%IRN(1), saverhs,
     &                        id%RHS(ibeg), c_y, r_w, keep(1),keep8(1))
                  ENDIF
                ENDIF
              ENDIF
            ELSE
C             ---------------------
C             Matrix is distributed
C             ---------------------
              CALL mpi_bcast( rhs_ir(ibeg), id%N,
     &              mpi_double_complex, master,
     &              id%COMM, ierr )
C             --------------------------------------
C             Compute Y = SAVERHS - A * RHS
C             Y, SAVERHS defined only on master
C             --------------------------------------
              IF ( i_am_slave .and.
     &           id%KEEP8(29) .NE. 0_8 ) THEN
                CALL zmumps_loc_mv8( id%N, id%KEEP8(29),
     &          id%IRN_loc(1), id%JCN_loc(1), id%A_loc(1),
     &          rhs_ir(ibeg), c_locwk54, keep(50), mtype )
              ELSE
                c_locwk54 = zero
              END IF
              IF ( id%MYID .eq. master ) THEN
                CALL mpi_reduce( c_locwk54, c_y,
     &          id%N, mpi_double_complex,
     &          mpi_sum,master,id%COMM, ierr)
C              ===========================
                c_y = saverhs - c_y
C              ===========================
              ELSE
                CALL mpi_reduce( c_locwk54, c_dummy,
     &          id%N, mpi_double_complex,
     &          mpi_sum,master,id%COMM, ierr)
              END IF
C             --------------------------------------
C             Compute
C             * If MTYPE = 1
C                   W(i) = Sum | Aij | | RHSj |
C                           j
C             * If MTYPE = 0
C                   W(j) = Sum | Aij | | RHSi |
C                           i
C             R_LOCWK54 used as local array for W
C             RHS has been broadcasted
C             --------------------------------------
              IF ( i_am_slave .and. id%KEEP8(29) .NE. 0_8 ) THEN
                CALL zmumps_loc_omega1( id%N, id%KEEP8(29),
     &          id%IRN_loc(1), id%JCN_loc(1), id%A_loc(1),
     &          rhs_ir(ibeg), r_locwk54, keep(50), mtype )
              ELSE
                r_locwk54 = rzero
              END IF
              IF ( id%MYID .eq. master ) THEN
                CALL mpi_reduce( r_locwk54, r_w,
     &          id%N, mpi_double_precision,
     &          mpi_sum,master,id%COMM, ierr)
              ELSE
                CALL mpi_reduce( r_locwk54, r_dummy,
     &          id%N, mpi_double_precision,
     &          mpi_sum, master, id%COMM, ierr)
              ENDIF
            ENDIF
C  =====================================
C   END computation RESIDUAL and |A||x|
C  =====================================
            IF ( id%MYID .eq. master ) THEN
C
              IF ((icntl11.GT.0).OR.(icntl10.GT.0)) THEN
C             --------------
C             Error analysis and test of convergence,
C             Compute the sparse componentwise backward error:
C               - at each step if test of convergence of IR is
C                 requested (ICNTL(10)>0)
C               - at step 1 and NITREF+1 if error analysis
C                 to be computed (ICNTL(11)>0) and if ICNTL(10)< 0
                IF (((icntl11.GT.0).OR.((icntl10.LT.0).AND.
     &               ((irstep.EQ.1).OR.(irstep.EQ.nitref+1)))
     &               .OR.((icntl10.EQ.0).AND.(irstep.EQ.1)))
     &                      .OR.(icntl10.GT.0)) THEN
C                 Compute w1 and w2
C                 always if ICNTL10>0 in the other case if ICNTL11>0
C                 -----------------
                  IF (icntl10.LT.0) CALL mumps_secdeb(timeea1)
                  CALL zmumps_sol_omega(id%N,saverhs,
     &                id%RHS(ibeg), c_y, r_w, c_w, iw1, iflag_ir,
     &                rinfog(7), noiter, testconv,
     &                mp, arret, keep(361) )
                  IF (icntl10.LT.0) THEN
                    CALL mumps_secfin(timeea1)
                    id%DKEEP(120)=id%DKEEP(120)+timeea1
                  ENDIF
                ENDIF
                IF ((icntl11.GT.0).AND.(
     &          (icntl10.LT.0.AND.(irstep.EQ.1.OR.irstep.EQ.nitref+1))
     &          .OR.((icntl10.GE.0).AND.(irstep.EQ.1))
     &                           )) THEN
C                 Error analysis before iterative refinement
C                 or for last if icntl10<0
C                 ------------------------------------------
                  CALL mumps_secdeb(timeea)
                  IF (icntl10.EQ.0) THEN
C                   No IR : there will be only the EA of the 1st sol.
                    IF ( mpg .GT. 0 ) WRITE( mpg, 170 )
                  ELSEIF (irstep.EQ.1) THEN
C                   IR :  we print the EA of the 1st sol.
                    IF ( mpg .GT. 0 ) WRITE( mpg, 55 )
                  ELSEIF ((icntl10.LT.0).AND.(irstep.EQ.nitref+1)) THEN
C                   IR with fixed n. of steps:  we print the EA
C                                               of the last sol.
                    IF ( mpg .GT. 0 ) THEN
                      WRITE( mpg, 81 )
                      WRITE( mpg, * )
                      WRITE( mpg, 141 ) 
     &          'NUMBER OF STEPS OF ITERATIVE REFINEMENT REQUESTED =', 
     &          noiter
                    ENDIF
                  ENDIF
                  givsol = .true.
                  CALL zmumps_sol_q(mtype,info(1),id%N,
     &            id%RHS(ibeg),
     &            saverhs,r_w(id%N+1),c_y,givsol,
     &            rinfog(4),rinfog(5),rinfog(6),mpg,icntl(1),
     &            keep(1),keep8(1))
                  IF ( mpg .GT. 0 ) THEN
C                   Error analysis before iterative refinement
                    WRITE( mpg, 115 )
     &              'RINFOG(7):COMPONENTWISE SCALED RESIDUAL(W1)=',
     &              rinfog(7)
                    WRITE( mpg, 115 )
     &              '------(8):---------------------------- (W2)=',
     &              rinfog(8)
                  END IF
                  CALL mumps_secfin(timeea)
                  id%DKEEP(120)=id%DKEEP(120)+timeea
C                 end EA of the first solution
                END IF
              END IF
C             --------------
              IF (irstep.EQ.nitref +1) THEN
C               If we are at the NITREF+1 step , we have refined the
C               solution NITREF times so we have to stop.
                kase = 0
C               If we test the convergence (ICNTL10.GT.0) and 
C               IFLAG_IR = 0 we set a warning : more than NITREF steps
C               needed
                IF ((icntl10.GT.0).AND.(iflag_ir.EQ.0))
     &             id%INFO(1) = id%INFO(1) + 8
              ELSE
                IF (icntl10.GT.0) THEN
C                 -------------------
C                 Results of the test of convergence.
C                 IFLAG_IR =  0 we should try to improve the solution
C                          =  1 the stopping criterium is satisfied
C                          =  2 the method is diverging, we go back
C                               to the previous iterate
C                          =  3 the convergence is too slow
                  IF (iflag_ir.GT.0) THEN
C                   If the convergence criterion is satisfied
C                   or the convergence too slow
C                   we set KASE=0 (end of the Iterative refinement)
                    kase = 0
C                   If the convergence is not improved,
C                   we go back to the previous iterate.
C                   IFLAG_IR can be equal to 2 only if IRStep >= 2
                    IF (iflag_ir.EQ.2)  noiter = noiter - 1
                  ELSE
C                   IFLAG_IR=0, try to improve the solution
                    kase = 2
                  ENDIF
                ELSEIF (icntl10.LT.0) THEN
C                 -------------------
                  kase = 2
                ELSE
C                 ICNTL10 = 0, we want to perform only EA and not IR.
C                 -----------------
                  kase = 0
                END IF
              ENDIF
C           End Master
            ENDIF
C           --------------
C           Broadcast KASE
C           --------------
            CALL mpi_bcast( kase, 1, mpi_integer, master,
     &                      id%COMM, ierr )
C           If Kase= 0 we quit the IR process
            IF (kase.LE.0) GOTO 666
            IF (kase.LT.0) THEN
              WRITE(*,*) "Internal error 17 in ZMUMPS_SOL_DRIVER"
            ENDIF
C  =========================================================
C   COMPUTE the solution of Ay = r
C  =========================================================
C           Call internal routine to avoid code duplication
            CALL zmumps_pp_solve()
            IF (info(1) .LT. 0) GOTO 90
C           -----------------------
C           Go back to beginning of
C           loop to apply next step
C           of iterative refinement
C           -----------------------
  22      CONTINUE
 666      CONTINUE
C         ************************************************
C
C         End of the iterative refinement procedure
C
C         ************************************************
          CALL mumps_secfin(timeit)
          IF ( id%MYID .EQ. master ) THEN
            IF ( nitref .GT. 0 ) THEN
              id%INFOG(15) = noiter
            END IF
C           id%DKEEP(114) time for the iterative refinement
C           id%DKEEP(120) time for the error analysis
C           id%DKEEP(121) time for condition number
C           these values are meaningful only on the host.
            IF (icntl10.EQ.0) THEN
C           No IR has been requested. All the time is needed
C           for computing EA
               id%DKEEP(120)=timeit
            ELSE
C            IR has been requested
             id%DKEEP(114)=timeit - id%DKEEP(120)
            ENDIF
          END IF
          IF ( prokg ) THEN
              IF (icntl10.GT.0) THEN
                WRITE( mpg, 81 )
                WRITE( mpg, * )
                WRITE( mpg, 141 )
     &          'NUMBER OF STEPS OF ITERATIVE REFINEMENTS PERFORMED  =',
     &          noiter
              ENDIF
          ENDIF
C
C         ==================================================
C         BEGIN
C         Perform error analysis after iterative refinement
C         ==================================================
          IF ((icntl11 .GT. 0).AND.(icntl10.GT.0)) THEN
C           If IR is requested with test of convergence,
C           the EA of the last step of IR is done here, 
C           otherwise EA of the last step is done at the 
C           end of IR
            CALL mumps_secdeb(timeea)
            kase = 0
            IF (id%MYID .eq. master ) THEN
C             Test if IFLAG_IR = 2, that is if the the IR was diverging,
C             we went back to the previous iterate
C             We have to do EA on the last computed solution.
              IF (iflag_ir.EQ.2) kase = 2
            ENDIF
C           --------------
C           Broadcast KASE
C           --------------
            CALL mpi_bcast( kase, 1, mpi_integer, master,
     &      id%COMM, ierr )
            IF (kase.EQ.2) THEN
C             We went back to the previous iterate
C             We have to do EA on the last computed solution.
C             Compute the residual in C_Y using IRN, JCN, ASPK
C             and the solution  RHS(IBEG)
C             The norm of the ith row in R_Y(I).
              IF ( keep(54) .eq. 0 ) THEN
C               ---------------------
C               Matrix is centralized
C               ---------------------
                IF (id%MYID .EQ. master) THEN
                  IF (keep(55).EQ.0) THEN
                    CALL zmumps_qd2( mtype, id%N, id%KEEP8(28), id%A(1),
     &              id%IRN(1), id%JCN(1),
     &              id%RHS(ibeg), saverhs, r_y, c_y, keep(1),keep8(1))
                  ELSE
                    CALL zmumps_eltqd2( mtype, id%N,
     &              id%NELT, id%ELTPTR(1),
     &              id%LELTVAR, id%ELTVAR(1),
     &              id%KEEP8(30), id%A_ELT(1),
     &              id%RHS(ibeg), saverhs, r_y, c_y, keep(1),keep8(1))
                  ENDIF
                ENDIF
              ELSE
C               ---------------------
C               Matrix is distributed
C               ---------------------
                CALL mpi_bcast( rhs_ir(ibeg), id%N,
     &              mpi_double_complex, master,
     &              id%COMM, ierr )
C               ----------------
C               Compute residual
C               ----------------
                IF ( i_am_slave .and.
     &            id%KEEP8(29) .NE. 0_8 ) THEN
                  CALL zmumps_loc_mv8( id%N, id%KEEP8(29),
     &            id%IRN_loc(1), id%JCN_loc(1), id%A_loc(1),
     &            rhs_ir(ibeg), c_locwk54, keep(50), mtype )
                ELSE
                  c_locwk54 = zero
                END IF
                IF ( id%MYID .eq. master ) THEN
                  CALL mpi_reduce( c_locwk54, c_y,
     &            id%N, mpi_double_complex,
     &            mpi_sum,master,id%COMM, ierr)
                  c_y = saverhs - c_y
                ELSE
                  CALL mpi_reduce( c_locwk54, c_dummy,
     &            id%N, mpi_double_complex,
     &            mpi_sum,master,id%COMM, ierr)
                END IF
              ENDIF
            ENDIF  ! KASE.EQ.2
            IF (id%MYID .EQ. master) THEN
C             Compute which equations are associated to w1 and which
C             ones are associated to w2 in case of IFLAG_IR=2.
C             If IFLAG_IR = 0 or 1 IW1 should be correct
              IF (iflag_ir.EQ.2) THEN
                testconv = .false.
                CALL zmumps_sol_omega(id%N,saverhs,
     &              id%RHS(ibeg), c_y, r_w, c_w, iw1, iflag_ir,
     &              rinfog(7), 0, testconv,
     &              mp, arret, keep(361) )
              ENDIF ! (IFLAG_IR.EQ.2)
c             Compute some statistics for
              givsol = .true.
              CALL zmumps_sol_q(mtype,info(1),id%N,
     &        id%RHS(ibeg),
     &        saverhs,r_w(id%N+1),c_y,givsol,
     &        rinfog(4),rinfog(5),rinfog(6),mpg,icntl(1),
     &        keep(1),keep8(1))
            ENDIF ! Master
            CALL mumps_secfin(timeea)
            id%DKEEP(120)=id%DKEEP(120)+timeea
          ENDIF ! ICNTL11>0 and ICNTL10>0
C  =========================================================
C   Compute the Condition number associated if requested.
C  =========================================================
          CALL mumps_secdeb(timelcond)
          IF (icntl11 .EQ. 1) THEN
            IF ( id%MYID .eq. master ) THEN
C             Notice that D is always the identity
              ALLOCATE( d(id%N),stat =allocok )
              IF ( allocok .GT. 0 ) THEN
                info(1)=-13
                info(2)=id%N
                GOTO 777
              ENDIF
              nb_bytes = nb_bytes + int(id%N,8)*k16_8
              DO i = 1, id%N
                d( i ) = rone
              END DO
            ENDIF
            kase = 0
 222        CONTINUE
            IF ( id%MYID .EQ. master ) THEN
              CALL zmumps_sol_lcond(id%N, saverhs,
     &        id%RHS(ibeg), c_y, d, r_w, c_w, iw1, kase,
     &        rinfog(7), rinfog(9), rinfog(10),
     &        mp, keep(1),keep8(1))
            ENDIF
C           --------------
C           Broadcast KASE
C           --------------
            CALL mpi_bcast( kase, 1, mpi_integer, master,
     &                      id%COMM, ierr )
C           KASE <= 0
C           We reach the end of iterative method to compute
C           LCOND1 and LCOND2
            IF (kase.LE.0) GOTO 224
            CALL zmumps_pp_solve()
            IF (info(1) .LT. 0) GOTO 90
C           ---------------------------
C           Go back to beginning of
C           loop to apply next step
C           of iterative method
C           -----------------------
            GO TO 222
C    End ICNTL11 = 1
          ENDIF
 224      CONTINUE
          CALL mumps_secfin(timelcond)
          id%DKEEP(121)=id%DKEEP(121)+timelcond
          IF ((id%MYID .EQ. master).AND.(icntl11.GT.0)) THEN
            IF (icntl10.GT.0) THEN
C             If ICNTL10<0 these stats have been printed before IR
              IF ( mpg .GT. 0 ) THEN
                WRITE( mpg, 115 )
     &          'RINFOG(7):COMPONENTWISE SCALED RESIDUAL(W1)=',
     &          rinfog(7)
                WRITE( mpg, 115 )
     &          '------(8):---------------------------- (W2)=',
     &          rinfog(8)
              ENDIF
            END IF
            IF (icntl11.EQ.1) THEN
C           If ICNTL11/=1 these stats haven't been computed
              IF (mpg.GT.0) THEN
               WRITE( mpg, 115 )
     &         '------(9):Upper bound ERROR ...............=',
     &         rinfog(9)
               WRITE( mpg, 115 )
     &         '-----(10):CONDITION NUMBER (1) ............=',
     &         rinfog(10)
               WRITE( mpg, 115 )
     &         '-----(11):CONDITION NUMBER (2) ............=',
     &         rinfog(11)
              END IF
            END IF
          END IF ! MASTER && ICNTL11.GT.0
          IF ( prokg .AND. abs(icntl10) .GT.0 ) WRITE( mpg, 131 )
C===================================================
C Perform error analysis after iterative refinements
C END
C===================================================
C
          IF (id%MYID == master) THEN
            nb_bytes = nb_bytes - int(size(c_w),8)*k35_8
            DEALLOCATE(c_w)
            nb_bytes = nb_bytes - int(size(r_w),8)*k16_8
     &                        - int(size(iw1),8)*k34_8
            DEALLOCATE(r_w)
            DEALLOCATE(iw1)
            IF (icntl11 .EQ. 1) THEN
C             We have used D only for LCOND1,2
              nb_bytes = nb_bytes - int(size(d  ),8)*k16_8
              DEALLOCATE(d)
            ENDIF
          ENDIF
          nb_bytes = nb_bytes -
     &     (int(size(r_y),8)+int(size(r_locwk54),8))*k16_8
          nb_bytes = nb_bytes -
     &     (int(size(c_y),8)+int(size(c_locwk54),8))*k35_8
          DEALLOCATE(r_y)
          DEALLOCATE(c_y)
          DEALLOCATE(r_locwk54)
          DEALLOCATE(c_locwk54)
C End POSTPros
        END IF
C============================================
C
C  ITERATIVE REFINEMENT AND/OR ERROR ANALYSIS
C
C  END
C
C============================================
C       ==========================
C       Begin reordering on master
C       corresponding to maximum transversal permutation
C       in case of centralized solution
C       (ICNTL21==0)
C
        IF ( id%MYID .EQ. master .AND. icntl21==0
     &     .AND. keep(23) .NE. 0.AND.keep(237).EQ.0) THEN
C         ((No transpose and backward performed and NO A-1)
C         or null space computation): permutation
C         must be done on solution.
          IF ((keep(221).NE.1 .AND. mtype .EQ. 1)
     &       .OR. keep(111) .NE.0 .OR. keep(252).NE.0 ) THEN
C           Permute the solution RHS according to the column
C           permutation held in UNS_PERM
C           Column J of the permuted matrix corresponds to
C           column UNS_PERM(J) of the original matrix.
C           RHS holds the permuted solution
C           Note that id%N>1 since KEEP(23)=0 when id%N=1
C
            ALLOCATE( c_rw1( id%N ),stat =allocok )
!           temporary not in NB_BYTES
            IF ( allocok .GT. 0 ) THEN
              info(1)=-13
              info(2)=id%N
              WRITE(*,*) 'could not allocate ', id%N, 'integers.'
              CALL mumps_abort()
            END IF
            DO k = 1, nbrhs_eff
              IF (keep(242).EQ.0) THEN
                kdec = (k-1)*ld_rhs+ibeg-1
              ELSE
C               -------------------------------
C               Columns just computed might not
C               be contiguous in original RHS
C               -------------------------------
                kdec = int(perm_rhs(k-1+jbeg_rhs)-1,8)*int(ld_rhs,8)
              ENDIF
              DO i = 1, id%N
                c_rw1(i) = id%RHS(kdec+i)
              ENDDO
              DO i = 1, id%N
               jperm = id%UNS_PERM(i)
               id%RHS( kdec+jperm ) = c_rw1( i )
              ENDDO
            ENDDO
            DEALLOCATE( c_rw1 ) !temporary not in NB_BYTES
          END IF
        END IF
C
C  End reordering on master
C  ========================
        IF (id%MYID.EQ.master .and.icntl21==0.and.keep(221).NE.1.AND.
     &     (keep(237).EQ.0) ) THEN
*         print out the solution
          IF ( info(1) .GE. 0 .AND. icntl(4).GE.3 .AND. icntl(3).GT.0)
     &    THEN
            k = min0(10, id%N)
            IF (icntl(4) .eq. 4 ) k = id%N
            j = min0(10,nbrhs_eff)
            IF (icntl(4) .eq. 4 ) j = nbrhs_eff
            DO ii=1, j
              WRITE(icntl(3),110) beg_rhs+ii-1
              WRITE(icntl(3),160)
     &      (id%RHS(ibeg+(ii-1)*ld_rhs+i-1),i=1,k)
            ENDDO
          END IF
        END IF
C     ==========================
C     blocking for multiple RHS (END OF DO WHILE (BEG_RHS.LE.NBRHS)
        IF ((keep(248).EQ.1).AND.(keep(237).EQ.0)) THEN
          ! case of general sparse: in case of empty columns
          ! nbrhs_eff might has been updated and broadcasted
          ! and holds the effective size of a contiguous block of
          ! non empty columns
          beg_rhs = beg_rhs + nbrhs_eff ! nb of nonempty columns
        ELSE
          beg_rhs = beg_rhs + nbrhs
        ENDIF
C     }
      ENDDO
C     END DO WHILE (BEG_RHS.LE.id%NRHS)
C     =================================
      IF (keep(400) .GT. 0) THEN
        CALL zmumps_sol_l0omp_ld(keep(400))
      ENDIF
C
C     ========================================================
C     Reset RHS to zero for all remaining columns that
C     have not been processed because they were emtpy
C     ========================================================
      IF (   (id%MYID.EQ.master)
     &       .AND. ( keep(248).NE.0 )  ! sparse RHS on input
     &       .AND. ( keep(237).EQ.0 )  ! No A-1
     &       .AND. ( icntl21.EQ.0 )    ! Centralized solution
     &       .AND. ( keep(221) .NE.1 ) ! Not Reduced RHS step of Schur
     &       .AND. ( jend_rhs .LT. id%NRHS )
     &   )
     &         THEN
        jbeg_new = jend_rhs + 1
        IF (do_permute_rhs.OR.interleave_par) THEN
           DO WHILE ( jbeg_new.LE. id%NRHS)
              DO i=1, id%N
                 id%RHS(int(perm_rhs(jbeg_new) -1,8)*int(ld_rhs,8)+i)
     &                     = zero
              ENDDO
              jbeg_new = jbeg_new +1
           ENDDO
        ELSE
          DO WHILE ( jbeg_new.LE. id%NRHS)
            DO i=1, id%N
                id%RHS(int(jbeg_new -1,8)*int(ld_rhs,8) + i) = zero
            ENDDO
            jbeg_new = jbeg_new +1
          ENDDO
        ENDIF                   ! End DO_PERMUTE_RHS.OR.INTERLEAVE_PAR
      ENDIF
C     ========================================================
C     Reset id%SOL_loc to zero for all remaining columns that
C     have not been processed because they were emtpy
C     ========================================================
      IF ( i_am_slave .AND. (icntl21.NE.0) .AND.
     &      ( jend_rhs .LT. id%NRHS ) .AND. keep(221).NE.1 ) THEN
        jbeg_new = jend_rhs + 1
        IF (do_permute_rhs.OR.interleave_par) THEN
           DO WHILE ( jbeg_new.LE. id%NRHS)
              DO i=1, keep(89)
                 id%SOL_loc(int(perm_rhs(jbeg_new) -1,8)*
     &                      int(id%LSOL_loc,8)+int(i,8)) = zero
              ENDDO
              jbeg_new = jbeg_new +1
           ENDDO
        ELSE
C
           DO WHILE ( jbeg_new.LE. id%NRHS)
            DO i=1, keep(89)
                id%SOL_loc((jbeg_new -1)*id%LSOL_loc + i) = zero
            ENDDO
            jbeg_new = jbeg_new +1
           ENDDO
        ENDIF
      ENDIF
C
C     ================================================================
C     Reset id%RHSCOMP and id%REDRHS to zero for all remaining columns
C     that have not been processed because they were emtpy
C     ================================================================
      IF ((keep(221).EQ.1) .AND.
     &        ( jend_rhs .LT. id%NRHS ) ) THEN
       IF (id%MYID .EQ. master) THEN
           jbeg_new = jend_rhs + 1
           DO WHILE ( jbeg_new.LE. id%NRHS)
            DO i=1,  id%SIZE_SCHUR
              id%REDRHS(int(jbeg_new -1,8)*int(ld_redrhs,8) +
     &                  int(i,8)) =  zero
            ENDDO
            jbeg_new = jbeg_new +1
           ENDDO
       ENDIF
       IF (i_am_slave) THEN
           jbeg_new = jend_rhs + 1
           DO WHILE ( jbeg_new.LE. id%NRHS)
            DO i=1,nbent_rhscomp
              id%RHSCOMP(int(jbeg_new -1,8)*int(ld_rhscomp,8) +
     &                   int(i,8)) =  zero
            ENDDO
            jbeg_new = jbeg_new +1
           ENDDO
       ENDIF
      ENDIF
C
C
C     ! maximum size used on that proc
      id%INFO(26) = int(nb_bytes_max / 1000000_8)
C     Centralize memory statistics on the host
C
C       INFOG(30) = size of mem in bytes for solve
C                   for the processor using largest memory
C       INFOG(31) = size of mem in bytes for solve
C                   sum over all processors
C     ----------------------------------------------------
      CALL mumps_mem_centralize( id%MYID, id%COMM,
     &                           id%INFO(26), id%INFOG(30), irank )
      IF ( prokg ) THEN
       IF (print_maxavg) THEN
        WRITE( mpg,'(A,I10) ')
     &  ' ** Rank of processor needing largest memory in solve     :',
     &  irank
        WRITE( mpg,'(A,I10) ')
     &  ' ** Space in MBYTES used by this processor for solve      :',
     &  id%INFOG(30)
        IF ( keep(46) .eq. 0 ) THEN
        WRITE( mpg,'(A,I10) ')
     &  ' ** Avg. Space in MBYTES per working proc during solve    :',
     &  ( id%INFOG(31)-id%INFO(26) ) / id%NSLAVES
        ELSE
        WRITE( mpg,'(A,I10) ')
     &  ' ** Avg. Space in MBYTES per working proc during solve    :',
     &  id%INFOG(31) / id%NSLAVES
        END IF
       ELSE
        WRITE( mpg,'(A,I10) ')
     &  ' ** Space in MBYTES used for solve                        :',
     &  id%INFOG(30)
       ENDIF
      END IF
*===============================
*End of Solve Phase
*===============================
C  Store and print timings
      CALL mumps_secfin(time3)
      id%DKEEP(112)=time3
      id%DKEEP(113)=timec2
      id%DKEEP(115)=timescatter2
      id%DKEEP(116)=timegather2
      id%DKEEP(122)=timecopyscale2
C     Reductions of DKEEP(115,116,117,118,119,122):
      CALL mpi_reduce( id%DKEEP(115), id%DKEEP(160),1,
     &mpi_double_precision, mpi_max, master, id%COMM, ierr )
      CALL mpi_reduce( id%DKEEP(116), id%DKEEP(161),1,
     &mpi_double_precision, mpi_max, master, id%COMM, ierr )
      CALL mpi_reduce( id%DKEEP(117), id%DKEEP(162),1,
     &mpi_double_precision, mpi_max, master, id%COMM, ierr )
      CALL mpi_reduce( id%DKEEP(118), id%DKEEP(163),1,
     &mpi_double_precision, mpi_max, master, id%COMM, ierr )
      CALL mpi_reduce( id%DKEEP(119), id%DKEEP(164),1,
     &mpi_double_precision, mpi_max, master, id%COMM, ierr )
      CALL mpi_reduce( id%DKEEP(122), id%DKEEP(165),1,
     &mpi_double_precision, mpi_max, master, id%COMM, ierr )
C
      IF (prokg) THEN
        WRITE ( mpg, *)
        WRITE ( mpg, *) "Leaving solve with ..."
        WRITE( mpg, 434 ) id%DKEEP(160) ! max id%DKEEP(115)
        WRITE( mpg, 432 ) id%DKEEP(113) ! ok without reduction
        WRITE( mpg, 435 ) id%DKEEP(162) ! max id%DKEEP(117)
        IF ((keep(38).NE.0).OR.(keep(20).NE.0))
     &     WRITE( mpg, 437 ) id%DKEEP(164) ! id%DKEEP(119)
        WRITE( mpg, 436 ) id%DKEEP(163) ! id%DKEEP(118)
        WRITE( mpg, 433 ) id%DKEEP(161) ! max(DKEEP(116)) -- Gather
        WRITE( mpg, 431 ) id%DKEEP(165) ! max(DKEEP(122)) -- Dist. sol.
      ENDIF
      IF ( prok ) THEN
        WRITE ( mp, *)
        WRITE ( mp, *) "Local statistics"
        WRITE( mp, 434 ) id%DKEEP(115)
        WRITE( mp, 432 ) id%DKEEP(113)
        WRITE( mp, 435 ) id%DKEEP(117)
        IF ((keep(38).NE.0).OR.(keep(20).NE.0))
     &     WRITE( mp, 437 ) id%DKEEP(119)
        WRITE( mp, 436 ) id%DKEEP(118)
        WRITE( mp, 433 ) id%DKEEP(116)
        WRITE( mp, 431 ) id%DKEEP(122)
      END IF
 90   CONTINUE
      IF (info(1) .LT.0 ) THEN
      ENDIF
      IF (keep(485) .EQ. 1) THEN
        keep(350) = keep350_save
        IF (is_lr_mod_to_struc_done) THEN
          CALL zmumps_blr_mod_to_struc(id%BLRARRAY_ENCODING)
          CALL mumps_fdm_mod_to_struc('F',id%FDM_F_ENCODING,
     &             id%INFO(1))
        ENDIF
      ENDIF
      IF (keep(201).GT.0)THEN
        IF (is_init_ooc_done) THEN
          CALL zmumps_ooc_end_solve(ierr)
          IF (ierr.LT.0 .AND. info(1) .GE. 0) info(1) = ierr
        ENDIF
        CALL mumps_propinfo( icntl(1), info(1),
     &         id%COMM,id%MYID)
      ENDIF
C     ------------------------
C     Check allocation before
C     to deallocate (cases of
C     errors that could happen
C     before or after allocate
C     statement)
C
C       Sparse RHS
C       Free space and reset pointers if needed
        IF (irhs_sparse_copy_allocated) THEN
             nb_bytes =  nb_bytes -
     &        int(size(irhs_sparse_copy),8)*k34_8
             DEALLOCATE(irhs_sparse_copy)
             irhs_sparse_copy_allocated=.false.
             NULLIFY(irhs_sparse_copy)
       ENDIF
       IF (irhs_ptr_copy_allocated) THEN
             nb_bytes =  nb_bytes -
     &        int(size(irhs_ptr_copy),8)*k34_8
             DEALLOCATE(irhs_ptr_copy)
             irhs_ptr_copy_allocated=.false.
             NULLIFY(irhs_ptr_copy)
       ENDIF
       IF (rhs_sparse_copy_allocated) THEN
             nb_bytes =  nb_bytes -
     &        int(size(rhs_sparse_copy),8)*k35_8
             DEALLOCATE(rhs_sparse_copy)
             rhs_sparse_copy_allocated=.false.
             NULLIFY(rhs_sparse_copy)
       ENDIF
       IF (allocated(map_rhs_loc)) THEN
         nb_bytes = nb_bytes - int(size(map_rhs_loc),8)*k34_8
         DEALLOCATE(map_rhs_loc)
       ENDIF
       IF (irhs_loc_ptr_allocated ) THEN
         nb_bytes = nb_bytes - int(size(irhs_loc_ptr),8)*k34_8
         DEALLOCATE(irhs_loc_ptr)
         NULLIFY(irhs_loc_ptr)
         irhs_loc_ptr_allocated = .false.
       ENDIF
       IF (i_am_slave.AND.lscal.AND.keep(248).EQ.-1) THEN
         nb_bytes = nb_bytes -
     &   int(size(scaling_data_dr%SCALING_LOC),8)*k16_8
         DEALLOCATE(scaling_data_dr%SCALING_LOC)
         NULLIFY   (scaling_data_dr%SCALING_LOC)
       ENDIF
      IF (allocated(perm_rhs)) THEN
        nb_bytes = nb_bytes - int(size(perm_rhs),8)*k34_8
        DEALLOCATE(perm_rhs)
      ENDIF
C     END A-1
      IF (allocated(uns_perm_inv)) THEN
        nb_bytes = nb_bytes - int(size(uns_perm_inv),8)*k34_8
        DEALLOCATE(uns_perm_inv)
      ENDIF
      IF (allocated(bufr)) THEN
          nb_bytes = nb_bytes - int(size(bufr),8)*k34_8
          DEALLOCATE(bufr)
      ENDIF
      IF ( i_am_slave ) THEN
        IF (allocated(rhs_bounds)) THEN
          nb_bytes = nb_bytes -
     &          int(size(rhs_bounds),8)*k34_8
          DEALLOCATE(rhs_bounds)
        ENDIF
        IF (allocated(iwk_solve)) THEN
          nb_bytes = nb_bytes - int(size(iwk_solve),8)*k34_8
          DEALLOCATE( iwk_solve )
        ENDIF
        IF (allocated(ptracb)) THEN
          nb_bytes = nb_bytes - int(size(ptracb),8)*k34_8*
     &                          int(keep(10),8)
          DEALLOCATE( ptracb )
        ENDIF
        IF (allocated(iwcb)) THEN
          nb_bytes = nb_bytes - int(size(iwcb),8)*k34_8
          DEALLOCATE( iwcb )
        ENDIF
C       ------------------------
C       SLAVE CODE
C       -----------------------
C       Deallocate send buffers
C       -----------------------
        IF (id%NSLAVES .GT. 1) THEN
          CALL zmumps_buf_deall_cb( ierr )
          CALL zmumps_buf_deall_small_buf( ierr )
        ENDIF
      END IF
C
      IF ( id%MYID .eq. master ) THEN
C       ------------------------
C       SAVERHS may  have been
C       allocated only on master
C       ------------------------
        IF (allocated(saverhs)) THEN
         nb_bytes = nb_bytes - int(size(saverhs),8)*k35_8
         DEALLOCATE( saverhs)
        ENDIF
C       Nullify RHS_IR might have been pointing to id%RHS
        NULLIFY(rhs_ir)
      ELSE
C       --------------------
C       Free right-hand-side
C       on slave processors
C       --------------------
        IF (associated(rhs_ir)) THEN
          nb_bytes = nb_bytes - int(size(rhs_ir),8)*k35_8
          DEALLOCATE(rhs_ir)
          NULLIFY(rhs_ir)
        END IF
      END IF
      IF (i_am_slave) THEN
C       Deallocate temporary workspace SRW3
        IF (allocated(srw3)) THEN
          nb_bytes = nb_bytes - int(size(srw3),8)*k35_8
          DEALLOCATE(srw3)
        ENDIF
        IF (lscal .AND. icntl21==1) THEN
C         Free local scaling arrays
          nb_bytes = nb_bytes -
     &              int(size(scaling_data_sol%SCALING_LOC),8)*k16_8
          DEALLOCATE(scaling_data_sol%SCALING_LOC)
          NULLIFY(scaling_data_sol%SCALING_LOC)
        ENDIF
C       Free memory until next call to ZMUMPS
        IF (wk_user_provided) THEN
C         S points to WK_USER provided by user
C         KEEP8(24) holds size of WK_USER
C         it should be saved and is used
C         in incore to check that size provided is consistent
C         (see error -41)
          NULLIFY(id%S)
        ELSE IF (associated(id%S).AND.keep(201).GT.0) THEN
C         OOC: free space for S that was allocated
          nb_bytes = nb_bytes - keep8(23)*k35_8
          id%KEEP8(23)=0_8
          DEALLOCATE(id%S)
          NULLIFY(id%S)
        ENDIF
        IF (keep(221).NE.1) THEN
C       -- After reduction of RHS to Schur variables
C       -- keep compressed RHS generated during FWD step
C       -- to be used for future expansion
         IF (associated(id%RHSCOMP)) THEN
            nb_bytes = nb_bytes - id%KEEP8(25)*k35_8
            DEALLOCATE(id%RHSCOMP)
            NULLIFY(id%RHSCOMP)
            id%KEEP8(25)=0_8
         ENDIF
         IF (associated(id%POSINRHSCOMP_ROW)) THEN
            nb_bytes = nb_bytes -
     &                 int(size(id%POSINRHSCOMP_ROW),8)*k34_8
            DEALLOCATE(id%POSINRHSCOMP_ROW)
            NULLIFY(id%POSINRHSCOMP_ROW)
         ENDIF
         IF (id%POSINRHSCOMP_COL_ALLOC) THEN
            nb_bytes = nb_bytes -
     &                 int(size(id%POSINRHSCOMP_COL),8)*k34_8
            DEALLOCATE(id%POSINRHSCOMP_COL)
            NULLIFY(id%POSINRHSCOMP_COL)
            id%POSINRHSCOMP_COL_ALLOC = .false.
         ENDIF
        ENDIF
        IF ( work_wcb_allocated ) THEN
          nb_bytes = nb_bytes - int(size(work_wcb),8)*k35_8
          DEALLOCATE( work_wcb )
        ENDIF
C       Otherwise, WORK_WCB may point to some
C       position inside id%S, nullify it
        NULLIFY( work_wcb )
      ENDIF
      RETURN
 55   FORMAT (//' ERROR ANALYSIS BEFORE ITERATIVE REFINEMENT')
 100  FORMAT(//' ****** SOLVE & CHECK STEP ********'/)
 110  FORMAT (//' Vector solution for column ',i12)
 115  FORMAT(1x, a44,1p,d9.2)
 434  FORMAT(' Time to build/scatter RHS        =',f15.6)
 432  FORMAT(' Time in solution step (fwd/bwd)  =',f15.6)
 435  FORMAT('  .. Time in forward (fwd) step   =   ',f15.6)
 437  FORMAT('  .. Time in ScaLAPACK root       =   ',f15.6)
 436  FORMAT('  .. Time in backward (bwd) step  =   ',f15.6)
 433  FORMAT(' Time to gather solution(cent.sol)=',f15.6)
 431  FORMAT(' Time to copy/scale dist. solution=',f15.6)
 150  FORMAT(' GLOBAL STATISTICS PRIOR SOLVE PHASE ...........'/
     &        ' Number of right-hand-sides                    =',i12/
     &        ' Blocking factor for multiple rhs              =',i12/
     &        ' ICNTL (9)                                     =',i12/
     &        '  --- (10)                                     =',i12/
     &        '  --- (11)                                     =',i12/
     &        '  --- (20)                                     =',i12/
     &        '  --- (21)                                     =',i12/
     &        '  --- (30)                                     =',i12/
     &        '  --- (35)                                     =',i12
     &        )
 151  FORMAT ('  --- (25)                                     =',i12)
 152  FORMAT ('  --- (26)                                     =',i12)
 153  FORMAT ('  --- (32)                                     =',i12)
 160  FORMAT (' RHS'/(1x,1p,5d14.6))
 170  FORMAT (/' ERROR ANALYSIS' )
 240  FORMAT (1x, a42,i4)
 270  FORMAT (//' BEGIN ITERATIVE REFINEMENT' )
  81  FORMAT (/' STATISTICS AFTER ITERATIVE REFINEMENT ')
 131  FORMAT (/' END   ITERATIVE REFINEMENT ')
 141  FORMAT(1x, a52,i4)
      CONTAINS
        SUBROUTINE zmumps_check_distrhs(
     &       idNloc_RHS,
     &       idLRHS_loc,
     &       NRHS,
     &       idIRHS_loc,
     &       idRHS_loc,
     &       INFO)
C
C  Purpose:
C  =======
C
C       Check distributed RHS format. We assume that
C       the user has indicated that he/she provided
C       a distributed RHS (KEEP(248)=-1). We also
C       assume that the nb of RHS columns NRHS has
C       been broadcasted to all processes. This
C       routine should then be called on the workers.
C
C  Arguments:
C  =========
C
        INTEGER, INTENT( IN ) :: idNloc_RHS
        INTEGER, INTENT( IN ) :: idLRHS_loc
        INTEGER, INTENT( IN ) :: NRHS
#if defined(MUMPS_F2003)
        INTEGER, INTENT( IN ), POINTER :: idIRHS_loc (:)
        COMPLEX(kind=8), INTENT( IN ), POINTER :: idRHS_loc  (:)
#else
        INTEGER, POINTER :: idIRHS_loc (:)
        COMPLEX(kind=8), POINTER :: idRHS_loc  (:)
#endif
        INTEGER, INTENT( INOUT ) :: INFO(80)
C
C  Local declarations:
C  ==================
C
        INTEGER(8) :: REQSIZE8
C
C  Executable statements:
C  =====================
C
C       Quick return if nothing on this proc
        IF (idnloc_rhs .LE. 0) RETURN
C       Check for leading dimension
        IF (nrhs.NE.1) THEN
          IF ( idlrhs_loc .LT. idnloc_rhs) THEN
            info(1)=-55
            info(2)=idlrhs_loc
            RETURN
          ENDIF
        ENDIF
        IF (idnloc_rhs .GT. 0) THEN
C         Check association and size of index array idIRHS_loc
          IF (.NOT. associated(idirhs_loc)) THEN
            id%INFO(1)=-22
            id%INFO(2)=17
            RETURN
          ELSE IF (size(idirhs_loc) .LT. idnloc_rhs) THEN
            info(1)=-22
            info(2)= 17
            RETURN
          ENDIF
C         Check association and size of value array idRHS_loc
          IF (.NOT. associated(idrhs_loc)) THEN
            id%INFO(1)=-22
            id%INFO(2)=18
            RETURN
          ELSE
C           Check size of array of values idRHS_loc
            reqsize8 = int(idlrhs_loc,8)*int(nrhs,8)
     &                  + int(-idlrhs_loc+idnloc_rhs,8)
#if defined(MUMPS_F2003)
            IF (size(idrhs_loc,kind=8) .LT. reqsize8) THEN
#else
            IF ( reqsize8 .LE. int(huge(idnloc_rhs),8) .AND.
     &        size(idrhs_loc) .LT. int(reqsize8) ) THEN
C             (Warning: this assumes that size(idRHS_loc)
C             does not overflow)
#endif
              info(1)=-22
              info(2)=18
              RETURN
            ENDIF
          ENDIF
        ENDIF
        RETURN
        END SUBROUTINE zmumps_check_distrhs
        SUBROUTINE zmumps_pp_solve()
        IMPLICIT NONE
C
C       Purpose:
C       =======
C       Scatter right-hand side, solve the system,
C       and gather the solution on the host during
C       post-processing.
C       We use an internal subroutine to avoid code
C       duplication without the complication of adding
C       new parameters or local variables. All variables
C       in this routine have the scope of ZMUMPS_SOL_DRIVER.
C
C
        IF (kase .NE. 1 .AND. kase .NE. 2) THEN
          WRITE(*,*) "Internal error 1 in ZMUMPS_PP_SOLVE"
          CALL mumps_abort()
        ENDIF
        IF ( id%MYID .eq. master ) THEN
C         Define matrix B as follows:
C            MTYPE=1 => B=A other values B=At
C         The user asked to solve the system Bx=b
C
C         THEN
C           KASE = 1........ RW1 = INV(TRANSPOSE(B)) * RW1
C           KASE = 2........ RW1 = INV(B) * RW1
          IF ( mtype .EQ. 1 ) THEN
            solvet = kase - 1
          ELSE
            solvet = kase
          END IF
C         SOLVET= 1 -> solve A x = B, other values solve Atx=b
C         We force SOLVET to have value either 0 or 1, in order
C         to be able to test both values, and also, be able to
C         test whether SOLVET = MTYPE or not.
          IF ( solvet.EQ.2 ) solvet = 0
          IF ( lscal ) THEN
            IF ( solvet .EQ. 1 ) THEN
C             Apply rowscaling
              DO k = 1, id%N
                c_y( k ) = c_y( k ) * id%ROWSCA( k )
              END DO
            ELSE
C             Apply column scaling
              DO k = 1, id%N
                c_y( k ) = c_y( k ) * id%COLSCA( k )
              END DO
            END IF
          END IF
        END IF ! MYID.EQ.MASTER
C       ------------------------------
C       Broadcast SOLVET to the slaves
C       ------------------------------
        CALL mpi_bcast( solvet, 1, mpi_integer, master,
     &                  id%COMM, ierr)
C       --------------------------------------------
C       Scatter the right hand side C_Y on all procs
C       --------------------------------------------
        IF ( .NOT.i_am_slave ) THEN
C         -- Master not working
          CALL zmumps_scatter_rhs(id%NSLAVES,id%N, id%MYID,
     &      id%COMM,
     &      solvet, c_y(1), id%N, 1,
     &      1,
     &      c_dummy, 1, 1,
     &      idummy, 0,
     &      jdummy, id%KEEP(1), id%KEEP8(1), id%PROCNODE_STEPS(1),
     &      idummy, 1,
     &      id%STEP(1),
     &      id%ICNTL(1),id%INFO(1))
        ELSE
          IF (solvet.EQ.mtype) THEN
C           POSINRHSCOMP_ROW is with respect to the
C           original linear system (transposed or not)
            ptr_posinrhscomp_fwd => id%POSINRHSCOMP_ROW
          ELSE
C           Transposed, use column indices of original
C           system (ie, col indices of A or A^T)
            ptr_posinrhscomp_fwd => id%POSINRHSCOMP_COL
          ENDIF
          liw_passed = max( liw, 1 )
          CALL zmumps_scatter_rhs(id%NSLAVES,id%N, id%MYID,
     &      id%COMM,
     &      solvet,  c_y(1), id%N, 1,
     &      1,
     &      id%RHSCOMP(ibeg_rhscomp), ld_rhscomp, 1,
     &      ptr_posinrhscomp_fwd(1), nb_fs_rhscomp_f,
C
     &      id%PTLUST_S(1), id%KEEP(1), id%KEEP8(1),
     &      id%PROCNODE_STEPS(1),
     &      is(1), liw_passed,
     &      id%STEP(1),
     &      id%ICNTL(1),id%INFO(1))
        ENDIF
        IF (info(1).LT.0) GOTO 89
C
C       Solve the system
C
        IF ( i_am_slave ) THEN
          liw_passed = max( liw, 1 )
          la_passed = max( la, 1_8 )
          IF (solvet.EQ.mtype) THEN
            ptr_posinrhscomp_fwd => id%POSINRHSCOMP_ROW
            ptr_posinrhscomp_bwd => id%POSINRHSCOMP_COL
          ELSE
            ptr_posinrhscomp_fwd => id%POSINRHSCOMP_COL
            ptr_posinrhscomp_bwd => id%POSINRHSCOMP_ROW
          ENDIF
          from_pp=.true.
          nbsparse_loc = .false.
          CALL zmumps_sol_c(id%root, id%N, id%S(1), la_passed, id%IS(1),
     & liw_passed,work_wcb(1),lwcb8_sol_c,iwcb,liwcb,nbrhs_eff,id%NA(1),
     & id%LNA,id%NE_STEPS(1),srw3,solvet,icntl(1),from_pp,id%STEP(1),
     & id%FRERE_STEPS(1),id%DAD_STEPS(1),id%FILS(1),id%PTLUST_S(1),
     & id%PTRFAC(1), iwk_solve(1), liwk_solve, ptracb, liwk_ptracb,
     & id%PROCNODE_STEPS(1), id%NSLAVES, info(1), keep(1), keep8(1),
     & id%DKEEP(1),id%COMM_NODES,id%MYID,id%MYID_NODES, bufr(1), lbufr,
     & lbufr_bytes, id%ISTEP_TO_INIV2(1), id%TAB_POS_IN_PERE(1,1),
C      Next 3 arguments are not used in this call
     & ibeg_root_def,iend_root_def,iroot_def_rhs_col1, ptr_rhs_root(1),
     & lptr_rhs_root, size_root, master_root, id%RHSCOMP(ibeg_rhscomp),
     & ld_rhscomp,ptr_posinrhscomp_fwd(1),ptr_posinrhscomp_bwd(1),
     & 1,1,1,1, idummy, 1, jdummy, kdummy, 1, ldummy, 1, mdummy, 1,1,
     & nbsparse_loc, ptr_rhs_bounds(1), lptr_rhs_bounds
     & , 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%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%L0_OMP_FACTORS(1), id%LL0_OMP_FACTORS
     &    )
        END IF
C       ------------------
C       Change error codes
C       ------------------
        IF (info(1).eq.-2) info(1)=-12
        IF (info(1).eq.-3) info(1)=-15
C
        IF (info(1) .GE. 0) THEN
C         We need a workspace of minimal size KEEP(247)
C         in order to unpack pieces of the solution during
C         ZMUMPS_GATHER_SOLUTION below
C         - Avoid allocation if error already occurred.
C         - DEALLOCATE called after GATHER_SOLUTION
C         CWORK not needed for AM1
          ALLOCATE( cwork(max(max(keep(247),keep(246)),1)),
     &           stat=allocok)
          IF (allocok > 0) THEN
            info(1)=-13
            info(2)=max(max(keep(247),keep(246)),1)
          ENDIF
        ENDIF
C       -------------------------
C       Propagate possible errors
C       -------------------------
 89     CALL mumps_propinfo( icntl(1), info(1),
     &                   id%COMM,id%MYID)
C
C       Return in case of error.
        IF (info(1).LT.0) RETURN
C       -------------------------------
C       Assemble the solution on master
C       -------------------------------
C       (Note: currently, if this part of code is executed,
C       then necessarily NBRHS_EFF = 1)
C
C       === GATHER and SCALE solution ==============
C
        IF ((id%MYID.NE.master).OR. .NOT.lscal) THEN
          pt_scaling => dummy_scal
        ELSE
          IF (solvet.EQ.1) THEN
            pt_scaling => id%COLSCA
          ELSE
            pt_scaling => id%ROWSCA
          ENDIF
        ENDIF
        liw_passed = max( liw, 1 )
C       Solution computed during ZMUMPS_SOL_C has been stored
C       in id%RHSCOMP and is gathered on the master in C_Y
        IF ( .NOT. i_am_slave ) THEN
C         I did not participate to computing part of the solution
C         (id%RHSCOMP not set/allocate) : receive solution, store
C         it and scale it.
          CALL zmumps_gather_solution(id%NSLAVES,id%N,
     &      id%MYID, id%COMM, nbrhs_eff,
     &      solvet, c_y, id%N, nbrhs_eff, 1,
     &      jdummy, id%KEEP(1),id%KEEP8(1), id%PROCNODE_STEPS(1),
     &      idummy, 1,
     &      id%STEP(1), bufr(1), lbufr, lbufr_bytes,
     &      cwork(1), size(cwork),
     &      lscal, pt_scaling(1), size(pt_scaling),
!     RHSCOMP not on non-working master
     &      c_dummy, 1 , 1, idummy, 1,
!     for sparse permuted RHS on host
     &      perm_rhs, size(perm_rhs)
     &      )
        ELSE
          CALL zmumps_gather_solution(id%NSLAVES,id%N,
     &      id%MYID, id%COMM, nbrhs_eff,
     &      solvet, c_y, id%N, nbrhs_eff, 1,
     &      id%PTLUST_S(1), id%KEEP(1),id%KEEP8(1),
     &      id%PROCNODE_STEPS(1),
     &      is(1), liw_passed,
     &      id%STEP(1), bufr(1), lbufr, lbufr_bytes,
     &      cwork(1), size(cwork),
     &      lscal, pt_scaling(1), size(pt_scaling),
     &      id%RHSCOMP(ibeg_rhscomp), ld_rhscomp, nbrhs_eff,
     &      ptr_posinrhscomp_bwd(1), id%N,
     &      perm_rhs, size(perm_rhs)) ! for sparse permuted RHS on host
        ENDIF
        DEALLOCATE( cwork )
        END SUBROUTINE zmumps_pp_solve
      END SUBROUTINE zmumps_solve_driver