dldu21.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\dldu21.f
00002 !
00089                      SUBROUTINE DLDU21
00090 !                    *****************
00091 !
00092      &(DB,XB,TYPDIA,XA,TYPEXA,IKLE,NELEM,NELMAX,NPOIN,W,COPY,LV)
00093 !
00094 !***********************************************************************
00095 ! BIEF   V6P1                                   21/08/2010
00096 !***********************************************************************
00097 !
00098 !
00099 !
00100 !
00101 !
00102 !
00103 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00104 !| COPY           |-->| IF .TRUE. A IS COPIED INTO B.
00105 !|                |   | IF .FALSE. B IS CONSIDERED ALREADY INITIALISED
00106 !| DB             |<--| DIAGONAL OF MATRIX B
00107 !| IKLE           |-->| CONNECTIVITY TABLE
00108 !| LV             |-->| VECTOR LENGTH OF THE COMPUTER
00109 !| NELEM          |-->| NUMBER OF ELEMENTS
00110 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00111 !| NPOIN          |-->| NUMBER OF POINTS
00112 !| TYPDIA         |<--| TYPE OF DIAGONAL ( 'Q', 'I' , OR '0' )
00113 !| TYPEXA         |<--| TYPE OF OFF-DIAGONAL TERMS ('Q','S',OR '0')
00114 !| W              |-->| WORK ARRAY OF DIMENSION (NELMAX,3)
00115 !| XA             |<--| OFF-DIAGONAL TERMS OF MATRIX A
00116 !| XB             |<--| OFF-DIAGONAL TERMS OF MATRIX B
00117 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00118 !
00119       USE BIEF, EX_DLDU21 => DLDU21
00120 !
00121       IMPLICIT NONE
00122       INTEGER LNG,LU
00123       COMMON/INFO/LNG,LU
00124 !
00125 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00126 !
00127       INTEGER, INTENT(IN)           :: NELEM,NELMAX,LV,NPOIN
00128       DOUBLE PRECISION, INTENT(OUT) :: DB(NPOIN),XB(NELMAX,*)
00129       DOUBLE PRECISION, INTENT(IN)  :: XA(NELMAX,*)
00130       CHARACTER(LEN=1), INTENT(IN)  :: TYPDIA,TYPEXA
00131       INTEGER, INTENT(IN)           :: IKLE(NELMAX,*)
00132       DOUBLE PRECISION, INTENT(OUT) :: W(NELMAX,4)
00133       LOGICAL, INTENT(IN)           :: COPY
00134 !
00135 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00136 !
00137       INTEGER IELEM
00138 !
00139       DOUBLE PRECISION Z(1),C
00140 !
00141       DOUBLE PRECISION     A12,A13,A14
00142       DOUBLE PRECISION A21,A22,A23,A24
00143       DOUBLE PRECISION A31,A32,A33,A34
00144       DOUBLE PRECISION A41,A42,A43,A44
00145       DOUBLE PRECISION        BETA12,BETA13,BETA14
00146       DOUBLE PRECISION ALFA21,BETA22,BETA23,BETA24
00147       DOUBLE PRECISION ALFA31,ALFA32,BETA33,BETA34
00148       DOUBLE PRECISION ALFA41,ALFA42,ALFA43,BETA44
00149 !
00150 !-----------------------------------------------------------------------
00151 !
00152 ! REQUIRES THAT THE DIAGONAL OF A BE THE IDENTITY
00153 !
00154       IF(TYPDIA(1:1).NE.'I'.AND.NCSIZE.LE.1) THEN
00155         IF (LNG.EQ.1) WRITE(LU,1000) TYPDIA(1:1)
00156         IF (LNG.EQ.2) WRITE(LU,1001) TYPDIA(1:1)
00157 1000    FORMAT(1X,'DLDU21 (BIEF) : DIAGONALE DE A NON EGALE A I :',A1)
00158 1001    FORMAT(1X,'DLDU21 (BIEF) : DIAGONAL OF A NOT EQUAL TO I :',A1)
00159         CALL PLANTE(0)
00160         STOP
00161       ENDIF
00162 !
00163 !-----------------------------------------------------------------------
00164 !
00165       IF(TYPEXA(1:1).EQ.'S') THEN
00166 !
00167         IF(COPY) CALL OV('X=Y     ' , XB , XA , Z , C  , NELMAX*6 )
00168 !
00169         DO IELEM=1,NELEM
00170 !
00171 ! MATRIX TO FACTORISE (SYMMETRICAL WITH 1S ON THE DIAGONAL)
00172 !
00173 ! LINE 1
00174 !          A11 = 1.D0
00175            A12 = XA(IELEM,1)
00176            A13 = XA(IELEM,2)
00177            A14 = XA(IELEM,3)
00178 ! LINE 2
00179            A22 = 1.D0
00180            A23 = XA(IELEM,4)
00181            A24 = XA(IELEM,5)
00182 ! LINE 3
00183            A33 = 1.D0
00184            A34 = XA(IELEM,6)
00185 ! LINE 4
00186            A44 = 1.D0
00187 !
00188 ! CROUT L*U FACTORISATION
00189 !
00190            ALFA21 = A12
00191            ALFA31 = A13
00192            ALFA41 = A14
00193 !
00194            BETA12 =  A12
00195            BETA22 =  A22 - ALFA21*BETA12
00196            ALFA32 = (A23 - ALFA31*BETA12)/BETA22
00197            ALFA42 = (A24 - ALFA41*BETA12)/BETA22
00198 !
00199            BETA13 =  A13
00200            BETA23 =  A23 - ALFA21*BETA13
00201            BETA33 =  A33 - ALFA31*BETA13 - ALFA32*BETA23
00202            ALFA43 = (A34 - ALFA41*BETA13 - ALFA42*BETA23)/BETA33
00203 !
00204            BETA14 =  A14
00205            BETA24 =  A24 - ALFA21*BETA14
00206            BETA34 =  A34 - ALFA31*BETA14 - ALFA32*BETA24
00207            BETA44 =  A44 - ALFA41*BETA14 - ALFA42*BETA24 - ALFA43*BETA34
00208 !
00209 ! STORES IN XB AND W2,...,W4
00210 !
00211            XB(IELEM,1 ) = ALFA21
00212            XB(IELEM,2 ) = ALFA31
00213            XB(IELEM,3 ) = ALFA41
00214            XB(IELEM,4 ) = ALFA32
00215            XB(IELEM,5 ) = ALFA42
00216            XB(IELEM,6 ) = ALFA43
00217 !
00218            W(IELEM,2)    = BETA22
00219            W(IELEM,3)    = BETA33
00220            W(IELEM,4)    = BETA44
00221 !
00222         ENDDO ! IELEM
00223 !
00224 !-----------------------------------------------------------------------
00225 !
00226       ELSEIF(TYPEXA(1:1).EQ.'Q') THEN
00227 !
00228         IF(COPY) CALL OV('X=Y     ' , XB , XA , Z , C  , NELMAX*12 )
00229 !
00230         DO IELEM=1,NELEM
00231 !
00232 ! MATRIX TO FACTORISE (WITH 1S ON THE DIAGONAL)
00233 !
00234 !          A11 = 1.D0
00235            A22 = 1.D0
00236            A33 = 1.D0
00237            A44 = 1.D0
00238 !
00239            A12 = XA(IELEM,1 )
00240            A13 = XA(IELEM,2 )
00241            A14 = XA(IELEM,3 )
00242            A23 = XA(IELEM,4 )
00243            A24 = XA(IELEM,5 )
00244            A34 = XA(IELEM,6 )
00245 !
00246            A21 = XA(IELEM,7 )
00247            A31 = XA(IELEM,8 )
00248            A41 = XA(IELEM,9 )
00249            A32 = XA(IELEM,10)
00250            A42 = XA(IELEM,11)
00251            A43 = XA(IELEM,12)
00252 !
00253 ! CROUT L*U FACTORISATION
00254 !
00255            ALFA21 = A21
00256            ALFA31 = A31
00257            ALFA41 = A41
00258 !
00259            BETA12 =  A12
00260            BETA22 =  A22 - ALFA21*BETA12
00261            ALFA32 = (A32 - ALFA31*BETA12)/BETA22
00262            ALFA42 = (A42 - ALFA41*BETA12)/BETA22
00263 !
00264            BETA13 =  A13
00265            BETA23 =  A23 - ALFA21*BETA13
00266            BETA33 =  A33 - ALFA31*BETA13 - ALFA32*BETA23
00267            ALFA43 = (A43 - ALFA41*BETA13 - ALFA42*BETA23)/BETA33
00268 !
00269            BETA14 =  A14
00270            BETA24 =  A24 - ALFA21*BETA14
00271            BETA34 =  A34 - ALFA31*BETA14 - ALFA32*BETA24
00272            BETA44 =  A44 - ALFA41*BETA14 - ALFA42*BETA24 - ALFA43*BETA34
00273 !
00274 ! STORES IN XB AND W2,...,W4
00275 ! L D U FACTORISATION AT THE SAME TIME
00276 !
00277            XB(IELEM,1 ) = BETA12
00278            XB(IELEM,2 ) = BETA13
00279            XB(IELEM,3 ) = BETA14
00280            XB(IELEM,4 ) = BETA23/BETA22
00281            XB(IELEM,5 ) = BETA24/BETA22
00282            XB(IELEM,6 ) = BETA34/BETA33
00283 !
00284            XB(IELEM,07) = ALFA21
00285            XB(IELEM,08) = ALFA31
00286            XB(IELEM,09) = ALFA41
00287            XB(IELEM,10) = ALFA32
00288            XB(IELEM,11) = ALFA42
00289            XB(IELEM,12) = ALFA43
00290 !
00291            W(IELEM,2)    = BETA22
00292            W(IELEM,3)    = BETA33
00293            W(IELEM,4)    = BETA44
00294 !
00295         ENDDO ! IELEM
00296 !
00297 !-----------------------------------------------------------------------
00298 !
00299       ELSE
00300         IF (LNG.EQ.1) WRITE(LU,2000) TYPEXA(1:1)
00301         IF (LNG.EQ.2) WRITE(LU,2001) TYPEXA(1:1)
00302 2000    FORMAT(1X,'DLDU21 (BIEF) : TYPE DE MATRICE NON PREVU :',A1)
00303 2001    FORMAT(1X,'DLDU21 (BIEF) : TYPE OF MATRIX NOT AVAILABLE :',A1)
00304         CALL PLANTE(0)
00305         STOP
00306       ENDIF
00307 !
00308 !-----------------------------------------------------------------------
00309 !
00310 !  MULTIPLICATIVE ASSEMBLY OF THE DIAGONAL WITH INITIALISATION OF DB TO 1
00311 !  SKIPS IKLE1 BECAUSE W1 = 1
00312 !
00313       CALL ASMVEC(DB,IKLE(1,2),NPOIN,NELEM,NELMAX,3,W(1,2),.TRUE.,LV)
00314 !
00315 !  INVERTS DB
00316 !
00317       CALL OV( 'X=1/Y   ' , DB , DB , Z , C , NPOIN )
00318 !
00319 !-----------------------------------------------------------------------
00320 !
00321       RETURN
00322       END

Generated on Fri Aug 31 2013 18:12:58 by S.E.Bourban (HRW) using doxygen 1.7.0