dldu11.f

Go to the documentation of this file.

!                   *****************
                    SUBROUTINE dldu11
!                   *****************
!
     &(db,xb,typdia,xa,typexa,ikle,nelem,nelmax,npoin,w,copy,lv)
!
!***********************************************************************
! BIEF   V6P1                                   21/08/2010
!***********************************************************************
!
!brief    L D U FACTORISATION OF THE ELEMENTARY MATRICES
!+                IN MATRIX A
!+                FOR P1 TRIANGLES.
!+
!+            REQUIRES THAT THE DIAGONAL OF A BE THE IDENTITY.
!code
!+            EACH ELEMENTARY MATRIX IS DECOMPOSED 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
!+
!+            "DE" MATRICES ARE CONSIDERED LIKE DIAGONALS OF SIZE
!+            NPOIN X NPOIN, WHICH ARE FILLED WITH 1S FOR THE POINTS
!+            WHICH DO NOT BELONG TO THE CONSIDERED ELEMENT
!+
!+            THEN PERFORMS THE PRODUCT OF ALL THESE DIAGONALS
!+            YIELDING DIAGONAL DB
!
!history  J-M HERVOUET (LNH)    ; F LEPEINTRE (LNH)
!+        24/04/97
!+        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
!
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!| COPY           |-->| IF .TRUE. A IS COPIED INTO B.
!|                |   | IF .FALSE. B IS CONSIDERED ALREADY INITIALISED
!| DB             |<--| DIAGONAL OF MATRIX B
!| IKLE           |-->| CONNECTIVITY TABLE
!| LV             |-->| VECTOR LENGTH OF THE COMPUTER
!| NELEM          |-->| NUMBER OF ELEMENTS
!| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
!| NPOIN          |-->| NUMBER OF POINTS
!| TYPDIA         |<--| TYPE OF DIAGONAL ( 'Q', 'I' , OR '0' )
!| TYPEXA         |<--| TYPE OF OFF-DIAGONAL TERMS ('Q','S',OR '0')
!| W              |-->| WORK ARRAY OF DIMENSION (NELMAX,3)
!| XA             |<--| OFF-DIAGONAL TERMS OF MATRIX A
!| XB             |<--| OFF-DIAGONAL TERMS OF MATRIX B
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!
      USE bief, ex_dldu11 => dldu11
!
      USE declarations_special
      IMPLICIT NONE
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER, INTENT(IN)           :: NELEM,NELMAX,LV,NPOIN
!     NO DEFAULT INITIALISATION FOR USER TYPE COMPONENTS ALLOWED
      DOUBLE PRECISION, INTENT(INOUT) :: DB(npoin),XB(nelmax,*)
      DOUBLE PRECISION, INTENT(IN)  :: XA(nelmax,*)
      CHARACTER(LEN=1), INTENT(IN)  :: TYPDIA,TYPEXA
      INTEGER, INTENT(IN)           :: IKLE(nelmax,*)
      DOUBLE PRECISION, INTENT(OUT) :: W(nelmax,3)
      LOGICAL, INTENT(IN)           :: COPY
!
!+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
!
      INTEGER IELEM
!
!-----------------------------------------------------------------------
!
! REQUIRES THAT THE DIAGONAL OF A BE THE IDENTITY (EXCEPT IN PARALLEL MODE)
!
      IF(typdia(1:1).NE.'I'.AND.ncsize.LE.1) THEN
        WRITE(lu,101) typdia(1:1)
101     FORMAT(1x,'DLDU11 (BIEF) : DIAGONAL OF A NOT EQUAL TO I :',a1)
        CALL plante(0)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
      IF(typexa(1:1).EQ.'S') THEN
!
        IF(copy) CALL ov('X=Y     ', x=xb, y=xa, dim1=nelmax*3)
!
        DO ielem = 1 , nelem
          w(ielem,2) = 1.d0 - xb(ielem,1)**2
          xb(ielem,3) = (xb(ielem,3)-xb(ielem,1)*xb(ielem,2))/w(ielem,2)
          w(ielem,3) = 1.d0 - xb(ielem,2)**2 -xb(ielem,3)**2
        ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
      ELSEIF(typexa(1:1).EQ.'Q') THEN
!
        IF(copy) CALL ov('X=Y     ', x=xb, y=xa, dim1=nelmax*6)
!
        DO ielem = 1 , nelem
! L U FACTORISATION
          w(ielem,2)=1.d0 - xb(ielem,1)*xb(ielem,4)
          xb(ielem,6) = (xb(ielem,6)-xb(ielem,1)*xb(ielem,5))/w(ielem,2)
          xb(ielem,3) =  xb(ielem,3)-xb(ielem,4)*xb(ielem,2)
          w(ielem,3)  = 1.d0-xb(ielem,2)*xb(ielem,5)
     &                - xb(ielem,3)*xb(ielem,6)
! L D U FACTORISATION
          xb(ielem,3) = xb(ielem,3) / w(ielem,2)
        ENDDO ! IELEM
!
!-----------------------------------------------------------------------
!
      ELSE
        WRITE(lu,201) typexa(1:1)
201     FORMAT(1x,'DLDU11 (BIEF) : TYPE OF MATRIX NOT AVAILABLE :',a1)
        CALL plante(0)
        stop
      ENDIF
!
!-----------------------------------------------------------------------
!
!  MULTIPLICATIVE ASSEMBLY OF THE DIAGONAL WITH INITIALISATION OF DB TO 1
!  SKIPS IKLE1 BECAUSE W1 = 1
!
      CALL asmvec(db,ikle(1,2),npoin,nelem,nelmax,2,w(1,2),.true.,lv)
!
!  INVERTS DB
!
      CALL ov('X=1/Y   ', x=db, y=db, dim1=npoin)
!
!-----------------------------------------------------------------------
!
      RETURN
      END