stype3_root.F

Go to the documentation of this file.

C
C  This file is part of MUMPS 5.5.1, released
C  on Tue Jul 12 13:17:24 UTC 2022
C
C
C  Copyright 1991-2022 CERFACS, CNRS, ENS Lyon, INP Toulouse, Inria,
C  Mumps Technologies, University of Bordeaux.
C
C  This version of MUMPS is provided to you free of charge. It is
C  released under the CeCILL-C license 
C  (see doc/CeCILL-C_V1-en.txt, doc/CeCILL-C_V1-fr.txt, and
C  https://cecill.info/licences/Licence_CeCILL-C_V1-en.html)
C
      SUBROUTINE smumps_ass_root( root, KEEP50,
     &                     NROW_SON, NCOL_SON, INDROW_SON, 
     &                     INDCOL_SON, NSUPCOL, VAL_SON, VAL_ROOT,
     &                     LOCAL_M, LOCAL_N,
     &                     RHS_ROOT, NLOC_ROOT, CBP )
      USE smumps_struc_def, ONLY : smumps_root_struc
      IMPLICIT NONE
      TYPE (SMUMPS_ROOT_STRUC) :: root
      INTEGER, INTENT(IN) :: KEEP50
      INTEGER NCOL_SON, NROW_SON, NSUPCOL
      INTEGER, intent(in) :: CBP   
      INTEGER INDROW_SON( NROW_SON ), INDCOL_SON( NCOL_SON )
      INTEGER LOCAL_M, LOCAL_N
      REAL VAL_SON( NCOL_SON, NROW_SON )
      REAL VAL_ROOT( LOCAL_M, LOCAL_N )
      INTEGER NLOC_ROOT
      REAL RHS_ROOT( LOCAL_M, NLOC_ROOT )
      INTEGER I, J, INDROW, INDCOL, IPOSROOT, JPOSROOT
      IF (cbp .EQ. 0) THEN
        DO i = 1, nrow_son
          indrow = indrow_son(i)
          iposroot = (root%NPROW*((indrow-1)/root%MBLOCK)+root%MYROW)
     &             * root%MBLOCK + mod(indrow-1,root%MBLOCK) + 1
          DO j = 1, ncol_son-nsupcol
          indcol = indcol_son(j)
          IF (keep50.NE.0) THEN
            jposroot = (root%NPCOL*((indcol-1)/root%NBLOCK)+root%MYCOL)
     &               * root%NBLOCK + mod(indcol-1,root%NBLOCK) + 1
            IF (iposroot < jposroot) THEN
              cycle
            ENDIF
          ENDIF
          val_root( indrow, indcol ) =
     &    val_root( indrow, indcol ) + val_son(j,i)
          END DO
          DO j = ncol_son-nsupcol+1, ncol_son
            indcol = indcol_son(j)
            rhs_root( indrow, indcol ) =
     &      rhs_root( indrow, indcol ) + val_son(j,i)
          ENDDO
        END DO
      ELSE
        DO i=1, nrow_son  
          DO j = 1, ncol_son  
           rhs_root( indrow_son( i ), indcol_son(j)) =
     &     rhs_root(indrow_son(i),indcol_son(j)) + val_son(j,i)
          ENDDO
        ENDDO
      ENDIF
      RETURN
      END SUBROUTINE smumps_ass_root
      RECURSIVE SUBROUTINE smumps_build_and_send_cb_root
     &  ( comm_load, ass_irecv, n, ison, iroot,
     &    ptri, ptrr,
     &    root, nbrow, nbcol, shift_list_row_son,
     &    shift_list_col_son, shift_val_son_arg, lda_arg, tag,
     &    myid, comm, bufr, lbufr, lbufr_bytes, procnode_steps, posfac,
     &    iwpos, iwposcb, iptrlu, lrlu, lrlus, iw, liw, a, la,
     &    ptrist, ptlust_s, ptrfac,
     &    ptrast, step, pimaster, pamaster,
     &    nstk, comp, iflag, ierror, perm,
     &    ipool, lpool, leaf, nbfin, slavef,
     &    opassw, opeliw, itloc, rhs_mumps,
     &    fils, dad, ptrarw, ptraiw,
     &    intarr,dblarr,icntl,keep,keep8,dkeep,transpose_asm,
     &    nd, frere,
     &    lptrar, nelt, frtptr, frtelt, 
     &    istep_to_iniv2, tab_pos_in_pere  
     &               , lrgroups
     &     )
      USE smumps_ooc        
      USE smumps_buf
      USE smumps_load
      USE smumps_struc_def, ONLY : smumps_root_struc
      USE smumps_dynamic_memory_m, ONLY : smumps_dm_set_dynptr
      IMPLICIT NONE
      INTEGER keep(500), icntl(60)
      INTEGER(8) keep8(150)
      REAL  dkeep(230)
      TYPE (smumps_root_struc) :: root
      INTEGER comm_load, ass_irecv
      INTEGER n, ison, iroot, tag
      INTEGER ptri( keep(28) )
      INTEGER(8) :: ptrr( keep(28) )
      INTEGER nbrow, nbcol
      INTEGER, INTENT(IN):: lda_arg
      INTEGER(8), INTENT(IN) :: shift_val_son_arg
      INTEGER shift_list_row_son, shift_list_col_son
      INTEGER myid, comm
      LOGICAL transpose_asm
      include 'mpif.h'
      INTEGER lbufr, lbufr_bytes
      INTEGER bufr( lbufr )
      INTEGER(8) :: posfac, iptrlu, lrlu, lrlus, la
      INTEGER iwpos, iwposcb
      INTEGER liw
      INTEGER IW( liw )
      REAL a( la )
      INTEGER, intent(in) :: lrgroups(n)
      INTEGER lptrar, nelt
      INTEGER frtptr( n+1 ), frtelt( nelt )
      INTEGER(8) :: ptrast(keep(28))
      INTEGER(8) :: ptrfac(keep(28))
      INTEGER(8) :: pamaster(keep(28))
      INTEGER ptrist( keep(28) ), ptlust_s(keep(28))
      INTEGER step(n), pimaster(keep(28)), nstk( n )
      INTEGER comp, iflag, ierror
      INTEGER PERM(n)
      INTEGER lpool, LEAF
      INTEGER ipool( lpool )
      INTEGER nbfin, slavef
      DOUBLE PRECISION opassw, opeliw
      INTEGER procnode_steps( keep(28) )
      INTEGER itloc( n + keep(253) ), fils( n ), dad(keep(28))
      REAL :: RHS_MUMPS(keep(255))
      INTEGER nd( keep(28) ), FRERE( keep(28) )
      INTEGER(8), INTENT(IN) :: ptrarw( lptrar ), ptraiw( lptrar )
      INTEGER intarr( keep8(27) )
      REAL dblarr( keep8(26) )
      INTEGER istep_to_iniv2(keep(71)), 
     &        tab_pos_in_pere(slavef+2,max(1,keep(56)))
      REAL, DIMENSION(:), POINTER :: sona_ptr
      INTEGER(8) :: lsona_ptr, possona_ptr
      INTEGER allocok
      INTEGER, ALLOCATABLE, DIMENSION(:) :: ptrrow,  ptrcol
      INTEGER, ALLOCATABLE, DIMENSION(:) :: nsuprow, nsupcol
      INTEGER, ALLOCATABLE, DIMENSION(:) :: row_index_list
      INTEGER, ALLOCATABLE, DIMENSION(:) :: col_index_list
      INTEGER :: status(mpi_status_size)
      INTEGER i, pos_in_root, irow, jcol, iglob, jglob
      INTEGER pdest, ierr
      INTEGER local_m, local_n
      INTEGER(8) :: posroot
      INTEGER nsubset_row, NSUBSET_COL
      INTEGER nrlocal, nclocal
      INTEGER :: lda
      INTEGER(8) :: shift_val_son
      LOGICAL set_irecv, blocking, message_received
      INTEGER nbrows_already_sent
      INTEGER size_msg
      INTEGER lp
      include 'mumps_headers.h'
      LOGICAL skiplast_rhs_rows, bcp_sym_nonempty
      INTEGER BBPCBP
      INTEGER mumps_procnode
      EXTERNAL mumps_procnode
      bbpcbp  = 0   
      lp = icntl(1)
      IF ( icntl(4) .LE. 0 ) lp = -1
      IF (lda_arg < 0) THEN 
        CALL smumps_set_lda_shift_val_son(iw, liw, ptri(step(ison)),
     &                             lda, shift_val_son)
      ELSE
        lda = lda_arg
        shift_val_son = shift_val_son_arg
      ENDIF
      ALLOCATE(ptrrow(root%NPROW + 1 ),  stat=allocok)
      if (allocok .GT. 0) THEN
       iflag  =-13
       ierror = root%NPROW + 1
      endif
      ALLOCATE(ptrcol(root%NPCOL + 1 ),  stat=allocok)
      if (allocok .GT. 0) THEN
       iflag  =-13
       ierror = root%NPCOL + 1
      endif
      ALLOCATE(nsuprow(root%NPROW + 1 ),  stat=allocok)
      if (allocok .GT. 0) THEN
       iflag  =-13
       ierror = root%NPROW + 1
      endif
      ALLOCATE(nsupcol(root%NPCOL + 1 ),  stat=allocok)
      if (allocok .GT. 0) THEN
       iflag  =-13
       ierror = root%NPCOL + 1
      endif
      IF (iflag.LT.0) THEN
         IF (lp > 0) write(6,*) myid, ' : MEMORY ALLOCATION ',
     &     'FAILURE in SMUMPS_BUILD_AND_SEND_CB_ROOT'
         CALL smumps_bdc_error( myid, slavef, comm, keep )
         RETURN
      ENDIF
      skiplast_rhs_rows = ((keep(253).GT.0).AND.(keep(50).EQ.0))
      bcp_sym_nonempty = .false.
      ptrrow = 0
      ptrcol = 0
      nsuprow = 0
      nsupcol = 0
      DO i = 1, nbrow                   
        iglob  =  iw( ptri(step(ison)) +
     &                          shift_list_row_son + i - 1 )
        IF (skiplast_rhs_rows.AND.(iglob.GT.n)) cycle
        IF ( .NOT. transpose_asm ) THEN 
          IF (iglob.GT.n) THEN
            bcp_sym_nonempty = .true.
            pos_in_root = iglob - n
            jcol =  mod((pos_in_root-1)/root%NBLOCK,root%NPCOL) 
            nsupcol(jcol+1) = nsupcol(jcol+1) + 1
            ptrcol( jcol + 2 ) = ptrcol( jcol + 2 ) + 1
          ELSE  
            pos_in_root = root%RG2L_ROW( iglob ) 
            irow  = mod((pos_in_root-1)/root%MBLOCK,root%NPROW)
            ptrrow( irow + 2 ) = ptrrow( irow + 2 ) + 1
          ENDIF
        ELSE        
          IF (iglob .GT. n) THEN 
            pos_in_root = iglob - n
          ELSE  
            pos_in_root = root%RG2L_COL( iglob ) 
          ENDIF
          jcol =  mod( ( pos_in_root - 1 ) / root%NBLOCK, root%NPCOL ) 
          IF (iglob.GT.n)  
     &               nsupcol(jcol+1) = nsupcol(jcol+1) + 1
          ptrcol( jcol + 2 ) = ptrcol( jcol + 2 ) + 1
        END IF
      END DO
      IF (keep(50).NE.0 .AND.(.NOT.transpose_asm).AND.bcp_sym_nonempty)
     &             bbpcbp = 1
      DO i = 1, nbcol                   
        jglob   =  iw( ptri(step(ison)) +
     &                shift_list_col_son + i - 1 ) 
        IF ((keep(50).GT.0) .AND. (jglob.GT.n)) cycle  
        IF ( .NOT. transpose_asm ) THEN
          IF (keep(50).EQ.0) THEN
            IF (jglob.LE.n) THEN
              pos_in_root = root%RG2L_COL(jglob)  
            ELSE
              pos_in_root = jglob-n
            ENDIF
            jcol =  mod((pos_in_root-1) / root%NBLOCK, root%NPCOL ) 
            IF (jglob.GT.n) THEN
             nsupcol(jcol+1) = nsupcol(jcol+1)  + 1  
            ENDIF
            ptrcol( jcol + 2 ) = ptrcol( jcol + 2 ) + 1
          ELSE 
            pos_in_root = root%RG2L_COL(jglob) 
            jcol =  mod((pos_in_root-1) / root%NBLOCK, root%NPCOL )
            ptrcol( jcol + 2 ) = ptrcol( jcol + 2 ) + 1
            IF (bcp_sym_nonempty) THEN
             pos_in_root = root%RG2L_ROW(jglob) 
             irow  = mod((pos_in_root-1)/root%MBLOCK,root%NPROW)
             nsuprow(irow+1) = nsuprow(irow+1)+1
             ptrrow( irow + 2 ) = ptrrow( irow + 2 ) + 1
            ENDIF
          ENDIF
        ELSE  
          IF (jglob.LE.n) THEN
           pos_in_root = root%RG2L_ROW( jglob ) 
          ELSE
           pos_in_root = jglob-n
          ENDIF
          irow        = mod( ( pos_in_root - 1 ) /
     &                  root%MBLOCK, root%NPROW )
          ptrrow( irow + 2 ) = ptrrow( irow + 2 ) + 1
        END IF
      END DO
      ptrrow( 1 ) = 1
      DO irow = 2, root%NPROW + 1
        ptrrow( irow ) = ptrrow( irow ) + ptrrow( irow - 1 )
      END DO
      ptrcol( 1 ) = 1
      DO jcol = 2, root%NPCOL + 1
        ptrcol( jcol ) = ptrcol( jcol ) + ptrcol( jcol - 1 )
      END DO
      ALLOCATE(row_index_list(ptrrow(root%NPROW+1)-1+1),
     &         stat=allocok)
      if (allocok .GT. 0) THEN
       iflag  =-13
       ierror = ptrrow(root%NPROW+1)-1+1
      endif
      ALLOCATE(col_index_list(ptrcol(root%NPCOL+1)-1+1),
     &         stat=allocok)
      if (allocok .GT. 0) THEN
       iflag  =-13
       ierror = ptrcol(root%NPCOL+1)-1+1
      endif
      DO i = 1, nbrow
        iglob  =  iw( ptri(step(ison)) +
     &                          shift_list_row_son + i - 1 )
        IF (skiplast_rhs_rows.AND.(iglob.GT.n)) cycle
        IF ( .NOT. transpose_asm ) THEN
          IF (iglob.GT.n) cycle   
          pos_in_root = root%RG2L_ROW( iglob )
          irow        = mod( ( pos_in_root - 1 ) / root%MBLOCK,
     &                       root%NPROW )
          row_index_list( ptrrow( irow + 1 ) ) = i 
          ptrrow( irow + 1 ) = ptrrow( irow + 1 ) + 1
        ELSE
          IF (iglob.LE.n) THEN
           pos_in_root = root%RG2L_COL( iglob )
          ELSE
           pos_in_root = iglob - n  
          ENDIF
          jcol        = mod( ( pos_in_root - 1 ) / root%NBLOCK,
     &                       root%NPCOL )
          col_index_list( ptrcol( jcol + 1 ) ) = i 
          ptrcol( jcol + 1 ) = ptrcol( jcol + 1 ) + 1 
        END IF
      END DO
      DO i = 1, nbcol 
        jglob =  iw( ptri(step(ison))+shift_list_col_son+i - 1 ) 
        IF ((keep(50).GT.0) .AND. (jglob.GT.n)) cycle  
        IF ( .NOT. transpose_asm ) THEN
          IF ( jglob.LE.n ) THEN
           pos_in_root = root%RG2L_COL( jglob )
          ELSE
           pos_in_root = jglob - n
          ENDIF
          jcol        = mod( ( pos_in_root - 1 ) /
     &               root%NBLOCK, root%NPCOL )
          col_index_list( ptrcol( jcol + 1 ) ) = i 
          ptrcol( jcol + 1 ) = ptrcol( jcol + 1 ) + 1
        ELSE
          IF ( jglob.LE.n ) THEN
           pos_in_root = root%RG2L_ROW( jglob )
          ELSE
           pos_in_root = jglob - n
          ENDIF
          irow        = mod( ( pos_in_root - 1 ) /
     &                root%MBLOCK, root%NPROW )
          row_index_list( ptrrow( irow + 1 ) ) = i    
          ptrrow( irow + 1 ) = ptrrow( irow + 1 ) + 1 
        END IF
      END DO
      IF (bcp_sym_nonempty) THEN
        DO i = 1, nbrow
          iglob  =  iw( ptri(step(ison)) +
     &                         shift_list_row_son + i - 1 )
          IF (iglob.LE.n) cycle  
          pos_in_root = iglob - n
          jcol =  mod((pos_in_root-1)/root%NBLOCK,root%NPCOL)
          col_index_list( ptrcol( jcol + 1 ) ) = i 
          ptrcol( jcol + 1 ) = ptrcol( jcol + 1 ) + 1
        ENDDO
        DO i=1, nbcol
         jglob =  iw( ptri(step(ison))+shift_list_col_son+i - 1 ) 
         IF (jglob.GT.n) THEN
           EXIT
         ELSE
           pos_in_root = root%RG2L_ROW(jglob) 
         ENDIF
         irow  = mod((pos_in_root-1)/root%MBLOCK,root%NPROW)
         row_index_list( ptrrow( irow + 1 ) ) = i    
         ptrrow( irow + 1 ) = ptrrow( irow + 1 ) + 1 
        ENDDO
      ENDIF
      DO irow = root%NPROW, 2, -1
        ptrrow( irow ) = ptrrow( irow - 1 )
      END DO
      ptrrow( 1 ) = 1
      DO jcol = root%NPCOL, 2, -1
        ptrcol( jcol ) = ptrcol( jcol - 1 )
      END DO
      ptrcol( 1 ) = 1
      jcol  = root%MYCOL
      irow  = root%MYROW
      IF ( root%yes ) THEN
         if (irow .ne. root%MYROW .or. jcol.ne.root%MYCOL) then
        write(*,*) ' error in grid position buildandsendcbroot'
        CALL mumps_abort()
        end if
        IF ( ptrist(step(iroot)).EQ.0.AND.
     &       ptlust_s(step(iroot)).EQ.0) THEN
           CALL smumps_root_alloc_static(root, iroot, n, iw, liw,
     &               a, la,
     &               fils, dad, myid, slavef, procnode_steps,
     &               lptrar, nelt, frtptr, frtelt,
     &               ptraiw, ptrarw, intarr, dblarr,
     &               lrlu, iptrlu,
     &               iwpos, iwposcb, ptrist, ptrast,
     &               step, pimaster, pamaster, itloc, rhs_mumps,
     &               comp, lrlus, iflag, keep,keep8,dkeep, ierror )
           keep(121) = -1 
           IF (iflag.LT.0) THEN
                CALL smumps_bdc_error( myid, slavef, comm, keep )
                RETURN
           ENDIF
        ELSE
           keep(121) = keep(121) - 1
           IF ( keep(121) .eq. 0 ) THEN 
              IF (keep(201).EQ.1) THEN
                 CALL smumps_ooc_force_wrt_buf_panel(ierr)
              ELSE IF (keep(201).EQ.2) THEN
                 CALL smumps_force_write_buf(ierr)              
              ENDIF
              CALL smumps_insert_pool_n(n, ipool, lpool, procnode_steps,
     &        slavef, keep(199), keep(28), keep(76), keep(80), keep(47),
     &        step, iroot+n )
              IF (keep(47) .GE. 3) THEN
                 CALL smumps_load_pool_upd_new_pool(
     &                ipool, lpool, 
     &                procnode_steps, keep,keep8, slavef, comm_load,
     &                myid, step, n, nd, fils )
              ENDIF
          END IF
        END IF
        CALL smumps_dm_set_dynptr( iw(ptri(step(ison))+xxs), a, la,
     &      ptrr(step(ison)), iw(ptri(step(ison))+xxd),
     &      iw(ptri(step(ison))+xxr),
     &      sona_ptr, possona_ptr, lsona_ptr )
       IF (keep(60) .NE. 0 ) THEN
         local_m = root%SCHUR_LLD
         local_n = root%SCHUR_NLOC
            nrlocal = ptrrow( irow + 2 ) - ptrrow( irow + 1 )
            nclocal = ptrcol( jcol + 2 ) - ptrcol( jcol + 1 )
            CALL smumps_root_local_assembly( n,
     &        root%SCHUR_POINTER(1),
     &        local_m, local_n,
     &        root%NPCOL, root%NPROW, root%MBLOCK, root%NBLOCK,
     &        nbcol, nbrow,
     &        iw( ptri(step(ison)) + shift_list_col_son ),
     &        iw( ptri(step(ison)) + shift_list_row_son ),
     &        lda, sona_ptr( possona_ptr + shift_val_son ),
     &        row_index_list( ptrrow( irow + 1 ) ),
     &        col_index_list( ptrcol( jcol + 1 ) ),
     &        nrlocal,
     &        nclocal,
     &        nsuprow(irow+1), nsupcol(jcol+1),
     &        root%RG2L_ROW(1), root%RG2L_COL(1), transpose_asm,
     &        keep,
     &        root%RHS_ROOT(1,1), root%RHS_NLOC )
       ELSE
        IF ( ptrist(step( iroot )) .GE. 0 ) THEN
          IF ( ptrist(step( iroot )) .EQ. 0 ) THEN
            local_n = iw( ptlust_s(step(iroot)) + 1 + keep(ixsz))
            local_m = iw( ptlust_s(step(iroot)) + 2 + keep(ixsz))
            posroot = ptrfac(iw( ptlust_s(step(iroot)) +4+keep(ixsz) ))
          ELSE
            local_n = - iw( ptrist(step(iroot)) +keep(ixsz))
            local_m = iw( ptrist(step(iroot)) + 1 +keep(ixsz))
            posroot = pamaster(step( iroot ))
          ENDIF
          nclocal = ptrcol( jcol + 2 ) - ptrcol( jcol + 1 )
          nrlocal = ptrrow( irow + 2 ) - ptrrow( irow + 1 )
          CALL smumps_root_local_assembly( n, a( posroot ),
     &        local_m, local_n,
     &        root%NPCOL, root%NPROW, root%MBLOCK, root%NBLOCK,
     &        nbcol, nbrow,
     &        iw( ptri(step(ison)) + shift_list_col_son ),
     &        iw( ptri(step(ison)) + shift_list_row_son ),
     &        lda, sona_ptr( possona_ptr + shift_val_son ),
     &        row_index_list( ptrrow( irow + 1 ) ),
     &        col_index_list( ptrcol( jcol + 1 ) ),
     &        nrlocal,
     &        nclocal,
     &        nsuprow(irow+1), nsupcol(jcol+1),
     &        root%RG2L_ROW(1), root%RG2L_COL(1), transpose_asm,
     &        keep,
     &        root%RHS_ROOT(1,1), root%RHS_NLOC )
        END IF
       ENDIF
      END IF
      DO irow = 0, root%NPROW - 1
        DO jcol = 0, root%NPCOL - 1
          pdest = irow * root%NPCOL + jcol
          IF ( (root%MYROW.eq.irow.and.root%MYCOL.eq.jcol) .and.
     &         myid.ne.pdest) THEN
            write(*,*) 'error: myrow,mycol=',root%MYROW,root%MYCOL
            write(*,*) ' MYID,PDEST=',myid,pdest
            CALL mumps_abort()
          END IF
          IF ( root%MYROW .NE. irow .OR. root%MYCOL .NE. jcol) THEN
            nbrows_already_sent = 0
            ierr = -1
            DO WHILE ( ierr .EQ. -1 )
              nsubset_row = ptrrow( irow + 2 ) - ptrrow( irow + 1 )
              nsubset_col = ptrcol( jcol + 2 ) - ptrcol( jcol + 1 )
              IF ( lrlu .LT. int(nsubset_row,8) * int(nsubset_col,8)
     &        .AND. lrlus .GT. int(nsubset_row,8) * int(nsubset_col,8) )
     &        THEN
                CALL smumps_compre_new(n, keep,
     &          iw, liw, a, la,
     &          lrlu, iptrlu,
     &          iwpos, iwposcb, ptrist, ptrast,
     &          step, pimaster, pamaster, lrlus,
     &          keep(ixsz), comp, dkeep(97),
     &          myid, slavef, procnode_steps, dad)
                IF ( lrlu .NE. lrlus ) THEN
                  WRITE(*,*) myid,": pb compress in",
     &                            "SMUMPS_BUILD_AND_SEND_CB_ROOT"
                  WRITE(*,*) myid,': LRLU, LRLUS=',lrlu,lrlus
                  CALL mumps_abort()
                END IF
              END IF
              CALL smumps_dm_set_dynptr(
     &        iw(ptri(step(ison))+xxs), a, la,
     &        ptrr(step(ison)), iw(ptri(step(ison))+xxd),
     &        iw(ptri(step(ison))+xxr),
     &        sona_ptr, possona_ptr, lsona_ptr )
              CALL smumps_buf_send_contrib_type3_i( n, ison,
     &        nbcol, nbrow,
     &        iw( ptri(step(ison)) + shift_list_col_son ),
     &        iw( ptri(step(ison)) + shift_list_row_son ),
     &        lda, sona_ptr( possona_ptr + shift_val_son ),
     &        tag,
     &        row_index_list( ptrrow( irow + 1 ) ),
     &        col_index_list( ptrcol( jcol + 1 ) ),
     &        nsubset_row, nsubset_col,
     &        nsuprow(irow+1), nsupcol(jcol+1),
     &        root%NPROW, root%NPCOL, root%MBLOCK,
     &        root%RG2L_ROW(1), root%RG2L_COL(1),
     &        root%NBLOCK, pdest,
     &        comm, ierr, a( posfac ), lrlu, transpose_asm,
     &        size_msg, nbrows_already_sent, keep, bbpcbp )
              IF ( ierr .EQ. -1 ) THEN
                  blocking  = .false.
                  set_irecv = .true.
                  message_received = .false.
                  CALL smumps_try_recvtreat( comm_load, ass_irecv, 
     &            blocking, set_irecv, message_received,
     &            mpi_any_source, mpi_any_tag, 
     &            status, bufr, lbufr,
     &            lbufr_bytes, procnode_steps, posfac, iwpos, iwposcb,
     &            iptrlu, lrlu, lrlus, n, iw, liw, a, la,
     &            ptrist, ptlust_s, ptrfac, ptrast, step,
     &            pimaster, pamaster, nstk,
     &            comp, iflag, ierror, comm, perm, ipool, lpool,
     &            leaf, nbfin, myid, slavef, root,
     &            opassw, opeliw, itloc, rhs_mumps, fils, dad,
     &            ptrarw,ptraiw,intarr,dblarr,icntl,keep,keep8,dkeep,
     &            nd, frere, lptrar, nelt, frtptr, frtelt, 
     &            istep_to_iniv2, tab_pos_in_pere, .true.
     &               , lrgroups
     &             )
                  IF ( iflag .LT. 0 ) GOTO 500
                  IF (lda_arg < 0) THEN 
                    CALL smumps_set_lda_shift_val_son(
     &                    iw, liw, ptri(step(ison)),
     &                    lda, shift_val_son)
                  ENDIF
              END IF
            END DO
            IF ( ierr == -2 ) THEN
              iflag  = -17
              ierror = size_msg
              IF (lp > 0) WRITE(lp, *) "FAILURE, SEND BUFFER TOO
     & SMALL DURING SMUMPS_BUILD_AND_SEND_CB_ROOT"
              CALL smumps_bdc_error( myid, slavef, comm, keep )
              GOTO 500
            ENDIF
            IF ( ierr == -3 ) THEN
              IF (lp > 0) WRITE(lp, *) "FAILURE, RECV BUFFER TOO
     & SMALL DURING SMUMPS_BUILD_AND_SEND_CB_ROOT"
              iflag  = -20
              ierror = size_msg
              CALL smumps_bdc_error( myid, slavef, comm, keep )
              GOTO 500
            ENDIF
          END IF
        END DO
      END DO
 500  CONTINUE
      DEALLOCATE(ptrrow)
      DEALLOCATE(ptrcol)
      DEALLOCATE(row_index_list)
      DEALLOCATE(col_index_list)
      RETURN
      CONTAINS
        SUBROUTINE smumps_set_lda_shift_val_son(IW, LIW, IOLDPS,
     &                             LDA, SHIFT_VAL_SON)
        INTEGER,    INTENT(IN)  :: LIW, IOLDPS
        INTEGER,    INTENT(IN)  :: IW(LIW)
        INTEGER,    INTENT(OUT) :: LDA
        INTEGER(8), INTENT(OUT) :: SHIFT_VAL_SON
        INCLUDE 'mumps_headers.h'
        INTEGER :: LCONT, NROW, NPIV, NASS, NELIM
        lcont  = iw(ioldps+keep(ixsz))
        nrow   = iw(ioldps+2+keep(ixsz))
        npiv   = iw(ioldps+3+keep(ixsz))
        nass   = iw(ioldps+4+keep(ixsz))
        nelim  = nass-npiv
        IF (iw(ioldps+xxs).EQ.s_nolcbnocontig38.OR.
     &      iw(ioldps+xxs).EQ.s_all) THEN
          shift_val_son      = int(npiv,8)
          lda                = lcont + npiv
        ELSE IF (iw(ioldps+xxs).EQ.s_nolcbcontig38) THEN
          shift_val_son = int(nrow,8)*int(lcont+npiv-nelim,8)
          lda           = nelim
        ELSE IF (iw(ioldps+xxs).EQ.s_nolcleaned38) THEN
          shift_val_son=0_8
          lda = nelim
       ELSE
          WRITE(*,*) myid,
     &    ": internal error in SMUMPS_SET_LDA_SHIFT_VAL_SON",
     &    iw(ioldps+xxs), "ISON=",ison
          CALL mumps_abort()
        ENDIF
        RETURN
        END SUBROUTINE smumps_set_lda_shift_val_son
      END SUBROUTINE smumps_build_and_send_cb_root
      SUBROUTINE smumps_root_local_assembly( N, VAL_ROOT,
     &   LOCAL_M, LOCAL_N,
     &   NPCOL, NPROW, MBLOCK, NBLOCK, NBCOL_SON, NBROW_SON, INDCOL_SON,
     &   INDROW_SON, LD_SON, VAL_SON, SUBSET_ROW, SUBSET_COL,
     &   NSUBSET_ROW, NSUBSET_COL, NSUPROW, NSUPCOL,
     &   RG2L_ROW, RG2L_COL, TRANSPOSE_ASM,
     &   KEEP, RHS_ROOT, NLOC  )
      USE smumps_struc_def, ONLY : smumps_root_struc
      IMPLICIT NONE
      INTEGER N, LOCAL_M, LOCAL_N
      REAL VAL_ROOT( LOCAL_M, LOCAL_N )
      INTEGER NPCOL, NPROW, MBLOCK, NBLOCK
      INTEGER NBCOL_SON, NBROW_SON
      INTEGER INDCOL_SON( NBCOL_SON ), INDROW_SON( NBROW_SON )
      INTEGER LD_SON 
      INTEGER NSUPROW, NSUPCOL
      REAL VAL_SON( LD_SON, NBROW_SON )
      INTEGER KEEP(500)
      INTEGER NSUBSET_ROW, NSUBSET_COL
      INTEGER SUBSET_ROW( NSUBSET_ROW ), SUBSET_COL( NSUBSET_COL )
      INTEGER RG2L_ROW( N ), RG2L_COL( N )
      LOGICAL TRANSPOSE_ASM
      INTEGER NLOC
      REAL RHS_ROOT( LOCAL_M, NLOC)
      INTEGER ISUB, JSUB, I, J, IPOS_ROOT, JPOS_ROOT
      INTEGER ILOC_ROOT, JLOC_ROOT, IGLOB, JGLOB
      IF (keep(50).EQ.0) THEN
        DO isub = 1, nsubset_row
          i         = subset_row( isub )
          iglob     = indrow_son( i )
          ipos_root = rg2l_row( iglob )   
          iloc_root = mblock
     &            * ( ( ipos_root - 1 ) / ( mblock * nprow ) )
     &            + mod( ipos_root - 1, mblock ) + 1
          DO jsub = 1, nsubset_col-nsupcol
            j         = subset_col( jsub )
            jglob     = indcol_son( j )
            jpos_root = rg2l_col( jglob )
            jloc_root = nblock
     &              * ( ( jpos_root - 1 ) / ( nblock * npcol ) )
     &              + mod( jpos_root - 1, nblock ) + 1
            val_root( iloc_root, jloc_root ) =
     &           val_root( iloc_root, jloc_root ) + val_son( j, i )
          END DO
          DO jsub = nsubset_col-nsupcol+1, nsubset_col
            j         = subset_col( jsub )
            jglob     = indcol_son( j )
             jpos_root = jglob - n  
             jloc_root = nblock
     &                * ( ( jpos_root - 1 ) / ( nblock * npcol ) )
     &                + mod( jpos_root - 1, nblock ) + 1
             rhs_root(iloc_root, jloc_root) =  
     &            rhs_root(iloc_root, jloc_root) + val_son( j, i )
          ENDDO
        END DO
      ELSE
        IF ( .NOT. transpose_asm ) THEN
          DO isub = 1, nsubset_row - nsuprow 
            i         = subset_row( isub )
            iglob     = indrow_son( i )
            ipos_root = rg2l_row( iglob )
            iloc_root = mblock
     &            * ( ( ipos_root - 1 ) / ( mblock * nprow ) )
     &            + mod( ipos_root - 1, mblock ) + 1
            DO jsub = 1, nsubset_col -nsupcol
              j         = subset_col( jsub )
              jglob     = indcol_son( j )
              jpos_root = rg2l_col( jglob )
              IF (keep(50).NE.0. and. jpos_root .GT. ipos_root) cycle
              jloc_root = nblock
     &                * ( ( jpos_root - 1 ) / ( nblock * npcol ) )
     &                + mod( jpos_root - 1, nblock ) + 1
              val_root( iloc_root, jloc_root ) =
     &            val_root( iloc_root, jloc_root ) + val_son( j, i )
            END DO
          END DO
          DO jsub = nsubset_col -nsupcol+1, nsubset_col
            j         = subset_col( jsub )
            jglob     = indrow_son( j )  
            jpos_root = jglob - n  
            jloc_root = nblock
     &                * ( ( jpos_root - 1 ) / ( nblock * npcol ) )
     &                + mod( jpos_root - 1, nblock ) + 1
            DO isub = nsubset_row - nsuprow +1, nsubset_row
              i         = subset_row( isub )
              iglob     = indcol_son( i )  
              ipos_root = rg2l_row(iglob)
              iloc_root = mblock
     &            * ( ( ipos_root - 1 ) / ( mblock * nprow ) )
     &            + mod( ipos_root - 1, mblock ) + 1
              rhs_root(iloc_root, jloc_root) =  
     &            rhs_root(iloc_root, jloc_root) + val_son( i, j )
            END DO
          END DO
        ELSE
          DO isub = 1, nsubset_col-nsupcol 
            i         = subset_col( isub )
            iglob     = indrow_son( i )
            jpos_root = rg2l_col( iglob )
            jloc_root = nblock
     &                * ( ( jpos_root - 1 ) / ( nblock * npcol ) )
     &                + mod( jpos_root - 1, nblock ) + 1
            DO jsub = 1, nsubset_row
              j         = subset_row( jsub )
              jglob     = indcol_son( j )
              ipos_root = rg2l_row( jglob )  
              iloc_root = mblock
     &                * ( ( ipos_root - 1 ) / ( mblock * nprow ) )
     &                + mod( ipos_root - 1, mblock ) + 1
              val_root( iloc_root, jloc_root ) =
     &            val_root( iloc_root, jloc_root ) + val_son( j, i )
            END DO
           ENDDO
           DO isub = nsubset_col-nsupcol+1, nsubset_col
            i         = subset_col( isub )
            iglob     = indrow_son( i )
            jpos_root = iglob - n 
            jloc_root = nblock
     &                * ( ( jpos_root - 1 ) / ( nblock * npcol ) )
     &                + mod( jpos_root - 1, nblock ) + 1
            DO jsub = 1, nsubset_row
              j         = subset_row( jsub )
              jglob     = indcol_son( j )
              ipos_root = rg2l_row( jglob )  
              iloc_root = mblock
     &                * ( ( ipos_root - 1 ) / ( mblock * nprow ) )
     &                + mod( ipos_root - 1, mblock ) + 1
              rhs_root( iloc_root, jloc_root ) =
     &            rhs_root( iloc_root, jloc_root ) + val_son( j, i )
            END DO
           ENDDO
        END IF
      END IF
      RETURN
      END SUBROUTINE smumps_root_local_assembly
      SUBROUTINE smumps_init_root_ana
     &( myid, nprocs, n, root, comm_root, iroot, fils,
     &  k50, k46, k51
     &     , k60, idnprow, idnpcol, idmblock, idnblock
     & )
      USE smumps_struc_def, ONLY : smumps_root_struc
      IMPLICIT NONE
      INTEGER MYID, MYID_ROOT
      TYPE (SMUMPS_ROOT_STRUC)::root
      INTEGER COMM_ROOT
      INTEGER N, IROOT, NPROCS, K50, K46, K51
      INTEGER FILS( N )
      INTEGER K60, IDNPROW, IDNPCOL, IDMBLOCK, IDNBLOCK
      INTEGER INODE, NPROWtemp, NPCOLtemp
      LOGICAL SLAVE
      root%ROOT_SIZE     = 0
      root%TOT_ROOT_SIZE = 0
      slave = ( myid .ne. 0 .or.
     &        ( myid .eq. 0 .and. k46 .eq. 1 ) )
      inode = iroot
      DO WHILE ( inode .GT. 0 )
        inode = fils( inode )
        root%ROOT_SIZE = root%ROOT_SIZE + 1
      END DO
      IF ( ( k60 .NE. 2 .AND. k60 .NE. 3 ) .OR.
     &       idnprow .LE. 0 .OR. idnpcol .LE. 0
     &      .OR. idmblock .LE.0 .OR. idnblock.LE.0
     &      .OR. idnprow * idnpcol .GT. nprocs ) THEN
        root%MBLOCK = k51
        root%NBLOCK = k51
        CALL smumps_def_grid( nprocs, root%NPROW, root%NPCOL,
     &                         root%ROOT_SIZE, k50 )
        IF  ( k60 .EQ. 2 .OR. k60 .EQ. 3 ) THEN
          idnprow = root%NPROW
          idnpcol = root%NPCOL
          idmblock = root%MBLOCK
          idnblock = root%NBLOCK
        ENDIF
      ELSE IF  ( k60 .EQ. 2 .OR. k60 .EQ. 3 ) THEN
        root%NPROW = idnprow
        root%NPCOL = idnpcol
        root%MBLOCK = idmblock
        root%NBLOCK = idnblock
      ENDIF
      IF  ( k60 .EQ. 2 .OR. k60 .EQ. 3 ) THEN
        IF (slave) THEN
          root%LPIV = 0
          IF (k46.EQ.0) THEN
            myid_root=myid-1
          ELSE
            myid_root=myid
          ENDIF
          IF (myid_root < root%NPROW*root%NPCOL) THEN
            root%MYROW = myid_root / root%NPCOL
            root%MYCOL = mod(myid_root, root%NPCOL)
            root%yes  = .true.
          ELSE
            root%MYROW = -1
            root%MYCOL = -1
            root%yes  = .false.
          ENDIF
        ELSE
          root%yes  = .false.
        ENDIF
      ELSE IF ( slave ) THEN
        IF ( root%gridinit_done) THEN
           IF (root%yes) THEN
             CALL blacs_gridexit( root%CNTXT_BLACS )
             root%gridinit_done = .false.
           ENDIF
        END IF
        root%CNTXT_BLACS = comm_root
        CALL blacs_gridinit( root%CNTXT_BLACS, 'R',
     &                       root%NPROW, root%NPCOL )
        root%gridinit_done = .true.
        CALL blacs_gridinfo( root%CNTXT_BLACS,
     &                       nprowtemp, npcoltemp,
     &                       root%MYROW, root%MYCOL )
        IF ( root%MYROW .NE. -1 ) THEN
          root%yes = .true.
        ELSE
          root%yes = .false.
        END IF
        root%LPIV = 0
      ELSE
        root%yes = .false.
      ENDIF
      RETURN
      END SUBROUTINE smumps_init_root_ana
      SUBROUTINE smumps_init_root_fac( N, root, FILS, IROOT,
     &                                 KEEP, INFO )
      USE smumps_struc_def, ONLY : smumps_root_struc
      IMPLICIT NONE
      TYPE ( SMUMPS_ROOT_STRUC ):: root
      INTEGER N, IROOT, INFO(80), KEEP(500)
      INTEGER FILS( N )
      INTEGER INODE, I, allocok
      IF ( associated( root%RG2L_ROW ) ) THEN
        DEALLOCATE( root%RG2L_ROW )
        NULLIFY( root%RG2L_ROW )
      ENDIF
      IF ( associated( root%RG2L_COL ) ) THEN
        DEALLOCATE( root%RG2L_COL )
        NULLIFY( root%RG2L_COL )
      ENDIF
      ALLOCATE( root%RG2L_ROW( n ), stat = allocok )
      IF ( allocok .GT. 0 ) THEN
        info(1)=-13
        info(2)=n
        RETURN
      ENDIF
      ALLOCATE( root%RG2L_COL( n ), stat = allocok )
      IF ( allocok .GT. 0 ) THEN
        DEALLOCATE( root%RG2L_ROW ); NULLIFY( root%RG2L_ROW )
        info(1)=-13
        info(2)=n
        RETURN
      ENDIF
      inode = iroot
      i = 1
      DO WHILE ( inode .GT. 0 )
        root%RG2L_ROW( inode ) = i
        root%RG2L_COL( inode ) = i
        i = i + 1
        inode = fils( inode )
      END DO
      root%TOT_ROOT_SIZE=0
      RETURN
      END SUBROUTINE smumps_init_root_fac
      SUBROUTINE smumps_def_grid( NPROCS, NPROW, NPCOL, SIZE, K50 )
      IMPLICIT NONE
      INTEGER NPROCS, NPROW, NPCOL, SIZE, K50
      INTEGER NPROWtemp, NPCOLtemp, NPROCSused, FLATNESS
      LOGICAL KEEPIT
      IF ( k50 .EQ. 1 ) THEN
        flatness = 2
      ELSE
        flatness = 3
      ENDIF
      nprow  = int(sqrt(real(nprocs)))
      nprowtemp = nprow
      npcol  = int(nprocs / nprow)
      npcoltemp = npcol
      nprocsused = nprowtemp * npcoltemp
 10   CONTINUE
      IF ( nprowtemp >= npcoltemp/flatness .AND. nprowtemp > 1) THEN
        nprowtemp = nprowtemp - 1
        npcoltemp = int(nprocs / nprowtemp)
        keepit=.false.
        IF ( nprowtemp * npcoltemp .GE. nprocsused ) THEN
          IF ( ( k50 .NE. 1 .AND. nprowtemp >= npcoltemp/flatness)
     &         .OR. nprowtemp * npcoltemp .GT. nprocsused )
     &         keepit=.true.
        END IF
        IF ( keepit ) THEN
          nprow = nprowtemp
          npcol = npcoltemp
          nprocsused = nprow * npcol
        END IF
        GO TO 10
      END IF
      RETURN
      END SUBROUTINE smumps_def_grid
      SUBROUTINE smumps_scatter_root(MYID, M, N, ASEQ,
     &                    LOCAL_M, LOCAL_N,
     &                    MBLOCK, NBLOCK,
     &                    APAR,
     &                    MASTER_ROOT,
     &                    NPROW, NPCOL,
     &                    COMM)
      IMPLICIT NONE
      INTEGER MYID, MASTER_ROOT, COMM
      INTEGER M, N
      INTEGER NPROW, NPCOL
      INTEGER LOCAL_M, LOCAL_N
      INTEGER MBLOCK, NBLOCK
      REAL APAR( LOCAL_M, LOCAL_N )
      REAL ASEQ( M, N )
      include 'mpif.h'
      INTEGER I, J, SIZE_IBLOCK, SIZE_JBLOCK, IDEST, IROW, ICOL
      INTEGER IBLOCK, JBLOCK, II, JJ, KK
      INTEGER IAPAR, JAPAR, IERR, allocok
      INTEGER :: STATUS(MPI_STATUS_SIZE)
      REAL, DIMENSION(:), ALLOCATABLE :: WK
      LOGICAL JUPDATE
      ALLOCATE(WK( MBLOCK * NBLOCK ), stat=allocok)
      if(allocok.ne.0) then
         WRITE(6,*) ' Allocation error of WK in '
     &       // 'routine SMUMPS_SCATTER_ROOT '
         CALL mumps_abort()
      endif
        iapar = 1
        japar = 1
        DO j = 1, n, nblock
          size_jblock = nblock
          IF ( j + nblock > n ) THEN
            size_jblock = n - j + 1
          END IF
          jupdate = .false.
          DO i = 1, m, mblock
            size_iblock = mblock
            IF ( i + mblock > m ) THEN
              size_iblock = m - i + 1
            END IF
            iblock = i / mblock
            jblock = j / nblock
            irow = mod( iblock, nprow )
            icol = mod( jblock, npcol )
            idest = irow * npcol + icol
            IF ( idest .NE. master_root ) THEN
              IF ( myid .EQ. master_root ) THEN
                kk=1
                DO jj=j,j+size_jblock-1
                DO ii=i,i+size_iblock-1
                  wk(kk)=aseq(ii,jj)
                  kk=kk+1
                END DO
                END DO
                CALL mpi_ssend( wk, size_iblock*size_jblock,
     &                         mpi_real,
     &                         idest, 128, comm, ierr )
              ELSE IF ( myid .EQ. idest ) THEN
                CALL mpi_recv( wk(1),
     &                         size_iblock*size_jblock,
     &                         mpi_real,
     &                         master_root,128,comm,status,ierr)
                kk=1
                DO jj=japar,japar+size_jblock-1
                DO ii=iapar,iapar+size_iblock-1
                  apar(ii,jj)=wk(kk)
                  kk=kk+1
                END DO
                END DO
                jupdate = .true.
                iapar = iapar + size_iblock
              END IF
            ELSE IF ( myid.EQ. master_root ) THEN
              apar( iapar:iapar+size_iblock-1,
     &              japar:japar+size_jblock-1 )
     &        = aseq(i:i+size_iblock-1,j:j+size_jblock-1)
              jupdate = .true.
              iapar = iapar + size_iblock
            END IF
          END DO
          IF ( jupdate ) THEN
            iapar = 1
            japar = japar + size_jblock
          END IF
        END DO
        DEALLOCATE(wk)
      RETURN
      END SUBROUTINE smumps_scatter_root
      SUBROUTINE smumps_gather_root(MYID, M, N, ASEQ,
     &                    LOCAL_M, LOCAL_N,
     &                    MBLOCK, NBLOCK,
     &                    APAR,
     &                    MASTER_ROOT,
     &                    NPROW, NPCOL,
     &                    COMM)
      IMPLICIT NONE
      INTEGER MYID, MASTER_ROOT, COMM
      INTEGER M, N
      INTEGER NPROW, NPCOL
      INTEGER LOCAL_M, LOCAL_N
      INTEGER MBLOCK, NBLOCK
      REAL APAR( LOCAL_M, LOCAL_N )
      REAL ASEQ( M, N )
      include 'mpif.h'
      INTEGER I, J, SIZE_IBLOCK, SIZE_JBLOCK, ISOUR, IROW, ICOL
      INTEGER IBLOCK, JBLOCK, II, JJ, KK
      INTEGER IAPAR, JAPAR, IERR, allocok
      INTEGER :: STATUS(MPI_STATUS_SIZE)
      REAL,DIMENSION(:), ALLOCATABLE :: WK
      LOGICAL JUPDATE
      ALLOCATE(wk( mblock * nblock ), stat=allocok)
      if(allocok.ne.0) then
         WRITE(6,*) ' Allocation error of WK in '
     &       // 'routine SMUMPS_GATHER_ROOT '
         CALL mumps_abort()
      endif
        iapar = 1
        japar = 1
        DO j = 1, n, nblock
          size_jblock = nblock
          IF ( j + nblock > n ) THEN
            size_jblock = n - j + 1
          END IF
          jupdate = .false.
          DO i = 1, m, mblock
            size_iblock = mblock
            IF ( i + mblock > m ) THEN
              size_iblock = m - i + 1
            END IF
            iblock = i / mblock
            jblock = j / nblock
            irow = mod( iblock, nprow )
            icol = mod( jblock, npcol )
            isour = irow * npcol + icol
            IF ( isour .NE. master_root ) THEN
              IF ( myid .EQ. master_root ) THEN
                CALL mpi_recv( wk(1), size_iblock*size_jblock,
     &                         mpi_real,
     &                         isour, 128, comm, status, ierr )
                kk=1
                DO jj=j,j+size_jblock-1
                DO ii=i,i+size_iblock-1
                  aseq(ii,jj)=wk(kk)
                  kk=kk+1
                END DO
                END DO
              ELSE IF ( myid .EQ. isour ) THEN
                kk=1
                DO jj=japar,japar+size_jblock-1
                DO ii=iapar,iapar+size_iblock-1
                  wk(kk)=apar(ii,jj)
                  kk=kk+1
                END DO
                END DO
                CALL mpi_ssend( wk( 1 ),
     &                         size_iblock*size_jblock,
     &                         mpi_real,
     &                         master_root,128,comm,ierr)
                jupdate = .true.
                iapar = iapar + size_iblock
              END IF
            ELSE IF ( myid.EQ. master_root ) THEN
              aseq(i:i+size_iblock-1,j:j+size_jblock-1)
     &        = apar( iapar:iapar+size_iblock-1,
     &                japar:japar+size_jblock-1 )
              jupdate = .true.
              iapar = iapar + size_iblock
            END IF
          END DO
          IF ( jupdate ) THEN
            iapar = 1
            japar = japar + size_jblock
          END IF
        END DO
        DEALLOCATE(wk)
      RETURN
      END SUBROUTINE smumps_gather_root
      SUBROUTINE smumps_root_alloc_static(root, IROOT, N,
     &                  IW, LIW, A, LA,
     &                  FILS, DAD, MYID, SLAVEF, PROCNODE_STEPS,
     &                  LPTRAR, NELT, FRTPTR, FRTELT,
     &                  PTRAIW, PTRARW, INTARR, DBLARR,
     &                  LRLU, IPTRLU,
     &                  IWPOS, IWPOSCB, PTRIST, PTRAST,
     &                  STEP, PIMASTER, PAMASTER, ITLOC, RHS_MUMPS,
     &                  COMP, LRLUS, IFLAG, KEEP,KEEP8,DKEEP,IERROR )
      USE smumps_struc_def, ONLY : smumps_root_struc
      IMPLICIT NONE
      INTEGER MYID
      INTEGER KEEP(500)
      INTEGER(8) KEEP8(150)
      REAL DKEEP(230)
      TYPE (SMUMPS_ROOT_STRUC ) :: root
      INTEGER(8) :: LA, LRLU, IPTRLU, LRLUS
      INTEGER IROOT, LIW, N, IWPOS, IWPOSCB
      INTEGER IW( LIW )
      REAL A( LA )
      INTEGER, INTENT(IN) :: SLAVEF
      INTEGER, INTENT(IN) :: PROCNODE_STEPS(KEEP(28))
      INTEGER PTRIST(KEEP(28)), STEP(N)
      INTEGER(8) :: PTRAST(KEEP(28)), PAMASTER(KEEP(28))
      INTEGER PIMASTER(KEEP(28))
      INTEGER ITLOC( N + KEEP(253) )
      REAL :: RHS_MUMPS(KEEP(255))
      INTEGER COMP, IFLAG, IERROR
      include 'mumps_headers.h'
      INTEGER FILS( N ), DAD(KEEP(28))
      INTEGER LPTRAR, NELT
      INTEGER FRTPTR( N+1), FRTELT( NELT )
      INTEGER(8), INTENT(IN) :: PTRAIW( LPTRAR ), PTRARW( LPTRAR )
      INTEGER INTARR(KEEP8(27))
      REAL DBLARR(KEEP8(26))
      INTEGER numroc
      EXTERNAL numroc
      REAL ZERO
      PARAMETER( ZERO = 0.0e0 )
      INTEGER(8) :: LREQA_ROOT
      INTEGER LREQI_ROOT, LOCAL_M, LOCAL_N, allocok
      LOGICAL :: EARLYT3ROOTINS
          local_m = numroc( root%ROOT_SIZE, root%MBLOCK,
     &              root%MYROW, 0, root%NPROW )
          local_m = max( 1, local_m )
          local_n = numroc( root%ROOT_SIZE, root%NBLOCK,
     &              root%MYCOL, 0, root%NPCOL )
          IF (keep(253).GT.0) THEN
            root%RHS_NLOC = numroc( keep(253), root%NBLOCK,
     &              root%MYCOL, 0, root%NPCOL )
            root%RHS_NLOC = max(1, root%RHS_NLOC)
          ELSE
            root%RHS_NLOC = 1
          ENDIF
          IF (associated( root%RHS_ROOT) ) 
     &             DEALLOCATE (root%RHS_ROOT)
          ALLOCATE(root%RHS_ROOT(local_m,root%RHS_NLOC), 
     &              stat=allocok)
          IF ( allocok.GT.0) THEN
            iflag=-13
            ierror = local_m*root%RHS_NLOC
            RETURN
          ENDIF
          IF (keep(253).NE.0) THEN
            root%RHS_ROOT = zero
            CALL smumps_asm_rhs_root (  n, fils, 
     &           root, keep, rhs_mumps, 
     &           iflag, ierror )
            IF ( iflag .LT. 0 ) RETURN
          ENDIF
          IF (keep(60) .NE. 0) THEN
            ptrist(step(iroot)) = -6666666
          ELSE
            lreqi_root = 2 + keep(ixsz)
            lreqa_root = int(local_m,8) * int(local_n,8)
            IF (lreqa_root.EQ.0_8) THEN
              ptrist(step(iroot)) = -9999999
              RETURN
            ENDIF
            CALL smumps_alloc_cb(.false.,0_8,.false.,.false.,
     &                     myid,n,keep,keep8,dkeep,iw,liw,a,la,
     &                     lrlu, iptrlu,
     &                     iwpos, iwposcb, slavef, procnode_steps, dad,
     &                     ptrist, ptrast,
     &                     step, pimaster, pamaster, lreqi_root,
     &                     lreqa_root, iroot, s_notfree, .true., comp,
     &                     lrlus, keep8(67), iflag, ierror
     &           )
            IF ( iflag .LT. 0 ) RETURN
            ptrist( step(iroot) ) = iwposcb + 1
            pamaster( step(iroot) ) = iptrlu  + 1_8
            iw( iwposcb + 1 + keep(ixsz)) = - local_n
            iw( iwposcb + 2 + keep(ixsz)) =   local_m
          ENDIF
          earlyt3rootins = keep(200) .EQ.0
     &        .OR.  ( keep(200) .LT. 0 .AND. keep(400) .EQ. 0 )
          IF (local_n > 0 .AND. .NOT. earlyt3rootins ) THEN
            IF (keep(60) .EQ. 0) THEN
              CALL smumps_set_to_zero(a(iptrlu+1_8), local_m,
     &        local_m, local_n, keep)
            ELSE
              CALL smumps_set_to_zero(root%SCHUR_POINTER(1),
     &        root%SCHUR_LLD, local_m, local_n, keep)
            ENDIF
            IF (keep(55) .eq. 0) THEN
              IF (keep(60) .EQ. 0) THEN
                CALL smumps_asm_arr_root( n, root, iroot,
     &          a(iptrlu+1_8), local_m, local_m, local_n,
     &          fils, ptraiw, ptrarw, intarr, dblarr,
     &          keep8(27), keep8(26), myid )
              ELSE
                CALL smumps_asm_arr_root( n, root, iroot,
     &          root%SCHUR_POINTER(1), root%SCHUR_LLD, local_m, local_n,
     &          fils, ptraiw, ptrarw, intarr, dblarr,
     &          keep8(27), keep8(26), myid )
              ENDIF
            ELSE
              IF (keep(60) .EQ. 0) THEN
                CALL smumps_asm_elt_root( n, root,
     &          a(iptrlu+1_8), local_m, local_m, local_n,
     &          lptrar, nelt, frtptr, frtelt,
     &          ptraiw, ptrarw, intarr, dblarr,
     &          keep8(27), keep8(26), keep, keep8, myid )
              ELSE
                CALL smumps_asm_elt_root( n, root,
     &           root%SCHUR_POINTER(1), root%SCHUR_LLD,
     &           root%SCHUR_MLOC, root%SCHUR_NLOC,
     &           lptrar, nelt, frtptr, frtelt,
     &           ptraiw, ptrarw, intarr, dblarr,
     &           keep8(27), keep8(26), keep, keep8, myid )
              ENDIF
            ENDIF
          ENDIF
      RETURN
      END SUBROUTINE smumps_root_alloc_static
      SUBROUTINE smumps_asm_elt_root( N, root, 
     &       VALROOT, LOCAL_M_LLD, LOCAL_M, LOCAL_N,
     &       LPTRAR, NELT, FRTPTR, FRTELT,
     &       PTRAIW, PTRARW,
     &       INTARR, DBLARR, LINTARR, LDBLARR,
     &       KEEP, KEEP8,
     &       MYID)
      USE smumps_struc_def, ONLY : smumps_root_struc
      IMPLICIT NONE
      TYPE (SMUMPS_ROOT_STRUC) :: root
      INTEGER :: N, MYID, LOCAL_M, LOCAL_N, KEEP(500)
      INTEGER :: LOCAL_M_LLD
      INTEGER(8) KEEP8(150)
      REAL VALROOT(LOCAL_M_LLD,LOCAL_N)
      INTEGER LPTRAR, NELT
      INTEGER FRTPTR( N+1), FRTELT( NELT )
      INTEGER(8), INTENT(IN) :: PTRARW( LPTRAR ), PTRAIW( LPTRAR )
      INTEGER(8), INTENT(IN) :: LINTARR, LDBLARR
      INTEGER, INTENT(INOUT) :: INTARR(LINTARR)
      REAL DBLARR(LDBLARR)
      INTEGER(8) :: J1, J2, K8, IPTR
      INTEGER    :: IELT, I, J, IGLOB, JGLOB, SIZEI, IBEG
      INTEGER    :: ARROW_ROOT
      INTEGER    :: IPOSROOT, JPOSROOT, IROW_GRID, JCOL_GRID
      INTEGER    :: ILOCROOT, JLOCROOT
      ARROW_ROOT = 0
      DO iptr = frtptr(keep(38)), frtptr(keep(38)+1) - 1
        ielt = frtelt( iptr )
        j1 = ptraiw(ielt)
        j2 = ptraiw(ielt+1)-1
        k8 = ptrarw(ielt)
        sizei=int(j2-j1)+1
        DO j=1, sizei
          jglob          = intarr(j1+j-1)
          intarr(j1+j-1) = root%RG2L_ROW(jglob)
        ENDDO
        DO j = 1, sizei
          jglob         = intarr(j1+j-1)
          IF ( keep(50).eq. 0 ) THEN
            ibeg = 1
          ELSE
            ibeg = j
          END IF
          DO i = ibeg, sizei
            iglob = intarr(j1+i-1)
            IF ( keep(50).eq.0 ) THEN
              iposroot = intarr(j1+i-1)
              jposroot = intarr(j1+j-1)
            ELSE
              IF ( intarr(j1+i-1).GT. intarr(j1+j-1) ) THEN
                iposroot = intarr(j1+i-1)
                jposroot = intarr(j1+j-1)
              ELSE
                iposroot = intarr(j1+j-1)
                jposroot = intarr(j1+i-1)
              END IF
            END IF
            irow_grid = mod( ( iposroot - 1 )/root%MBLOCK,
     &                         root%NPROW )
            jcol_grid = mod( ( jposroot - 1 )/root%NBLOCK,
     &                         root%NPCOL )
            IF ( irow_grid.EQ.root%MYROW .AND.
     &           jcol_grid.EQ.root%MYCOL ) THEN
              ilocroot = root%MBLOCK * ( ( iposroot - 1 ) /
     &               ( root%MBLOCK * root%NPROW ) )
     &             + mod( iposroot - 1, root%MBLOCK ) + 1
              jlocroot = root%NBLOCK * ( ( jposroot - 1 ) /
     &               ( root%NBLOCK * root%NPCOL ) )
     &             + mod( jposroot - 1, root%NBLOCK ) + 1
              valroot( ilocroot, jlocroot ) = 
     &        valroot( ilocroot, jlocroot ) + dblarr(k8)
            ENDIF
            k8 = k8 + 1_8
          END DO
        END DO
        arrow_root = arrow_root + int(ptrarw(ielt+1_8)-ptrarw(ielt))
      END DO
      keep(49) = arrow_root
      RETURN
      END SUBROUTINE smumps_asm_elt_root
      SUBROUTINE smumps_asm_rhs_root
     &           ( n, fils, root, keep, rhs_mumps, 
     &             iflag, ierror )
      USE smumps_struc_def, ONLY : smumps_root_struc
      IMPLICIT NONE
      INTEGER N, KEEP(500), IFLAG, IERROR
      INTEGER FILS(N)
      TYPE (SMUMPS_ROOT_STRUC ) :: root
      REAL :: RHS_MUMPS(KEEP(255))
      INTEGER JCOL, IPOS_ROOT, JPOS_ROOT,
     &        irow_grid, jcol_grid, ilocrhs, jlocrhs,
     &        inode
      inode = keep(38)
      DO WHILE (inode.GT.0)
        ipos_root = root%RG2L_ROW( inode ) 
        irow_grid  = mod( ( ipos_root - 1 ) / root%MBLOCK, root%NPROW )
        IF (  irow_grid .NE. root%MYROW ) GOTO 100 
        ilocrhs = root%MBLOCK * ( ( ipos_root - 1 ) /
     &                 ( root%MBLOCK * root%NPROW ) )
     &               + mod( ipos_root - 1, root%MBLOCK ) + 1
        DO jcol = 1, keep(253) 
          jpos_root = jcol
          jcol_grid  = mod((jpos_root-1)/root%NBLOCK, root%NPCOL)
          IF (jcol_grid.NE.root%MYCOL ) cycle
           jlocrhs = root%NBLOCK * ( ( jpos_root - 1 ) /
     &                 ( root%NBLOCK * root%NPCOL ) )
     &               + mod( jpos_root - 1, root%NBLOCK ) + 1
          root%RHS_ROOT(ilocrhs, jlocrhs) =
     &                 rhs_mumps(inode+(jcol-1)*keep(254))
        ENDDO
 100    CONTINUE
        inode=fils(inode)
      ENDDO
      RETURN
      END SUBROUTINE smumps_asm_rhs_root
      SUBROUTINE smumps_asm_arr_root( N, root, IROOT,
     &   VALROOT, LOCAL_M_LLD, LOCAL_M, LOCAL_N, FILS,
     &       PTRAIW, PTRARW,
     &       INTARR, DBLARR, LINTARR, LDBLARR,
     &       MYID)
      USE smumps_struc_def, ONLY : smumps_root_struc
      IMPLICIT NONE
      TYPE (SMUMPS_ROOT_STRUC) :: root
      INTEGER :: N, MYID, IROOT, LOCAL_M, LOCAL_N
      INTEGER :: LOCAL_M_LLD
      INTEGER FILS( N )
      INTEGER(8), INTENT(IN) :: PTRARW( N ), PTRAIW( N )
      REAL VALROOT(LOCAL_M_LLD,LOCAL_N)
      INTEGER(8), INTENT(IN) :: LINTARR, LDBLARR
      INTEGER INTARR(LINTARR)
      REAL DBLARR(LDBLARR)
      REAL VAL
      INTEGER(8) :: JJ, J1,JK, J2,J3, J4, AINPUT
      INTEGER IORG, IBROT, NUMORG,
     &        irow, jcol
      INTEGER IPOSROOT, JPOSROOT, IROW_GRID, JCOL_GRID
      INTEGER ILOCROOT, JLOCROOT
      numorg = root%ROOT_SIZE
      ibrot  = iroot
      DO iorg = 1, numorg
        jk = ptraiw(ibrot)
        ainput = ptrarw(ibrot)
        ibrot = fils(ibrot)
        jj = jk + 1
        j1 = jj + 1
        j2 = j1 + intarr(jk)
        j3 = j2 + 1
        j4 = j2 - intarr(jj)
        jcol = intarr(j1)
        DO jj = j1, j2
         irow = intarr(jj)
         val  = dblarr(ainput)
         ainput = ainput + 1
         iposroot = root%RG2L_ROW( irow )
         jposroot = root%RG2L_COL( jcol )
         irow_grid  = mod( ( iposroot - 1 ) / root%MBLOCK, root%NPROW )
         jcol_grid  = mod( ( jposroot - 1 ) / root%NBLOCK, root%NPCOL )
         IF ( irow_grid .EQ. root%MYROW .AND.
     &        jcol_grid .EQ. root%MYCOL ) THEN
            ilocroot = root%MBLOCK * ( ( iposroot - 1 ) /
     &                 ( root%MBLOCK * root%NPROW ) )
     &               + mod( iposroot - 1, root%MBLOCK ) + 1
            jlocroot = root%NBLOCK * ( ( jposroot - 1 ) /
     &                 ( root%NBLOCK * root%NPCOL ) )
     &               + mod( jposroot - 1, root%NBLOCK ) + 1
            valroot( ilocroot, jlocroot ) =
     &      valroot( ilocroot, jlocroot ) + val
         END IF
        END DO
        IF (j3 .LE. j4) THEN
         irow =  intarr(j1)
         DO jj= j3,j4
          jcol = intarr(jj)
          val  = dblarr(ainput)
          ainput = ainput + 1
          iposroot = root%RG2L_ROW( irow )
          jposroot = root%RG2L_COL( jcol )
          irow_grid= mod( ( iposroot - 1 )/root%MBLOCK, root%NPROW)
          jcol_grid= mod( ( jposroot - 1 )/root%NBLOCK, root%NPCOL)
          IF ( irow_grid .EQ. root%MYROW .AND.
     &        jcol_grid .EQ. root%MYCOL ) THEN
            ilocroot = root%MBLOCK * ( ( iposroot - 1 ) /
     &                 ( root%MBLOCK * root%NPROW ) )
     &               + mod( iposroot - 1, root%MBLOCK ) + 1
            jlocroot = root%NBLOCK * ( ( jposroot - 1 ) /
     &                 ( root%NBLOCK * root%NPCOL ) )
     &               + mod( jposroot - 1, root%NBLOCK ) + 1
            valroot( ilocroot, jlocroot ) = 
     &      valroot( ilocroot, jlocroot ) + val
          END IF
         END DO
        ENDIF
      ENDDO
      RETURN
      END SUBROUTINE smumps_asm_arr_root