dldu41.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\dldu41.f
00002 !
00089                      SUBROUTINE DLDU41
00090 !                    *****************
00091 !
00092      &(DB,XB,TYPDIA,XA,TYPEXA,
00093      & IKLE,NELEM,NELMAX,NPOIN,W,COPY,LV)
00094 !
00095 !***********************************************************************
00096 ! BIEF   V6P1                                   21/08/2010
00097 !***********************************************************************
00098 !
00099 !
00100 !
00101 !
00102 !
00103 !
00104 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00105 !| COPY           |-->| IF .TRUE. A IS COPIED INTO B.
00106 !|                |   | IF .FALSE. B IS CONSIDERED ALREADY INITIALISED
00107 !| DB             |<--| DIAGONAL OF MATRIX B
00108 !| IKLE           |-->| CONNECTIVITY TABLE
00109 !| LV             |-->| VECTOR LENGTH OF THE COMPUTER
00110 !| NELEM          |-->| NUMBER OF ELEMENTS
00111 !| NELMAX         |-->| MAXIMUM NUMBER OF ELEMENTS
00112 !| NPOIN          |-->| NUMBER OF POINTS
00113 !| TYPDIA         |<--| TYPE OF DIAGONAL ( 'Q', 'I' , OR '0' )
00114 !| TYPEXA         |<--| TYPE OF OFF-DIAGONAL TERMS ('Q','S',OR '0')
00115 !| W              |-->| WORK ARRAY OF DIMENSION (NELMAX,3)
00116 !| XA             |<--| OFF-DIAGONAL TERMS OF MATRIX A
00117 !| XB             |<--| OFF-DIAGONAL TERMS OF MATRIX B
00118 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00119 !
00120       USE BIEF, EX_DLDU41 => DLDU41
00121 !
00122       IMPLICIT NONE
00123       INTEGER LNG,LU
00124       COMMON/INFO/LNG,LU
00125 !
00126 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00127 !
00128       INTEGER, INTENT(IN)           :: NELEM,NELMAX,LV,NPOIN
00129       DOUBLE PRECISION, INTENT(OUT) :: DB(NPOIN),XB(NELMAX,*)
00130       DOUBLE PRECISION, INTENT(IN)  :: XA(NELMAX,*)
00131       CHARACTER(LEN=1), INTENT(IN)  :: TYPDIA,TYPEXA
00132       INTEGER, INTENT(IN)           :: IKLE(NELMAX,*)
00133       DOUBLE PRECISION, INTENT(OUT) :: W(NELMAX,6)
00134       LOGICAL, INTENT(IN)           :: COPY
00135 !
00136 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00137 !
00138       INTEGER IELEM
00139 !
00140       DOUBLE PRECISION Z(1),C
00141 !
00142       DOUBLE PRECISION     A12,A13,A14,A15,A16
00143       DOUBLE PRECISION A21,A22,A23,A24,A25,A26
00144       DOUBLE PRECISION A31,A32,A33,A34,A35,A36
00145       DOUBLE PRECISION A41,A42,A43,A44,A45,A46
00146       DOUBLE PRECISION A51,A52,A53,A54,A55,A56
00147       DOUBLE PRECISION A61,A62,A63,A64,A65,A66
00148       DOUBLE PRECISION        BETA12,BETA13,BETA14,BETA15,BETA16
00149       DOUBLE PRECISION ALFA21,BETA22,BETA23,BETA24,BETA25,BETA26
00150       DOUBLE PRECISION ALFA31,ALFA32,BETA33,BETA34,BETA35,BETA36
00151       DOUBLE PRECISION ALFA41,ALFA42,ALFA43,BETA44,BETA45,BETA46
00152       DOUBLE PRECISION ALFA51,ALFA52,ALFA53,ALFA54,BETA55,BETA56
00153       DOUBLE PRECISION ALFA61,ALFA62,ALFA63,ALFA64,ALFA65,BETA66
00154 !
00155 !-----------------------------------------------------------------------
00156 !
00157 ! REQUIRES THAT THE DIAGONAL OF A BE THE IDENTITY
00158 !
00159       IF(TYPDIA(1:1).NE.'I'.AND.NCSIZE.LE.1) THEN
00160         IF (LNG.EQ.1) WRITE(LU,1000) TYPDIA(1:1)
00161         IF (LNG.EQ.2) WRITE(LU,1001) TYPDIA(1:1)
00162 1000    FORMAT(1X,'DLDU41 (BIEF) : DIAGONALE DE A NON EGALE A I :',A1)
00163 1001    FORMAT(1X,'DLDU41 (BIEF) : DIAGONAL OF A NOT EQUAL TO I :',A1)
00164         CALL PLANTE(0)
00165         STOP
00166       ENDIF
00167 !
00168 !-----------------------------------------------------------------------
00169 !
00170       IF(TYPEXA(1:1).EQ.'S') THEN
00171 !
00172         IF(COPY) CALL OV('X=Y     ' , XB , XA , Z , C  , NELMAX*15 )
00173 !
00174         DO IELEM=1,NELEM
00175 !
00176 ! MATRIX TO FACTORISE (SYMMETRICAL WITH 1S ON THE DIAGONAL)
00177 !
00178 ! LINE 1
00179 !          A11 = 1.D0
00180            A12 = XA(IELEM,1 )
00181            A13 = XA(IELEM,2 )
00182            A14 = XA(IELEM,3 )
00183            A15 = XA(IELEM,4 )
00184            A16 = XA(IELEM,5 )
00185 ! LINE 2
00186            A22 = 1.D0
00187            A23 = XA(IELEM,6 )
00188            A24 = XA(IELEM,7 )
00189            A25 = XA(IELEM,8 )
00190            A26 = XA(IELEM,9 )
00191 ! LINE 3
00192            A33 = 1.D0
00193            A34 = XA(IELEM,10)
00194            A35 = XA(IELEM,11)
00195            A36 = XA(IELEM,12)
00196 ! LINE 4
00197            A44 = 1.D0
00198            A45 = XA(IELEM,13)
00199            A46 = XA(IELEM,14)
00200 ! LINE 5
00201            A55 = 1.D0
00202            A56 = XA(IELEM,15)
00203 ! LINE 6
00204            A66 = 1.D0
00205 !
00206 ! CROUT L*U FACTORISATION
00207 !
00208 ! COLUMN 1
00209            ALFA21 = A12
00210            ALFA31 = A13
00211            ALFA41 = A14
00212            ALFA51 = A15
00213            ALFA61 = A16
00214 !
00215 ! COLUMN 2
00216            BETA12 =  A12
00217            BETA22 =  A22 - ALFA21*BETA12
00218            ALFA32 = (A23 - ALFA31*BETA12)/BETA22
00219            ALFA42 = (A24 - ALFA41*BETA12)/BETA22
00220            ALFA52 = (A25 - ALFA51*BETA12)/BETA22
00221            ALFA62 = (A26 - ALFA61*BETA12)/BETA22
00222 !
00223 ! COLUMN 3
00224            BETA13 =  A13
00225            BETA23 =  A23 - ALFA21*BETA13
00226            BETA33 =  A33 - ALFA31*BETA13 - ALFA32*BETA23
00227            ALFA43 = (A34 - ALFA41*BETA13 - ALFA42*BETA23)/BETA33
00228            ALFA53 = (A35 - ALFA51*BETA13 - ALFA52*BETA23)/BETA33
00229            ALFA63 = (A36 - ALFA61*BETA13 - ALFA62*BETA23)/BETA33
00230 !
00231 ! COLUMN 4
00232            BETA14 =  A14
00233            BETA24 =  A24 - ALFA21*BETA14
00234            BETA34 =  A34 - ALFA31*BETA14 - ALFA32*BETA24
00235            BETA44 =  A44 - ALFA41*BETA14 - ALFA42*BETA24 - ALFA43*BETA34
00236            ALFA54 = (A45 - ALFA51*BETA14 - ALFA52*BETA24 - ALFA53*BETA34
00237      &     )/BETA44
00238            ALFA64 = (A46 - ALFA61*BETA14 - ALFA62*BETA24 - ALFA63*BETA34
00239      &     )/BETA44
00240 !
00241 ! COLUMN 5
00242            BETA15 =  A15
00243            BETA25 =  A25 - ALFA21*BETA15
00244            BETA35 =  A35 - ALFA31*BETA15 - ALFA32*BETA25
00245            BETA45 =  A45 - ALFA41*BETA15 - ALFA42*BETA25 - ALFA43*BETA35
00246            BETA55 =  A55 - ALFA51*BETA15 - ALFA52*BETA25 - ALFA53*BETA35
00247      &                   - ALFA54*BETA45
00248            ALFA65 = (A56 - ALFA61*BETA15 - ALFA62*BETA25 - ALFA63*BETA35
00249      &                   - ALFA64*BETA45
00250      &     )/BETA55
00251 !
00252 ! COLUMN 6
00253            BETA16 =  A16
00254            BETA26 =  A26 - ALFA21*BETA16
00255            BETA36 =  A36 - ALFA31*BETA16 - ALFA32*BETA26
00256            BETA46 =  A46 - ALFA41*BETA16 - ALFA42*BETA26 - ALFA43*BETA36
00257            BETA56 =  A56 - ALFA51*BETA16 - ALFA52*BETA26 - ALFA53*BETA36
00258      &                   - ALFA54*BETA46
00259            BETA66 =  A66 - ALFA61*BETA16 - ALFA62*BETA26 - ALFA63*BETA36
00260      &                   - ALFA64*BETA46 - ALFA65*BETA56
00261 !
00262 ! STORES IN XB AND W2,...,W6
00263 !
00264            XB(IELEM,1 ) = ALFA21
00265            XB(IELEM,2 ) = ALFA31
00266            XB(IELEM,3 ) = ALFA41
00267            XB(IELEM,4 ) = ALFA51
00268            XB(IELEM,5 ) = ALFA61
00269 !
00270            XB(IELEM,6 ) = ALFA32
00271            XB(IELEM,7 ) = ALFA42
00272            XB(IELEM,8 ) = ALFA52
00273            XB(IELEM,9 ) = ALFA62
00274 !
00275            XB(IELEM,10) = ALFA43
00276            XB(IELEM,11) = ALFA53
00277            XB(IELEM,12) = ALFA63
00278 !
00279            XB(IELEM,13) = ALFA54
00280            XB(IELEM,14) = ALFA64
00281 !
00282            XB(IELEM,15) = ALFA65
00283 !
00284            W(IELEM,2)    = BETA22
00285            W(IELEM,3)    = BETA33
00286            W(IELEM,4)    = BETA44
00287            W(IELEM,5)    = BETA55
00288            W(IELEM,6)    = BETA66
00289 !
00290         ENDDO ! IELEM
00291 !
00292 !-----------------------------------------------------------------------
00293 !
00294       ELSEIF(TYPEXA(1:1).EQ.'Q') THEN
00295 !
00296         IF(COPY) CALL OV('X=Y     ' , XB , XA , Z , C  , NELMAX*30 )
00297 !
00298         DO IELEM=1,NELEM
00299 !
00300 ! MATRIX TO FACTORISE (WITH 1S ON THE DIAGONAL)
00301 !
00302 ! LINE 1
00303 !          A11 = 1.D0
00304            A12 = XA(IELEM,1 )
00305            A13 = XA(IELEM,2 )
00306            A14 = XA(IELEM,3 )
00307            A15 = XA(IELEM,4 )
00308            A16 = XA(IELEM,5 )
00309 ! LINE 2
00310            A21 = XA(IELEM,16)
00311            A22 = 1.D0
00312            A23 = XA(IELEM,6 )
00313            A24 = XA(IELEM,7 )
00314            A25 = XA(IELEM,8 )
00315            A26 = XA(IELEM,9 )
00316 ! LINE 3
00317            A31 = XA(IELEM,17)
00318            A32 = XA(IELEM,21)
00319            A33 = 1.D0
00320            A34 = XA(IELEM,10)
00321            A35 = XA(IELEM,11)
00322            A36 = XA(IELEM,12)
00323 ! LINE 4
00324            A41 = XA(IELEM,18)
00325            A42 = XA(IELEM,22)
00326            A43 = XA(IELEM,25)
00327            A44 = 1.D0
00328            A45 = XA(IELEM,13)
00329            A46 = XA(IELEM,14)
00330 ! LINE 5
00331            A51 = XA(IELEM,19)
00332            A52 = XA(IELEM,23)
00333            A53 = XA(IELEM,26)
00334            A54 = XA(IELEM,28)
00335            A55 = 1.D0
00336            A56 = XA(IELEM,15)
00337 ! LINE 6
00338            A61 = XA(IELEM,20)
00339            A62 = XA(IELEM,24)
00340            A63 = XA(IELEM,27)
00341            A64 = XA(IELEM,29)
00342            A65 = XA(IELEM,30)
00343            A66 = 1.D0
00344 !
00345 ! CROUT L*U FACTORISATION
00346 !
00347 ! COLUMN 1
00348            ALFA21 = A21
00349            ALFA31 = A31
00350            ALFA41 = A41
00351            ALFA51 = A51
00352            ALFA61 = A61
00353 !
00354 ! COLUMN 2
00355            BETA12 =  A12
00356            BETA22 =  A22 - ALFA21*BETA12
00357            ALFA32 = (A32 - ALFA31*BETA12)/BETA22
00358            ALFA42 = (A42 - ALFA41*BETA12)/BETA22
00359            ALFA52 = (A52 - ALFA51*BETA12)/BETA22
00360            ALFA62 = (A62 - ALFA61*BETA12)/BETA22
00361 !
00362 ! COLUMN 3
00363            BETA13 =  A13
00364            BETA23 =  A23 - ALFA21*BETA13
00365            BETA33 =  A33 - ALFA31*BETA13 - ALFA32*BETA23
00366            ALFA43 = (A43 - ALFA41*BETA13 - ALFA42*BETA23)/BETA33
00367            ALFA53 = (A53 - ALFA51*BETA13 - ALFA52*BETA23)/BETA33
00368            ALFA63 = (A63 - ALFA61*BETA13 - ALFA62*BETA23)/BETA33
00369 !
00370 ! COLUMN 4
00371            BETA14 =  A14
00372            BETA24 =  A24 - ALFA21*BETA14
00373            BETA34 =  A34 - ALFA31*BETA14 - ALFA32*BETA24
00374            BETA44 =  A44 - ALFA41*BETA14 - ALFA42*BETA24 - ALFA43*BETA34
00375            ALFA54 = (A54 - ALFA51*BETA14 - ALFA52*BETA24 - ALFA53*BETA34
00376      &     )/BETA44
00377            ALFA64 = (A64 - ALFA61*BETA14 - ALFA62*BETA24 - ALFA63*BETA34
00378      &     )/BETA44
00379 !
00380 ! COLUMN 5
00381            BETA15 =  A15
00382            BETA25 =  A25 - ALFA21*BETA15
00383            BETA35 =  A35 - ALFA31*BETA15 - ALFA32*BETA25
00384            BETA45 =  A45 - ALFA41*BETA15 - ALFA42*BETA25 - ALFA43*BETA35
00385            BETA55 =  A55 - ALFA51*BETA15 - ALFA52*BETA25 - ALFA53*BETA35
00386      &                   - ALFA54*BETA45
00387            ALFA65 = (A65 - ALFA61*BETA15 - ALFA62*BETA25 - ALFA63*BETA35
00388      &                   - ALFA64*BETA45
00389      &     )/BETA55
00390 !
00391 ! COLUMN 6
00392            BETA16 =  A16
00393            BETA26 =  A26 - ALFA21*BETA16
00394            BETA36 =  A36 - ALFA31*BETA16 - ALFA32*BETA26
00395            BETA46 =  A46 - ALFA41*BETA16 - ALFA42*BETA26 - ALFA43*BETA36
00396            BETA56 =  A56 - ALFA51*BETA16 - ALFA52*BETA26 - ALFA53*BETA36
00397      &                   - ALFA54*BETA46
00398            BETA66 =  A66 - ALFA61*BETA16 - ALFA62*BETA26 - ALFA63*BETA36
00399      &                   - ALFA64*BETA46 - ALFA65*BETA56
00400 !
00401 ! STORES IN XB AND W2,...,W6
00402 ! L D U FACTORISATION AT THE SAME TIME
00403 !
00404            XB(IELEM,1 ) = BETA12
00405            XB(IELEM,2 ) = BETA13
00406            XB(IELEM,3 ) = BETA14
00407            XB(IELEM,4 ) = BETA15
00408            XB(IELEM,5 ) = BETA16
00409 !
00410            XB(IELEM,6 ) = BETA23/BETA22
00411            XB(IELEM,7 ) = BETA24/BETA22
00412            XB(IELEM,8 ) = BETA25/BETA22
00413            XB(IELEM,9 ) = BETA26/BETA22
00414 !
00415            XB(IELEM,10) = BETA34/BETA33
00416            XB(IELEM,11) = BETA35/BETA33
00417            XB(IELEM,12) = BETA36/BETA33
00418 !
00419            XB(IELEM,13) = BETA45/BETA44
00420            XB(IELEM,14) = BETA46/BETA44
00421 !
00422            XB(IELEM,15) = BETA56/BETA55
00423 !
00424            XB(IELEM,16) = ALFA21
00425            XB(IELEM,17) = ALFA31
00426            XB(IELEM,18) = ALFA41
00427            XB(IELEM,19) = ALFA51
00428            XB(IELEM,20) = ALFA61
00429 !
00430            XB(IELEM,21) = ALFA32
00431            XB(IELEM,22) = ALFA42
00432            XB(IELEM,23) = ALFA52
00433            XB(IELEM,24) = ALFA62
00434 !
00435            XB(IELEM,25) = ALFA43
00436            XB(IELEM,26) = ALFA53
00437            XB(IELEM,27) = ALFA63
00438 !
00439            XB(IELEM,28) = ALFA54
00440            XB(IELEM,29) = ALFA64
00441 !
00442            XB(IELEM,30) = ALFA65
00443 !
00444            W(IELEM,2)    = BETA22
00445            W(IELEM,3)    = BETA33
00446            W(IELEM,4)    = BETA44
00447            W(IELEM,5)    = BETA55
00448            W(IELEM,6)    = BETA66
00449 !
00450         ENDDO ! IELEM
00451 !
00452 !-----------------------------------------------------------------------
00453 !
00454       ELSE
00455         IF (LNG.EQ.1) WRITE(LU,2000) TYPEXA(1:1)
00456         IF (LNG.EQ.2) WRITE(LU,2001) TYPEXA(1:1)
00457 2000    FORMAT(1X,'DLDU41 (BIEF) : TYPE DE MATRICE NON PREVU :',A1)
00458 2001    FORMAT(1X,'DLDU41 (BIEF) : TYPE OF MATRIX NOT AVAILABLE :',A1)
00459         CALL PLANTE(0)
00460         STOP
00461       ENDIF
00462 !
00463 !-----------------------------------------------------------------------
00464 !
00465 !  MULTIPLICATIVE ASSEMBLY OF THE DIAGONAL WITH INITIALISATION OF DB TO 1
00466 !  SKIPS IKLE1 BECAUSE W1 = 1
00467 !
00468       CALL ASMVEC(DB,IKLE(1,2),NPOIN,NELEM,NELMAX,5,W(1,2),.TRUE.,LV)
00469 !
00470 !  INVERTS DB
00471 !
00472       CALL OV( 'X=1/Y   ' , DB , DB , Z , C , NPOIN )
00473 !
00474 !-----------------------------------------------------------------------
00475 !
00476       RETURN
00477       END

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