zsol_matvec.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_mv_elt( N, NELT, ELTPTR, ELTVAR, A_ELT,
     &                          X, Y, K50, MTYPE )
      IMPLICIT NONE
C
C  Purpose
C  =======
C
C  To perform the matrix vector product
C      A_ELT X = Y    if MTYPE = 1
C      A_ELT^T X = Y  if MTYPE = 0
C
C  If K50 is different from 0, then the elements are
C  supposed to be in symmetric packed storage; the
C  lower part is stored by columns.
C  Otherwise, the element is square, stored by columns.
C
C  Note
C  ====
C
C  A_ELT is processed entry by entry and this code is not
C  optimized. In particular, one could gather/scatter
C  X / Y for each element to improve performance.
C
C  Arguments
C  =========
C
      INTEGER N, NELT, K50, MTYPE
      INTEGER ELTPTR( NELT + 1 ), ELTVAR( * )
      COMPLEX(kind=8) A_ELT( * ), X( N ), Y( N )
C
C  Local variables
C  ===============
C
      INTEGER IEL, I , J, SIZEI, IELPTR
      INTEGER(8) :: K8
      COMPLEX(kind=8) TEMP
      COMPLEX(kind=8) ZERO
      parameter( zero = (0.0d0,0.0d0) )
C
C
C     Executable statements
C     =====================
C
      y = zero
      k8 = 1_8
C     --------------------
C     Process the elements
C     --------------------
      DO iel = 1, nelt
        sizei  = eltptr( iel + 1 ) - eltptr( iel )
        ielptr = eltptr( iel ) - 1
        IF ( k50 .eq. 0 ) THEN
C         -------------------
C         Unsymmetric element
C         stored by columns
C         -------------------
          IF ( mtype .eq. 1 ) THEN
C           -----------------
C           Compute A_ELT x X
C           -----------------
            DO j = 1, sizei
              temp = x( eltvar( ielptr + j ) )
              DO i = 1, sizei
                y( eltvar( ielptr + i ) ) =
     &          y( eltvar( ielptr + i ) ) +
     &             a_elt( k8 ) * temp
                k8 = k8 + 1
              END DO
            END DO
          ELSE
C           -------------------
C           Compute A_ELT^T x X
C           -------------------
            DO j = 1, sizei
              temp = y( eltvar( ielptr + j ) )
              DO i = 1, sizei
                temp = temp + 
     &          a_elt( k8 ) * x( eltvar( ielptr + i ) )
                k8 = k8 + 1
              END DO
              y( eltvar( ielptr + j ) ) = temp
            END DO
          END IF
        ELSE
C         -----------------
C         Symmetric element
C         L stored by cols
C         -----------------
          DO j = 1, sizei
C           Diagonal counted once
            y( eltvar( ielptr + j ) ) =
     &      y( eltvar( ielptr + j ) ) +
     &           a_elt( k8 ) * x( eltvar( ielptr + j ) )
            k8 = k8 + 1
            DO i = j+1, sizei
C             Off diagonal + transpose
              y( eltvar( ielptr + i ) ) =
     &        y( eltvar( ielptr + i ) ) +
     &           a_elt( k8 ) * x( eltvar( ielptr + j ) )
              y( eltvar( ielptr + j ) ) =
     &        y( eltvar( ielptr + j ) ) +
     &           a_elt( k8 ) * x( eltvar( ielptr + i ) )
              k8 = k8 + 1
            END DO
          END DO
        END IF
      END DO
      RETURN
      END SUBROUTINE zmumps_mv_elt
      SUBROUTINE zmumps_loc_mv8
     &( n, nz_loc8, irn_loc, jcn_loc, a_loc, x, y_loc,
     &  ldlt, mtype)
      IMPLICIT NONE
C
C     Purpose:
C     =======
C
C     Perform a distributed matrix vector product.
C        Y_loc <- A X   if MTYPE = 1
C        Y_loc <- A^T X if MTYPE = 0
C
C     Notes:
C     =====
C
C     1) assembly of all Y_loc still has to be done on exit.
C     2) X should be available on all processors.
C
C     Arguments:
C     =========
C
      INTEGER N
      INTEGER(8) :: NZ_loc8
      INTEGER IRN_loc( NZ_loc8 ), JCN_loc( NZ_loc8 )
      COMPLEX(kind=8) A_loc( NZ_loc8 ), X( N ), Y_loc( N )
      INTEGER LDLT, MTYPE
C
C     Locals variables:
C     ================
C
      INTEGER I, J
      INTEGER(8) :: K8
      COMPLEX(kind=8) ZERO
      parameter( zero = (0.0d0,0.0d0) )
      y_loc = zero
      IF ( ldlt .eq. 0 ) THEN
C       Unsymmetric
        IF ( mtype .eq. 1 ) THEN
C         No transpose
          DO k8 = 1_8, nz_loc8
            i = irn_loc(k8)
            j = jcn_loc(k8)
            IF ((i .LE. 0) .OR. (i .GT. n) .OR.
     &          (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
          y_loc(i) = y_loc(i) + a_loc(k8) * x(j)
        ENDDO
        ELSE
C         Transpose
          DO k8 = 1_8, nz_loc8
            i = irn_loc(k8)
            j = jcn_loc(k8)
            IF ((i .LE. 0) .OR. (i .GT. n)
     &        .OR. (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
          y_loc(j) = y_loc(j) + a_loc(k8) * x(i)
        ENDDO
        END IF
      ELSE
C       Lower (or upper) part of symmetric
C       matrix was provided (LDLT facto)
        DO k8 = 1_8, nz_loc8
          i = irn_loc(k8)
          j = jcn_loc(k8)
          IF ((i .LE. 0) .OR. (i .GT. n) .OR.
     &        (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
          y_loc(i) = y_loc(i) + a_loc(k8) * x(j)
          IF (j.NE.i) THEN
            y_loc(j) = y_loc(j) + a_loc(k8) * x(i)
          ENDIF
        ENDDO
      END IF
      RETURN
      END SUBROUTINE zmumps_loc_mv8
      SUBROUTINE zmumps_mv8( N, NZ8, IRN, ICN, ASPK, X, Y,
     &                      LDLT, MTYPE, MAXTRANS, PERM,
     &                      IFLAG, IERROR )
C
C     Purpose:
C     =======
C
C     Perform matrix-vector product
C        Y <- A X if MTYPE = 1
C        Y <- A^T X if MTYPE = 0
C
C
C     Note:
C     ====
C
C     MAXTRANS should be set to 1 if a column permutation
C     was applied on A and we still want the matrix vector
C     product wrt the original matrix.
C
C     Arguments:
C     =========
C
      INTEGER N, LDLT, MTYPE, MAXTRANS
      INTEGER(8) :: NZ8
      INTEGER IRN( NZ8 ), ICN( NZ8 ) 
      INTEGER PERM( N )
      COMPLEX(kind=8) ASPK( NZ8 ), X( N ), Y( N )
      INTEGER, intent(inout) :: IFLAG, IERROR
C
C     Local variables
C     ===============
C
      INTEGER I, J
      INTEGER(8) :: K8
      COMPLEX(kind=8), DIMENSION(:), ALLOCATABLE :: PX
      COMPLEX(kind=8) ZERO
      INTEGER :: allocok
      parameter( zero = (0.0d0,0.0d0) )
      y = zero
      ALLOCATE(px(n), stat=allocok)
      IF (allocok < 0) THEN
        iflag  = -13
        ierror = n
        RETURN
      ENDIF
C
C     --------------------------------------
C     Permute X if A has been permuted
C     with some max-trans column permutation
C     --------------------------------------
      IF ( maxtrans .eq. 1 .and. mtype .eq. 1) THEN
        DO i = 1, n
          px(i) = x( perm( i ) )
        END DO
      ELSE
        px = x
      END IF
      IF ( ldlt .eq. 0 ) THEN
C
C     Complete unsymmetric matrix was provided (LU facto)
       IF (mtype .EQ. 1) THEN
        DO k8 = 1_8, nz8
          i = irn(k8)
          j = icn(k8)
          IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
          y(i) = y(i) + aspk(k8) * px(j)
        ENDDO
       ELSE
        DO k8 = 1_8, nz8
          i = irn(k8)
          j = icn(k8)
          IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
          y(j) = y(j) + aspk(k8) * px(i)
        ENDDO
       ENDIF
C
      ELSE
C
C       Lower (or upper) part of symmetric
C       matrix was provided (LDLT facto)
        DO k8 = 1_8, nz8
          i = irn(k8)
          j = icn(k8)
          IF ((i .LE. 0) .OR. (i .GT. n) .OR. (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
          y(i) = y(i) + aspk(k8) * px(j)
          IF (j.NE.i) THEN
            y(j) = y(j) + aspk(k8) * px(i)
          ENDIF
        ENDDO
      END IF
      IF ( maxtrans .EQ. 1 .AND. mtype .eq. 0 ) THEN
      px = y
      DO i = 1, n
        y( perm( i ) ) = px( i )
      END DO
      END IF
      DEALLOCATE(px)
      RETURN
      END SUBROUTINE zmumps_mv8
C
C
      SUBROUTINE zmumps_loc_omega1
     &( n, nz_loc8, irn_loc, jcn_loc, a_loc, x, y_loc,
     &  ldlt, mtype)
      IMPLICIT NONE
C
C     Purpose:
C     =======
C     Compute
C        * If MTYPE = 1
C            Y_loc(i) = Sum | Aij | | Xj |
C                        j
C        * If MTYPE = 0
C            Y_loc(j) = Sum | Aij | | Xi |
C
C
C     Notes:
C     =====
C
C     1) assembly of all Y_loc still has to be done.
C     2) X should be available on all processors.
C
C     Arguments:
C     =========
C
      INTEGER N
      INTEGER(8) :: NZ_loc8
      INTEGER IRN_loc( NZ_loc8 ), JCN_loc( NZ_loc8 )
      COMPLEX(kind=8) A_loc( NZ_loc8 ), X( N )
      DOUBLE PRECISION Y_loc( N )
      INTEGER LDLT, MTYPE
C
C     Local variables:
C     ===============
C
      INTEGER I, J
      INTEGER(8) :: K8
      DOUBLE PRECISION, PARAMETER :: RZERO=0.0d0
C
      y_loc = rzero
      IF ( ldlt .eq. 0 ) THEN
C       Unsymmetric
        IF ( mtype .eq. 1 ) THEN
C         No transpose
          DO k8 = 1_8, nz_loc8
            i = irn_loc(k8)
            j = jcn_loc(k8)
            IF ((i .LE. 0) .OR. (i .GT. n) .OR.
     &          (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
            y_loc(i) = y_loc(i) + abs( a_loc(k8) * x(j) )
          ENDDO
        ELSE
C         Transpose
          DO k8 = 1_8, nz_loc8
            i = irn_loc(k8)
            j = jcn_loc(k8)
            IF ((i .LE. 0) .OR. (i .GT. n)
     &        .OR. (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
          y_loc(j) = y_loc(j) + abs( a_loc(k8) * x(i) )
          ENDDO
        END IF
      ELSE
C       Lower (or upper) part of symmetric
C       matrix was provided (LDLT facto)
        DO k8 = 1_8, nz_loc8
          i = irn_loc(k8)
          j = jcn_loc(k8)
          IF ((i .LE. 0) .OR. (i .GT. n) .OR.
     &        (j .LE. 0) .OR. (j .GT. n)
     &        ) cycle
          y_loc(i) = y_loc(i) + abs( a_loc(k8) * x(j) )
          IF (j.NE.i) THEN
            y_loc(j) = y_loc(j) + abs( a_loc(k8) * x(i) )
          ENDIF
        ENDDO
      END IF
      RETURN
      END SUBROUTINE zmumps_loc_omega1