mer41.f

Go to the documentation of this file.

!                   ****************
                    SUBROUTINE mer41
!                   ****************
!
     &(x, xa1 ,xa2 ,xa3 ,xa4 ,xa5 ,
     &    xa6 ,xa7 ,xa8 ,xa9 ,xa10,
     &    xa11,xa12,xa13,xa14,xa15,
     &    ikle1,ikle2,ikle3,ikle4,ikle5,ikle6,
     &    nelem,nelmax,npoin,lv)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    PRODUCT X = U B (BEWARE: ELEMENT BY ELEMENT).
!+
!+            HERE: P1 PRISM ELEMENT OR ELEMENT WITH 6 POINTS.
!+
!+            REVERSE OPERATION FROM THAT IN SUBROUTINE REMONT,
!+                HENCE THE NAME.
!code
!+        THE MATRIX U IS HERE THE RESULT OF A DECOMPOSITION
!+        DONE IN SUBROUTINE DECLDU.
!+
!+        EACH ELEMENTARY MATRIX HAS BEEN FACTORISED IN THE FORM :
!+
!+        LE X DE X UE
!+
!+        LE : LOWER TRIANGULAR WITH 1S ON THE DIAGONAL
!+        DE : DIAGONAL
!+        UE : UPPER TRIANGULAR WITH 1S ON THE DIAGONAL
!+
!+                                            T
!+        IF THE MATRIX IS SYMMETRICAL : LE =  UE
!
!history  J-M HERVOUET (LNH)    ; F LEPEINTRE (LNH)
!+        05/02/91
!+        V5P1
!+
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        13/07/2010
!+        V6P0
!+   Translation of French comments within the FORTRAN sources into
!+   English comments
!
!history  N.DURAND (HRW), S.E.BOURBAN (HRW)
!+        21/08/2010
!+        V6P0
!+   Creation of DOXYGEN tags for automated documentation and
!+   cross-referencing of the FORTRAN sources
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| IKLE1          |-->| FIRST POINTS OF PRISMS
!| IKLE2          |-->| SECOND POINTS OF PRISMS
!| IKLE3          |-->| THIRD POINTS OF PRISMS
!| IKLE4          |-->| FOUTH POINTS OF PRISMS
!| IKLE5          |-->| FIFTH POINTS OF PRISMS
!| IKLE6          |-->| SIXTH POINTS OF PRISMS
!| LV             |-->| VECTOR LENGTH OF THE MACHINE
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| X              |<--| RESULTING VECTOR
!| XA10           |-->| OFF-DIAGONAL TERM
!| XA11           |-->| OFF-DIAGONAL TERM
!| XA12           |-->| OFF-DIAGONAL TERM
!| XA13           |-->| OFF-DIAGONAL TERM
!| XA14           |-->| OFF-DIAGONAL TERM
!| XA15           |-->| OFF-DIAGONAL TERM
!| XA2            |-->| OFF-DIAGONAL TERM
!| XA3            |-->| OFF-DIAGONAL TERM
!| XA4            |-->| OFF-DIAGONAL TERM
!| XA5            |-->| OFF-DIAGONAL TERM
!| XA6            |-->| OFF-DIAGONAL TERM
!| XA7            |-->| OFF-DIAGONAL TERM
!| XA8            |-->| OFF-DIAGONAL TERM
!| XA9            |-->| OFF-DIAGONAL TERM
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN) :: NPOIN,NELEM,NELMAX,LV
      INTEGER, INTENT(IN) :: IKLE1(nelmax),IKLE2(nelmax),IKLE3(nelmax)
      INTEGER, INTENT(IN) :: IKLE4(nelmax),IKLE5(nelmax),IKLE6(nelmax)
!
      DOUBLE PRECISION, INTENT(INOUT) :: X(npoin)
      DOUBLE PRECISION, INTENT(IN) :: XA1(nelmax),XA2(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA3(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA4(nelmax),XA5(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA6(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA7(nelmax),XA8(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA9(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA10(nelmax),XA11(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA12(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA13(nelmax),XA14(nelmax)
      DOUBLE PRECISION, INTENT(IN) :: XA15(nelmax)
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM,IB
!
      INTRINSIC min
!
!-----------------------------------------------------------------------
!
! RESUMES INVERSION OF THE LOWER TRIANGULAR MATRICES
!
!-----------------------------------------------------------------------
! LOOP IN SCALAR MODE (LV=1) OR WITH FORCED VECTORISATION
!-----------------------------------------------------------------------
!
      IF(lv.EQ.1) THEN
!
!  SCALAR MODE
!
      DO ielem = 1 , nelem
          x(ikle1(ielem))=x(ikle1(ielem))+xa5(ielem)*x(ikle6(ielem))
     &                                   +xa4(ielem)*x(ikle5(ielem))
     &                                   +xa3(ielem)*x(ikle4(ielem))
     &                                   +xa2(ielem)*x(ikle3(ielem))
     &                                   +xa1(ielem)*x(ikle2(ielem))
          x(ikle2(ielem))=x(ikle2(ielem))+xa9(ielem)*x(ikle6(ielem))
     &                                   +xa8(ielem)*x(ikle5(ielem))
     &                                   +xa7(ielem)*x(ikle4(ielem))
     &                                   +xa6(ielem)*x(ikle3(ielem))
          x(ikle3(ielem))=x(ikle3(ielem))+xa12(ielem)*x(ikle6(ielem))
     &                                   +xa11(ielem)*x(ikle5(ielem))
     &                                   +xa10(ielem)*x(ikle4(ielem))
          x(ikle4(ielem))=x(ikle4(ielem))+xa14(ielem)*x(ikle6(ielem))
     &                                   +xa13(ielem)*x(ikle5(ielem))
          x(ikle5(ielem))=x(ikle5(ielem))+xa15(ielem)*x(ikle6(ielem))
      ENDDO ! IELEM
!
      ELSE
!
!  VECTOR MODE
!
      DO ib = 1,(nelem+lv-1)/lv
!VOCL LOOP,NOVREC
!DIR$ IVDEP
      DO ielem = 1+(ib-1)*lv , min(nelem,ib*lv)
          x(ikle1(ielem))=x(ikle1(ielem))+xa5(ielem)*x(ikle6(ielem))
     &                                   +xa4(ielem)*x(ikle5(ielem))
     &                                   +xa3(ielem)*x(ikle4(ielem))
     &                                   +xa2(ielem)*x(ikle3(ielem))
     &                                   +xa1(ielem)*x(ikle2(ielem))
          x(ikle2(ielem))=x(ikle2(ielem))+xa9(ielem)*x(ikle6(ielem))
     &                                   +xa8(ielem)*x(ikle5(ielem))
     &                                   +xa7(ielem)*x(ikle4(ielem))
     &                                   +xa6(ielem)*x(ikle3(ielem))
          x(ikle3(ielem))=x(ikle3(ielem))+xa12(ielem)*x(ikle6(ielem))
     &                                   +xa11(ielem)*x(ikle5(ielem))
     &                                   +xa10(ielem)*x(ikle4(ielem))
          x(ikle4(ielem))=x(ikle4(ielem))+xa14(ielem)*x(ikle6(ielem))
     &                                   +xa13(ielem)*x(ikle5(ielem))
          x(ikle5(ielem))=x(ikle5(ielem))+xa15(ielem)*x(ikle6(ielem))
      ENDDO ! IELEM
      ENDDO ! IB
!
      ENDIF
!
!-----------------------------------------------------------------------
!
      RETURN
      END