v8p2r0/html/dldu21_8f_source.html

 !                   *****************
                     SUBROUTINE dldu21
 !                   *****************
 !
      &(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 Q1 QUADRILATERALS.
 !+
 !+            REQUIRES THAT THE DIAGONAL OF A BE THE IDENTITY.
 !code
 !+            EACH ELEMENTARY MATRIX IS DECOMPOSED IN THE FORM:
 !+
 !+            LE * DE * 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
 !
 !warning  FOR NONSYMMETRICAL MATRICES: UE  LE
 !+
 !+  UE (THE BETAS) IS STORED IN XB (. , 1 TO  6)
 !+
 !+  LE (THE ALFAS) IS STORED IN XB (. , 7 TO 12)
 !
 !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_dldu21 => dldu21
 !
       USE declarations_special
       IMPLICIT NONE
 !
 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
 !
       INTEGER, INTENT(IN)           :: NELEM,NELMAX,LV,NPOIN
       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,4)
       LOGICAL, INTENT(IN)           :: COPY
 !
 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
 !
       INTEGER IELEM
 !
       DOUBLE PRECISION     A12,A13,A14
       DOUBLE PRECISION A21,A22,A23,A24
       DOUBLE PRECISION A31,A32,A33,A34
       DOUBLE PRECISION A41,A42,A43,A44
       DOUBLE PRECISION        BETA12,BETA13,BETA14
       DOUBLE PRECISION ALFA21,BETA22,BETA23,BETA24
       DOUBLE PRECISION ALFA31,ALFA32,BETA33,BETA34
       DOUBLE PRECISION ALFA41,ALFA42,ALFA43,BETA44
 !
 !-----------------------------------------------------------------------
 !
 ! REQUIRES THAT THE DIAGONAL OF A BE THE IDENTITY
 !
       IF(typdia(1:1).NE.'I'.AND.ncsize.LE.1) THEN
         WRITE(lu,1001) typdia(1:1)
 1001    FORMAT(1x,'DLDU21 (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*6)
 !
         DO ielem=1,nelem
 !
 ! MATRIX TO FACTORISE (SYMMETRICAL WITH 1S ON THE DIAGONAL)
 !
 ! LINE 1
 !         A11 = 1.D0
           a12 = xa(ielem,1)
           a13 = xa(ielem,2)
           a14 = xa(ielem,3)
 ! LINE 2
           a22 = 1.d0
           a23 = xa(ielem,4)
           a24 = xa(ielem,5)
 ! LINE 3
           a33 = 1.d0
           a34 = xa(ielem,6)
 ! LINE 4
           a44 = 1.d0
 !
 ! CROUT L*U FACTORISATION
 !
           alfa21 = a12
           alfa31 = a13
           alfa41 = a14
 !
           beta12 =  a12
           beta22 =  a22 - alfa21*beta12
           alfa32 = (a23 - alfa31*beta12)/beta22
           alfa42 = (a24 - alfa41*beta12)/beta22
 !
           beta13 =  a13
           beta23 =  a23 - alfa21*beta13
           beta33 =  a33 - alfa31*beta13 - alfa32*beta23
           alfa43 = (a34 - alfa41*beta13 - alfa42*beta23)/beta33
 !
           beta14 =  a14
           beta24 =  a24 - alfa21*beta14
           beta34 =  a34 - alfa31*beta14 - alfa32*beta24
           beta44 =  a44 - alfa41*beta14 - alfa42*beta24 - alfa43*beta34
 !
 ! STORES IN XB AND W2,...,W4
 !
           xb(ielem,1 ) = alfa21
           xb(ielem,2 ) = alfa31
           xb(ielem,3 ) = alfa41
           xb(ielem,4 ) = alfa32
           xb(ielem,5 ) = alfa42
           xb(ielem,6 ) = alfa43
 !
           w(ielem,2)    = beta22
           w(ielem,3)    = beta33
           w(ielem,4)    = beta44
 !
         ENDDO ! IELEM
 !
 !-----------------------------------------------------------------------
 !
       ELSEIF(typexa(1:1).EQ.'Q') THEN
 !
         IF(copy) CALL ov('X=Y     ', x=xb, y=xa, dim1=nelmax*12)
 !
         DO ielem=1,nelem
 !
 ! MATRIX TO FACTORISE (WITH 1S ON THE DIAGONAL)
 !
 !         A11 = 1.D0
           a22 = 1.d0
           a33 = 1.d0
           a44 = 1.d0
 !
           a12 = xa(ielem,1 )
           a13 = xa(ielem,2 )
           a14 = xa(ielem,3 )
           a23 = xa(ielem,4 )
           a24 = xa(ielem,5 )
           a34 = xa(ielem,6 )
 !
           a21 = xa(ielem,7 )
           a31 = xa(ielem,8 )
           a41 = xa(ielem,9 )
           a32 = xa(ielem,10)
           a42 = xa(ielem,11)
           a43 = xa(ielem,12)
 !
 ! CROUT L*U FACTORISATION
 !
           alfa21 = a21
           alfa31 = a31
           alfa41 = a41
 !
           beta12 =  a12
           beta22 =  a22 - alfa21*beta12
           alfa32 = (a32 - alfa31*beta12)/beta22
           alfa42 = (a42 - alfa41*beta12)/beta22
 !
           beta13 =  a13
           beta23 =  a23 - alfa21*beta13
           beta33 =  a33 - alfa31*beta13 - alfa32*beta23
           alfa43 = (a43 - alfa41*beta13 - alfa42*beta23)/beta33
 !
           beta14 =  a14
           beta24 =  a24 - alfa21*beta14
           beta34 =  a34 - alfa31*beta14 - alfa32*beta24
           beta44 =  a44 - alfa41*beta14 - alfa42*beta24 - alfa43*beta34
 !
 ! STORES IN XB AND W2,...,W4
 ! L D U FACTORISATION AT THE SAME TIME
 !
           xb(ielem,1 ) = beta12
           xb(ielem,2 ) = beta13
           xb(ielem,3 ) = beta14
           xb(ielem,4 ) = beta23/beta22
           xb(ielem,5 ) = beta24/beta22
           xb(ielem,6 ) = beta34/beta33
 !
           xb(ielem,07) = alfa21
           xb(ielem,08) = alfa31
           xb(ielem,09) = alfa41
           xb(ielem,10) = alfa32
           xb(ielem,11) = alfa42
           xb(ielem,12) = alfa43
 !
           w(ielem,2)    = beta22
           w(ielem,3)    = beta33
           w(ielem,4)    = beta44
 !
         ENDDO ! IELEM
 !
 !-----------------------------------------------------------------------
 !
       ELSE
         WRITE(lu,2001) typexa(1:1)
 2001    FORMAT(1x,'DLDU21 (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,3,w(1,2),.true.,lv)
 !
 !  INVERTS DB
 !
       CALL ov('X=1/Y   ', x=db, y=db, dim1=npoin)
 !
 !-----------------------------------------------------------------------
 !
       RETURN
       END
ov
subroutine ov(OP, X, Y, Z, C, DIM1)
Definition: ov.f:7

declarations_special
Definition: declarations_special.F:3

asmvec
subroutine asmvec(X, IKLE, NPOIN, NELEM, NELMAX, NDP, W, INIT, LV)
Definition: asmvec.f:7

dldu21
subroutine dldu21(DB, XB, TYPDIA, XA, TYPEXA, IKLE, NELEM, NELMAX, NPOIN, W, COPY, LV)
Definition: dldu21.f:7

declarations_special::lu
integer lu
Definition: declarations_special.F:45

bief
Definition: bief.f:3