cfac_asm_ELT.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 cmumps_elt_asm_s_2_s_init(
     &    NELT, FRT_PTR, FRT_ELT,
     &    N, INODE, IW, LIW, A, LA, 
     &    NBROWS, NBCOLS,
     &    OPASSW, OPELIW, STEP, PTRIST, PTRAST, ITLOC,
     &    RHS_MUMPS,
     &    FILS, PTRARW, PTRAIW, INTARR, DBLARR, 
     &    ICNTL, KEEP, KEEP8, MYID, LRGROUPS)
      USE cmumps_dynamic_memory_m, ONLY : cmumps_dm_set_dynptr
      IMPLICIT NONE
      INTEGER NELT, N,LIW
      INTEGER(8) :: LA
      INTEGER KEEP(500), ICNTL(60)
      INTEGER(8) KEEP8(150)
      INTEGER INODE, MYID
      INTEGER NBROWS, NBCOLS 
      INTEGER(8) :: PTRAST(KEEP(28))
      INTEGER IW(LIW), ITLOC(N+KEEP(253)), STEP(N),
     &        ptrist(keep(28)), fils(n)
      INTEGER(8), INTENT(IN) :: PTRARW(NELT+1), PTRAIW(NELT+1)
      COMPLEX :: RHS_MUMPS(KEEP(255))
      INTEGER INTARR(KEEP8(27))
      INTEGER FRT_PTR(N+1), FRT_ELT(NELT)
      COMPLEX :: A(LA)
      COMPLEX :: DBLARR(KEEP8(26))
      DOUBLE PRECISION OPASSW, OPELIW
      INTEGER, INTENT(IN)    :: LRGROUPS(N)
      INTEGER(8) :: POSELT
      COMPLEX, DIMENSION(:), POINTER :: A_PTR
      INTEGER(8) :: LA_PTR
      INTEGER IOLDPS, NBCOLF, NBROWF, NSLAVES, HF,
     &        k1,k2,k,j,jpos,nass
      COMPLEX ZERO
      parameter( zero = (0.0e0,0.0e0) )
      include 'mumps_headers.h'
      ioldps  = ptrist(step(inode))
      CALL cmumps_dm_set_dynptr( iw(ioldps+xxs), a, la,
     &     ptrast(step(inode)), iw(ioldps+xxd), iw(ioldps+xxr),
     &     a_ptr, poselt, la_ptr )
      nbcolf  = iw(ioldps+keep(ixsz))
      nbrowf  = iw(ioldps+2+keep(ixsz))
      nass    = iw(ioldps+1+keep(ixsz))
      nslaves = iw(ioldps+5+keep(ixsz))
      hf      = 6 + nslaves + keep(ixsz)
      IF (nass.LT.0) THEN
          nass         = -nass
          iw(ioldps+1+keep(ixsz)) = nass
          CALL cmumps_asm_slave_elements( inode, n, nelt, iw, liw,
     &    ioldps, a_ptr(poselt), la_ptr, 1_8, keep, keep8, itloc, fils,
     &    ptraiw, ptrarw,
     &    intarr, dblarr, keep8(27), keep8(26), frt_ptr, frt_elt,
     &    rhs_mumps, lrgroups)
      ENDIF
      IF (nbrows.GT.0) THEN
          k1 = ioldps + hf + nbrowf
          k2 = k1 + nbcolf - 1
          jpos = 1
          DO k = k1, k2
           j        = iw(k)
           itloc(j) = jpos
           jpos     = jpos + 1
          ENDDO
      ENDIF
      RETURN
      END SUBROUTINE cmumps_elt_asm_s_2_s_init
      SUBROUTINE cmumps_asm_slave_elements( INODE, N, NELT, IW, LIW,
     &IOLDPS, A, LA, POSELT, KEEP, KEEP8, ITLOC, FILS, PTRAIW, PTRARW,
     &INTARR, DBLARR, LINTARR, LDBLARR, FRT_PTR, FRT_ELT, RHS_MUMPS,
     &LRGROUPS)
!$    USE OMP_LIB
      USE cmumps_ana_lr,    ONLY : get_cut
      USE cmumps_lr_core,   ONLY : max_cluster
      USE mumps_lr_common,  ONLY : compute_blr_vcs
      IMPLICIT NONE
      INTEGER, intent(in)    :: N, NELT, LIW, IOLDPS, INODE
      INTEGER(8), intent(in) :: LA, POSELT, LINTARR, LDBLARR
      INTEGER, intent(in)    :: IW(LIW)
      INTEGER, intent(in)    :: KEEP(500)
      INTEGER(8), intent(in) :: KEEP8(150)
      INTEGER, intent(inout) :: ITLOC(N+KEEP(253))
      COMPLEX, intent(inout) :: A(LA)
      COMPLEX, intent(in)    :: RHS_MUMPS(KEEP(255))
      INTEGER, intent(in)    :: INTARR(LINTARR)
      COMPLEX, intent(in)    :: DBLARR(LDBLARR)
      INTEGER, intent(in)    :: FRT_PTR(N+1), FRT_ELT(NELT)
      INTEGER, intent(in)    :: FILS(N)
      INTEGER(8), intent(in) :: PTRAIW(NELT+1), PTRARW(NELT+1)
      INTEGER, INTENT(IN)    :: LRGROUPS(N)
!$    INTEGER :: NOMP
!$    INTEGER(8) :: CHUNK8  
!$    INTEGER    :: CHUNK
      include 'mumps_headers.h'
      INTEGER    :: HF, NBROWF, NBCOLF, NASS, NSLAVES
      INTEGER    :: ILOC, IELL, ELTI, ELBEG, NUMELT
      INTEGER(8) :: SIZE_ELTI8 
      INTEGER    :: I, J, K, K1, K2
      INTEGER    :: IPOS, IPOS1, IPOS2, JPOS, IJROW
      INTEGER    :: IN
      INTEGER(8) :: II8, JJ8, J18, J28 
      INTEGER(8) :: AINPUT8
      INTEGER(8) :: AII8
      INTEGER(8) :: APOS, APOS2, ICT12
      INTEGER, POINTER, DIMENSION(:) :: BEGS_BLR_LS
      INTEGER :: NB_BLR_LS, NPARTSCB, NPARTSASS, MAXI_CLUSTER, 
     &           ibcksz2, minsize, topdiag
      INTEGER(8) :: JJ3
      INTEGER    :: K1RHS, K2RHS, JFirstRHS
      COMPLEX ZERO
      parameter( zero = (0.0e0,0.0e0) )
      nbcolf  = iw(ioldps+keep(ixsz))
      nbrowf  = iw(ioldps+2+keep(ixsz))
      nass    = iw(ioldps+1+keep(ixsz))
      nslaves= iw(ioldps+5 + keep(ixsz))
      hf      = 6 + nslaves + keep(ixsz)
!$    NOMP = OMP_GET_MAX_THREADS()
      IF (keep(50) .EQ. 0 .OR. nbrowf .LT. keep(63)) THEN
!$      CHUNK8 = int(KEEP(361),8)
!$OMP   PARALLEL DO PRIVATE(JJ8) SCHEDULE(STATIC, CHUNK8)
!$OMP&  IF (int(NBROWF,8)*int(NBCOLF,8) > int(KEEP(361),8)
!$OMP&    .AND. NOMP .GT. 1)
        DO jj8=poselt, poselt+int(nbrowf,8)*int(nbcolf,8)-1_8
          a(jj8) = zero
        ENDDO
!$omp   END PARALLEL do
      ELSE
        topdiag = 0
        IF (iw(ioldps+xxlr).GE.1) THEN
          CALL get_cut(iw(ioldps+hf:ioldps+hf+nbrowf-1), 0,
     &                    nbrowf, lrgroups, npartscb, 
     &                    npartsass, begs_blr_ls)
          nb_blr_ls = npartscb
          call max_cluster(begs_blr_ls,nb_blr_ls+1,maxi_cluster)
          DEALLOCATE(begs_blr_ls)
          CALL compute_blr_vcs(keep(472), ibcksz2, keep(488), nass)
          minsize = int(ibcksz2 / 2)
          topdiag = max(2*minsize + maxi_cluster-1, topdiag)
        ENDIF
!$      CHUNK = max( KEEP(360)/2,
!$   &               ((nbrowf+nomp-1)/nomp +2) / 3 )
!$OMP   PARALLEL DO PRIVATE(APOS,JJ3,JJ8) SCHEDULE(STATIC,CHUNK)
!$OMP&  IF (NBROWF .GT. KEEP(360) .AND. NOMP .GT. 1)
        DO jj8 = 0_8, int(nbrowf-1,8)
          apos = poselt+ jj8*int(nbcolf,8)
          jj3  = min( int(nbcolf,8)  - 1_8, 
     &           jj8 + int(nbcolf-nbrowf,8) + topdiag )
          a(apos: apos+jj3) = zero
        ENDDO
!$OMP   END PARALLEL DO
      ENDIF
          k1 = ioldps + hf + nbrowf
          k2 = k1 + nbcolf - 1
          jpos = 1
          DO k = k1, k2
           j        = iw(k)
           itloc(j) = -jpos
           jpos     = jpos + 1
          END DO
          k1 = ioldps + hf 
          k2 = k1 + nbrowf - 1
          jpos = 1
          IF ((keep(253).GT.0).AND.(keep(50).NE.0)) THEN
           k1rhs = 0
           k2rhs = -1
           DO k = k1, k2
            j        = iw(k)
            itloc(j) = -itloc(j)*nbcolf + jpos
            IF ((k1rhs.EQ.0).AND.(j.GT.n)) THEN
             k1rhs = k
             jfirstrhs=j-n 
            ENDIF
            jpos     = jpos + 1
           ENDDO
           IF (k1rhs.GT.0) k2rhs=k2
           IF ( k2rhs.GE.k1rhs ) THEN
             in = inode
             DO WHILE (in.GT.0) 
               ijrow = -itloc(in)  
               DO k = k1rhs, k2rhs
                j    = iw(k)       
                i    = itloc(j)    
                iloc = mod(i,nbcolf) 
              apos = poselt+int(iloc-1,8)*int(nbcolf,8) + 
     &               int(ijrow-1,8) 
              a(apos) = a(apos) + rhs_mumps(
     &                 (jfirstrhs+(k-k1rhs)-1)*keep(254)+ in)
             ENDDO
             in = fils(in)
            ENDDO
           ENDIF
          ELSE  
           DO k = k1, k2
            j        = iw(k)
            itloc(j) = -itloc(j)*nbcolf + jpos
            jpos     = jpos + 1
           END DO
          ENDIF
          elbeg  = frt_ptr(inode)
          numelt = frt_ptr(inode+1) - elbeg
          DO iell=elbeg,elbeg+numelt-1
           elti = frt_elt(iell)
           j18= ptraiw(elti)
           j28= ptraiw(elti+1)-1_8
           aii8 = ptrarw(elti)
           size_elti8 = j28 - j18 + 1_8
           DO ii8=j18,j28
            i = itloc(intarr(ii8))
            IF (keep(50).EQ.0) THEN
             IF (i.LE.0) cycle
             ainput8    = aii8 + ii8 - j18
             ipos = mod(i,nbcolf)
             ict12 = poselt + int(ipos-1,8) * int(nbcolf,8)
             DO jj8 = j18, j28
              jpos = itloc(intarr(jj8))
              IF (jpos.LE.0) THEN 
                   jpos = -jpos
              ELSE
                   jpos = jpos/nbcolf
              END IF
              apos2    = ict12 + int(jpos - 1,8)
              a(apos2) = a(apos2) +  dblarr(ainput8)
              ainput8   = ainput8 + size_elti8
             END DO
            ELSE
              IF ( i .EQ. 0 ) THEN 
               aii8 = aii8 + j28 - ii8 + 1_8
               cycle
              ENDIF
              IF ( i .LE. 0 ) THEN 
               ipos1 = -i
               ipos2 = 0
              ELSE 
               ipos1 = i/nbcolf
               ipos2 = mod(i,nbcolf)
              END IF
              ict12 =  poselt + int(ipos2-1,8)*int(nbcolf,8)
              DO jj8=ii8,j28
               aii8 = aii8 + 1_8
               j = itloc(intarr(jj8))
               IF ( j .EQ. 0 ) cycle
               IF ( ipos2.EQ.0 .AND. j.LE.0) cycle
               IF ( j .LE. 0 ) THEN
                jpos = -j
               ELSE
                jpos = j/nbcolf
               END IF
               IF ( (ipos1.GE.jpos) .AND. (ipos2.GT.0) ) THEN
                 apos2 = ict12  + int(jpos - 1,8)
                 a(apos2) = a(apos2) +  dblarr(aii8-1_8)
               END IF
               IF ( (ipos1.LT.jpos) .AND. (j.GT.0) ) THEN
                 ipos = mod(j,nbcolf)
                 jpos = ipos1
                 apos2 = poselt + int(ipos-1,8)*int(nbcolf,8)
     &                          + int(jpos - 1,8)
                 a(apos2) = a(apos2) +  dblarr(aii8-1_8)
               END IF
              END DO
            END IF
           END DO
          END DO
          k1 = ioldps + hf + nbrowf
          k2 = k1 + nbcolf - 1
          DO k = k1, k2
           j = iw(k)
           itloc(j) = 0
          END DO
      END SUBROUTINE cmumps_asm_slave_elements