The TELEMAC-MASCARET system  trunk
dldu41.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE dldu41
3 ! *****************
4 !
5  &(db,xb,typdia,xa,typexa,
6  & ikle,nelem,nelmax,npoin,w,copy,lv)
7 !
8 !***********************************************************************
9 ! BIEF V6P1 21/08/2010
10 !***********************************************************************
11 !
12 !brief L D U FACTORISATION OF THE ELEMENTARY MATRICES
13 !+ IN MATRIX A
14 !+ FOR P1 PRISMS.
15 !+
16 !+ REQUIRES THAT THE DIAGONAL OF A BE THE IDENTITY.
17 !code
18 !+ EACH ELEMENTARY MATRIX IS DECOMPOSED IN THE FORM:
19 !+
20 !+ LE * DE * UE
21 !+
22 !+ LE : LOWER TRIANGULAR WITH 1S ON THE DIAGONAL
23 !+ DE : DIAGONAL
24 !+ UE : UPPER TRIANGULAR WITH 1S ON THE DIAGONAL
25 !+
26 !+ T
27 !+ IF THE MATRIX IS SYMMETRICAL : LE = UE
28 !+
29 !+ "DE" MATRICES ARE CONSIDERED LIKE DIAGONALS OF SIZE
30 !+ NPOIN X NPOIN, WHICH ARE FILLED WITH 1S FOR THE POINTS
31 !+ WHICH DO NOT BELONG TO THE CONSIDERED ELEMENT
32 !+
33 !+ THEN PERFORMS THE PRODUCT OF ALL THESE DIAGONALS
34 !+ YIELDING DIAGONAL DB
35 !
36 !warning FOR THE NONSYMMETRICAL MATRICES: UE LE
37 !+
38 !+ UE (THE BETAS) IS STORED IN XB (. , 1 TO 15)
39 !+
40 !+ LE (THE ALFAS) IS STORED IN XB (. , 16 TO 30)
41 !
42 !history J-M HERVOUET (LNH) ; F LEPEINTRE (LNH)
43 !+ 24/04/97
44 !+ V5P1
45 !+
46 !
47 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
48 !+ 13/07/2010
49 !+ V6P0
50 !+ Translation of French comments within the FORTRAN sources into
51 !+ English comments
52 !
53 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
54 !+ 21/08/2010
55 !+ V6P0
56 !+ Creation of DOXYGEN tags for automated documentation and
57 !+ cross-referencing of the FORTRAN sources
58 !
59 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
60 !| COPY |-->| IF .TRUE. A IS COPIED INTO B.
61 !| | | IF .FALSE. B IS CONSIDERED ALREADY INITIALISED
62 !| DB |<--| DIAGONAL OF MATRIX B
63 !| IKLE |-->| CONNECTIVITY TABLE
64 !| LV |-->| VECTOR LENGTH OF THE COMPUTER
65 !| NELEM |-->| NUMBER OF ELEMENTS
66 !| NELMAX |-->| MAXIMUM NUMBER OF ELEMENTS
67 !| NPOIN |-->| NUMBER OF POINTS
68 !| TYPDIA |<--| TYPE OF DIAGONAL ( 'Q', 'I' , OR '0' )
69 !| TYPEXA |<--| TYPE OF OFF-DIAGONAL TERMS ('Q','S',OR '0')
70 !| W |-->| WORK ARRAY OF DIMENSION (NELMAX,3)
71 !| XA |<--| OFF-DIAGONAL TERMS OF MATRIX A
72 !| XB |<--| OFF-DIAGONAL TERMS OF MATRIX B
73 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
74 !
75  USE bief, ex_dldu41 => dldu41
76 !
77  USE declarations_special
78  IMPLICIT NONE
79 !
80 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
81 !
82  INTEGER, INTENT(IN) :: NELEM,NELMAX,LV,NPOIN
83  DOUBLE PRECISION, INTENT(INOUT) :: DB(npoin),XB(nelmax,*)
84  DOUBLE PRECISION, INTENT(IN) :: XA(nelmax,*)
85  CHARACTER(LEN=1), INTENT(IN) :: TYPDIA,TYPEXA
86  INTEGER, INTENT(IN) :: IKLE(nelmax,*)
87  DOUBLE PRECISION, INTENT(OUT) :: W(nelmax,6)
88  LOGICAL, INTENT(IN) :: COPY
89 !
90 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
91 !
92  INTEGER IELEM
93 !
94  DOUBLE PRECISION A12,A13,A14,A15,A16
95  DOUBLE PRECISION A21,A22,A23,A24,A25,A26
96  DOUBLE PRECISION A31,A32,A33,A34,A35,A36
97  DOUBLE PRECISION A41,A42,A43,A44,A45,A46
98  DOUBLE PRECISION A51,A52,A53,A54,A55,A56
99  DOUBLE PRECISION A61,A62,A63,A64,A65,A66
100  DOUBLE PRECISION BETA12,BETA13,BETA14,BETA15,BETA16
101  DOUBLE PRECISION ALFA21,BETA22,BETA23,BETA24,BETA25,BETA26
102  DOUBLE PRECISION ALFA31,ALFA32,BETA33,BETA34,BETA35,BETA36
103  DOUBLE PRECISION ALFA41,ALFA42,ALFA43,BETA44,BETA45,BETA46
104  DOUBLE PRECISION ALFA51,ALFA52,ALFA53,ALFA54,BETA55,BETA56
105  DOUBLE PRECISION ALFA61,ALFA62,ALFA63,ALFA64,ALFA65,BETA66
106 !
107 !-----------------------------------------------------------------------
108 !
109 ! REQUIRES THAT THE DIAGONAL OF A BE THE IDENTITY
110 !
111  IF(typdia(1:1).NE.'I'.AND.ncsize.LE.1) THEN
112  WRITE(lu,1001) typdia(1:1)
113 1001 FORMAT(1x,'DLDU41 (BIEF) : DIAGONAL OF A NOT EQUAL TO I :',a1)
114  CALL plante(0)
115  stop
116  ENDIF
117 !
118 !-----------------------------------------------------------------------
119 !
120  IF(typexa(1:1).EQ.'S') THEN
121 !
122  IF(copy) CALL ov('X=Y ', x=xb, y=xa, dim1=nelmax*15)
123 !
124  DO ielem=1,nelem
125 !
126 ! MATRIX TO FACTORISE (SYMMETRICAL WITH 1S ON THE DIAGONAL)
127 !
128 ! LINE 1
129 ! A11 = 1.D0
130  a12 = xa(ielem,1 )
131  a13 = xa(ielem,2 )
132  a14 = xa(ielem,3 )
133  a15 = xa(ielem,4 )
134  a16 = xa(ielem,5 )
135 ! LINE 2
136  a22 = 1.d0
137  a23 = xa(ielem,6 )
138  a24 = xa(ielem,7 )
139  a25 = xa(ielem,8 )
140  a26 = xa(ielem,9 )
141 ! LINE 3
142  a33 = 1.d0
143  a34 = xa(ielem,10)
144  a35 = xa(ielem,11)
145  a36 = xa(ielem,12)
146 ! LINE 4
147  a44 = 1.d0
148  a45 = xa(ielem,13)
149  a46 = xa(ielem,14)
150 ! LINE 5
151  a55 = 1.d0
152  a56 = xa(ielem,15)
153 ! LINE 6
154  a66 = 1.d0
155 !
156 ! CROUT L*U FACTORISATION
157 !
158 ! COLUMN 1
159  alfa21 = a12
160  alfa31 = a13
161  alfa41 = a14
162  alfa51 = a15
163  alfa61 = a16
164 !
165 ! COLUMN 2
166  beta12 = a12
167  beta22 = a22 - alfa21*beta12
168  alfa32 = (a23 - alfa31*beta12)/beta22
169  alfa42 = (a24 - alfa41*beta12)/beta22
170  alfa52 = (a25 - alfa51*beta12)/beta22
171  alfa62 = (a26 - alfa61*beta12)/beta22
172 !
173 ! COLUMN 3
174  beta13 = a13
175  beta23 = a23 - alfa21*beta13
176  beta33 = a33 - alfa31*beta13 - alfa32*beta23
177  alfa43 = (a34 - alfa41*beta13 - alfa42*beta23)/beta33
178  alfa53 = (a35 - alfa51*beta13 - alfa52*beta23)/beta33
179  alfa63 = (a36 - alfa61*beta13 - alfa62*beta23)/beta33
180 !
181 ! COLUMN 4
182  beta14 = a14
183  beta24 = a24 - alfa21*beta14
184  beta34 = a34 - alfa31*beta14 - alfa32*beta24
185  beta44 = a44 - alfa41*beta14 - alfa42*beta24 - alfa43*beta34
186  alfa54 = (a45 - alfa51*beta14 - alfa52*beta24 - alfa53*beta34
187  & )/beta44
188  alfa64 = (a46 - alfa61*beta14 - alfa62*beta24 - alfa63*beta34
189  & )/beta44
190 !
191 ! COLUMN 5
192  beta15 = a15
193  beta25 = a25 - alfa21*beta15
194  beta35 = a35 - alfa31*beta15 - alfa32*beta25
195  beta45 = a45 - alfa41*beta15 - alfa42*beta25 - alfa43*beta35
196  beta55 = a55 - alfa51*beta15 - alfa52*beta25 - alfa53*beta35
197  & - alfa54*beta45
198  alfa65 = (a56 - alfa61*beta15 - alfa62*beta25 - alfa63*beta35
199  & - alfa64*beta45
200  & )/beta55
201 !
202 ! COLUMN 6
203  beta16 = a16
204  beta26 = a26 - alfa21*beta16
205  beta36 = a36 - alfa31*beta16 - alfa32*beta26
206  beta46 = a46 - alfa41*beta16 - alfa42*beta26 - alfa43*beta36
207  beta56 = a56 - alfa51*beta16 - alfa52*beta26 - alfa53*beta36
208  & - alfa54*beta46
209  beta66 = a66 - alfa61*beta16 - alfa62*beta26 - alfa63*beta36
210  & - alfa64*beta46 - alfa65*beta56
211 !
212 ! STORES IN XB AND W2,...,W6
213 !
214  xb(ielem,1 ) = alfa21
215  xb(ielem,2 ) = alfa31
216  xb(ielem,3 ) = alfa41
217  xb(ielem,4 ) = alfa51
218  xb(ielem,5 ) = alfa61
219 !
220  xb(ielem,6 ) = alfa32
221  xb(ielem,7 ) = alfa42
222  xb(ielem,8 ) = alfa52
223  xb(ielem,9 ) = alfa62
224 !
225  xb(ielem,10) = alfa43
226  xb(ielem,11) = alfa53
227  xb(ielem,12) = alfa63
228 !
229  xb(ielem,13) = alfa54
230  xb(ielem,14) = alfa64
231 !
232  xb(ielem,15) = alfa65
233 !
234  w(ielem,2) = beta22
235  w(ielem,3) = beta33
236  w(ielem,4) = beta44
237  w(ielem,5) = beta55
238  w(ielem,6) = beta66
239 !
240  ENDDO ! IELEM
241 !
242 !-----------------------------------------------------------------------
243 !
244  ELSEIF(typexa(1:1).EQ.'Q') THEN
245 !
246  IF(copy) CALL ov('X=Y ', x=xb, y=xa, dim1=nelmax*30)
247 !
248  DO ielem=1,nelem
249 !
250 ! MATRIX TO FACTORISE (WITH 1S ON THE DIAGONAL)
251 !
252 ! LINE 1
253 ! A11 = 1.D0
254  a12 = xa(ielem,1 )
255  a13 = xa(ielem,2 )
256  a14 = xa(ielem,3 )
257  a15 = xa(ielem,4 )
258  a16 = xa(ielem,5 )
259 ! LINE 2
260  a21 = xa(ielem,16)
261  a22 = 1.d0
262  a23 = xa(ielem,6 )
263  a24 = xa(ielem,7 )
264  a25 = xa(ielem,8 )
265  a26 = xa(ielem,9 )
266 ! LINE 3
267  a31 = xa(ielem,17)
268  a32 = xa(ielem,21)
269  a33 = 1.d0
270  a34 = xa(ielem,10)
271  a35 = xa(ielem,11)
272  a36 = xa(ielem,12)
273 ! LINE 4
274  a41 = xa(ielem,18)
275  a42 = xa(ielem,22)
276  a43 = xa(ielem,25)
277  a44 = 1.d0
278  a45 = xa(ielem,13)
279  a46 = xa(ielem,14)
280 ! LINE 5
281  a51 = xa(ielem,19)
282  a52 = xa(ielem,23)
283  a53 = xa(ielem,26)
284  a54 = xa(ielem,28)
285  a55 = 1.d0
286  a56 = xa(ielem,15)
287 ! LINE 6
288  a61 = xa(ielem,20)
289  a62 = xa(ielem,24)
290  a63 = xa(ielem,27)
291  a64 = xa(ielem,29)
292  a65 = xa(ielem,30)
293  a66 = 1.d0
294 !
295 ! CROUT L*U FACTORISATION
296 !
297 ! COLUMN 1
298  alfa21 = a21
299  alfa31 = a31
300  alfa41 = a41
301  alfa51 = a51
302  alfa61 = a61
303 !
304 ! COLUMN 2
305  beta12 = a12
306  beta22 = a22 - alfa21*beta12
307  alfa32 = (a32 - alfa31*beta12)/beta22
308  alfa42 = (a42 - alfa41*beta12)/beta22
309  alfa52 = (a52 - alfa51*beta12)/beta22
310  alfa62 = (a62 - alfa61*beta12)/beta22
311 !
312 ! COLUMN 3
313  beta13 = a13
314  beta23 = a23 - alfa21*beta13
315  beta33 = a33 - alfa31*beta13 - alfa32*beta23
316  alfa43 = (a43 - alfa41*beta13 - alfa42*beta23)/beta33
317  alfa53 = (a53 - alfa51*beta13 - alfa52*beta23)/beta33
318  alfa63 = (a63 - alfa61*beta13 - alfa62*beta23)/beta33
319 !
320 ! COLUMN 4
321  beta14 = a14
322  beta24 = a24 - alfa21*beta14
323  beta34 = a34 - alfa31*beta14 - alfa32*beta24
324  beta44 = a44 - alfa41*beta14 - alfa42*beta24 - alfa43*beta34
325  alfa54 = (a54 - alfa51*beta14 - alfa52*beta24 - alfa53*beta34
326  & )/beta44
327  alfa64 = (a64 - alfa61*beta14 - alfa62*beta24 - alfa63*beta34
328  & )/beta44
329 !
330 ! COLUMN 5
331  beta15 = a15
332  beta25 = a25 - alfa21*beta15
333  beta35 = a35 - alfa31*beta15 - alfa32*beta25
334  beta45 = a45 - alfa41*beta15 - alfa42*beta25 - alfa43*beta35
335  beta55 = a55 - alfa51*beta15 - alfa52*beta25 - alfa53*beta35
336  & - alfa54*beta45
337  alfa65 = (a65 - alfa61*beta15 - alfa62*beta25 - alfa63*beta35
338  & - alfa64*beta45
339  & )/beta55
340 !
341 ! COLUMN 6
342  beta16 = a16
343  beta26 = a26 - alfa21*beta16
344  beta36 = a36 - alfa31*beta16 - alfa32*beta26
345  beta46 = a46 - alfa41*beta16 - alfa42*beta26 - alfa43*beta36
346  beta56 = a56 - alfa51*beta16 - alfa52*beta26 - alfa53*beta36
347  & - alfa54*beta46
348  beta66 = a66 - alfa61*beta16 - alfa62*beta26 - alfa63*beta36
349  & - alfa64*beta46 - alfa65*beta56
350 !
351 ! STORES IN XB AND W2,...,W6
352 ! L D U FACTORISATION AT THE SAME TIME
353 !
354  xb(ielem,1 ) = beta12
355  xb(ielem,2 ) = beta13
356  xb(ielem,3 ) = beta14
357  xb(ielem,4 ) = beta15
358  xb(ielem,5 ) = beta16
359 !
360  xb(ielem,6 ) = beta23/beta22
361  xb(ielem,7 ) = beta24/beta22
362  xb(ielem,8 ) = beta25/beta22
363  xb(ielem,9 ) = beta26/beta22
364 !
365  xb(ielem,10) = beta34/beta33
366  xb(ielem,11) = beta35/beta33
367  xb(ielem,12) = beta36/beta33
368 !
369  xb(ielem,13) = beta45/beta44
370  xb(ielem,14) = beta46/beta44
371 !
372  xb(ielem,15) = beta56/beta55
373 !
374  xb(ielem,16) = alfa21
375  xb(ielem,17) = alfa31
376  xb(ielem,18) = alfa41
377  xb(ielem,19) = alfa51
378  xb(ielem,20) = alfa61
379 !
380  xb(ielem,21) = alfa32
381  xb(ielem,22) = alfa42
382  xb(ielem,23) = alfa52
383  xb(ielem,24) = alfa62
384 !
385  xb(ielem,25) = alfa43
386  xb(ielem,26) = alfa53
387  xb(ielem,27) = alfa63
388 !
389  xb(ielem,28) = alfa54
390  xb(ielem,29) = alfa64
391 !
392  xb(ielem,30) = alfa65
393 !
394  w(ielem,2) = beta22
395  w(ielem,3) = beta33
396  w(ielem,4) = beta44
397  w(ielem,5) = beta55
398  w(ielem,6) = beta66
399 !
400  ENDDO ! IELEM
401 !
402 !-----------------------------------------------------------------------
403 !
404  ELSE
405  WRITE(lu,2001) typexa(1:1)
406 2001 FORMAT(1x,'DLDU41 (BIEF) : TYPE OF MATRIX NOT AVAILABLE :',a1)
407  CALL plante(0)
408  stop
409  ENDIF
410 !
411 !-----------------------------------------------------------------------
412 !
413 ! MULTIPLICATIVE ASSEMBLY OF THE DIAGONAL WITH INITIALISATION OF DB TO 1
414 ! SKIPS IKLE1 BECAUSE W1 = 1
415 !
416  CALL asmvec(db,ikle(1,2),npoin,nelem,nelmax,5,w(1,2),.true.,lv)
417 !
418 ! INVERTS DB
419 !
420  CALL ov('X=1/Y ', x=db, y=db, dim1=npoin)
421 !
422 !-----------------------------------------------------------------------
423 !
424  RETURN
425  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
declarations_special::lu
integer lu
Definition: declarations_special.F:45
dldu41
subroutine dldu41(DB, XB, TYPDIA, XA, TYPEXA, IKLE, NELEM, NELMAX, NPOIN, W, COPY, LV)
Definition: dldu41.f:8
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:36 for The TELEMAC-MASCARET system by doxygen 1.8.13