The TELEMAC-MASCARET system  trunk
matrbl.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE matrbl
3 ! *****************
4 !
5  &( op , x , a , y , c , mesh )
6 !
7 !***********************************************************************
8 ! BIEF V6P1 21/08/2010
9 !***********************************************************************
10 !
11 !brief MATRIX VECTOR OPERATIONS.
12 !code
13 !+ OP IS A STRING OF 8 CHARACTERS, WHICH INDICATES THE OPERATION TO BE
14 !+ PERFORMED ON VECTORS X,Y AND MATRIX M.
15 !+
16 !+ THE RESULT IS VECTOR X (A NON-ASSEMBLED PART OF WHICH CAN BE IN
17 !+ ARRAY W IF LEGO = .FALSE.)
18 !+
19 !+ THESE OPERATIONS ARE DIFFERENTS DEPENDING ON THE DIAGONAL TYPE
20 !+ AND THE OFF-DIAGONAL TERMS TYPE.
21 !+
22 !+ IMPLEMENTED OPERATIONS :
23 !+
24 !+ OP = 'X=AY ' : X = AY
25 !+ OP = 'X=X+AY ' : X = X + AY
26 !+ OP = 'X=X-AY ' : X = X - AY
27 !+ OP = 'X=X+CAY ' : X = X + C AY
28 !+ OP = 'X=TAY ' : X = TA Y (TA: TRANSPOSE OF A)
29 !+ OP = 'X=X+TAY ' : X = X + TA Y
30 !+ OP = 'X=X-TAY ' : X = X - TA Y
31 !+ OP = 'X=X+CTAY' : X = X + C TA Y
32 !
33 !history J-M HERVOUET (LNH) ; F LEPEINTRE (LNH)
34 !+ 24/04/97
35 !+ V5P1
36 !+
37 !
38 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
39 !+ 13/07/2010
40 !+ V6P0
41 !+ Translation of French comments within the FORTRAN sources into
42 !+ English comments
43 !
44 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
45 !+ 21/08/2010
46 !+ V6P0
47 !+ Creation of DOXYGEN tags for automated documentation and
48 !+ cross-referencing of the FORTRAN sources
49 !
50 !history R.NHEILI (Univerte de Perpignan, DALI)
51 !+ 24/02/2016
52 !+ V7P3
53 !+ ADD MODASS=3 FOR THE INTERFACE NODE ASSEMBLY
54 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
55 !| A |-->| MATRIX OR BLOCK OF MATRICES
56 !| C |-->| A GIVEN CONSTANT
57 !| MESH |-->| MESH STRUCTURE
58 !| OP |-->| THE OPERATION TO BE DONE
59 !| X |<--| RESULTING VECTOR OR BLOCK OF VECTORS
60 !| Y |-->| GIVEN VECTOR OR BLOCK OF VECTORS
61 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
62 !
63  USE bief, ex_matrbl => matrbl
64  USE declarations_telemac, ONLY : modass
65 !
67  IMPLICIT NONE
68 !
69 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
70 !
71  CHARACTER(LEN=8), INTENT(IN) :: OP
72  TYPE(bief_obj), INTENT(INOUT) :: X
73  TYPE(bief_obj), INTENT(IN) :: A,Y
74  DOUBLE PRECISION, INTENT(IN) :: C
75  TYPE(bief_mesh), INTENT(INOUT) :: MESH
76 !
77 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
78 !
79  INTEGER S
80 !
81 !-----------------------------------------------------------------------
82 !
83 ! CASE WHERE THE STRUCTURES ARE BLOCKS
84 !
85  IF(a%TYPE.EQ.4) THEN
86 !
87  s = x%N
88 !
89  IF(s.EQ.1) THEN
90 !
91  CALL matvec( op,x%ADR(1)%P,a%ADR(1)%P,y%ADR(1)%P,c,mesh)
92 !
93  ELSEIF(s.EQ.2) THEN
94 !
95  IF(op(1:8).EQ.'X=AY ') THEN
96  CALL matvec('X=AY ',
97  & x%ADR(1)%P,a%ADR(1)%P,y%ADR(1)%P,c,mesh,lego=.false.)
98  CALL matvec('X=X+AY ',
99  & x%ADR(1)%P,a%ADR(2)%P,y%ADR(2)%P,c,mesh,lego=.true.)
100  CALL matvec('X=AY ',
101  & x%ADR(2)%P,a%ADR(3)%P,y%ADR(1)%P,c,mesh,lego=.false.)
102  CALL matvec('X=X+AY ',
103  & x%ADR(2)%P,a%ADR(4)%P,y%ADR(2)%P,c,mesh,lego=.true.)
104  ELSEIF(op(1:8).EQ.'X=TAY ') THEN
105  CALL matvec('X=TAY ',
106  & x%ADR(1)%P,a%ADR(1)%P,y%ADR(1)%P,c,mesh,lego=.false.)
107  CALL matvec('X=X+TAY ',
108  & x%ADR(1)%P,a%ADR(3)%P,y%ADR(2)%P,c,mesh,lego=.true.)
109  CALL matvec('X=TAY ',
110  & x%ADR(2)%P,a%ADR(2)%P,y%ADR(1)%P,c,mesh,lego=.false.)
111  CALL matvec('X=X+TAY ',
112  & x%ADR(2)%P,a%ADR(4)%P,y%ADR(2)%P,c,mesh,lego=.true.)
113 !
114  ELSE
115  WRITE(lu,11) op
116  CALL plante(1)
117  stop
118  ENDIF
119 !
120  ELSEIF(s.EQ.3) THEN
121 !
122  IF(op(1:8).EQ.'X=AY ') THEN
123  CALL matvec('X=AY ',
124  & x%ADR(1)%P,a%ADR(1)%P,y%ADR(1)%P,c,mesh,lego=.false.)
125  CALL matvec('X=X+AY ',
126  & x%ADR(1)%P,a%ADR(2)%P,y%ADR(2)%P,c,mesh,lego=.false.)
127  CALL matvec('X=X+AY ',
128  & x%ADR(1)%P,a%ADR(3)%P,y%ADR(3)%P,c,mesh,lego=.true.)
129  CALL matvec('X=AY ',
130  & x%ADR(2)%P,a%ADR(4)%P,y%ADR(1)%P,c,mesh,lego=.false.)
131  CALL matvec('X=X+AY ',
132  & x%ADR(2)%P,a%ADR(5)%P,y%ADR(2)%P,c,mesh,lego=.false.)
133  CALL matvec('X=X+AY ',
134  & x%ADR(2)%P,a%ADR(6)%P,y%ADR(3)%P,c,mesh,lego=.true. )
135  CALL matvec('X=AY ',
136  & x%ADR(3)%P,a%ADR(7)%P,y%ADR(1)%P,c,mesh,lego=.false.)
137  CALL matvec('X=X+AY ',
138  & x%ADR(3)%P,a%ADR(8)%P,y%ADR(2)%P,c,mesh,lego=.false.)
139  CALL matvec('X=X+AY ',
140  & x%ADR(3)%P,a%ADR(9)%P,y%ADR(3)%P,c,mesh,lego=.true.)
141  ELSEIF(op(1:8).EQ.'X=TAY ') THEN
142  CALL matvec('X=TAY ',
143  & x%ADR(1)%P,a%ADR(1)%P,y%ADR(1)%P,c,mesh,lego=.false.)
144  CALL matvec('X=X+TAY ',
145  & x%ADR(1)%P,a%ADR(4)%P,y%ADR(2)%P,c,mesh,lego=.false.)
146  CALL matvec('X=X+TAY ',
147  & x%ADR(1)%P,a%ADR(7)%P,y%ADR(3)%P,c,mesh,lego=.true.)
148  CALL matvec('X=TAY ',
149  & x%ADR(2)%P,a%ADR(2)%P,y%ADR(1)%P,c,mesh,lego=.false.)
150  CALL matvec('X=X+TAY ',
151  & x%ADR(2)%P,a%ADR(5)%P,y%ADR(2)%P,c,mesh,lego=.false.)
152  CALL matvec('X=X+TAY ',
153  & x%ADR(2)%P,a%ADR(8)%P,y%ADR(3)%P,c,mesh,lego=.true.)
154  CALL matvec('X=TAY ',
155  & x%ADR(3)%P,a%ADR(3)%P,y%ADR(1)%P,c,mesh,lego=.false.)
156  CALL matvec('X=X+TAY ',
157  & x%ADR(3)%P,a%ADR(6)%P,y%ADR(2)%P,c,mesh,lego=.false.)
158  CALL matvec('X=X+TAY ',
159  & x%ADR(3)%P,a%ADR(9)%P,y%ADR(3)%P,c,mesh,lego=.true.)
160 !
161  ELSE
162  WRITE(lu,11) op
163 11 FORMAT(1x,'MATRBL (BIEF) : UNKNOWN OPERATION : ',a8)
164  CALL plante(0)
165  stop
166  ENDIF
167 !
168  ELSE
169 !
170  WRITE(lu,151) s
171  WRITE(lu,60) x%NAME,x%TYPE
172  WRITE(lu,61) y%NAME,y%TYPE
173  WRITE(lu,62) a%NAME,a%TYPE
174  WRITE(lu,63)
175 151 FORMAT(1x,'MATRBL (BIEF) : TOO MANY VECTORS :',1i6)
176  CALL plante(1)
177  stop
178 !
179  ENDIF
180 !
181 !-----------------------------------------------------------------------
182 !
183 ! CASE WHERE THE STRUCTURES ARE NOT BLOCKS
184 !
185  ELSEIF(a%TYPE.EQ.3.AND.x%TYPE.EQ.4.AND.y%TYPE.EQ.4) THEN
186 !
187  CALL matvec( op , x%ADR(1)%P , a , y%ADR(1)%P , c , mesh )
188 !
189 !-----------------------------------------------------------------------
190 !
191  ELSEIF(a%TYPE.EQ.3.AND.x%TYPE.EQ.2.AND.y%TYPE.EQ.2) THEN
192 !
193  CALL matvec( op , x , a , y , c , mesh )
194 !
195 !-----------------------------------------------------------------------
196 !
197 ! ERROR
198 !
199  ELSE
200 !
201  WRITE(lu,60) x%NAME,x%TYPE
202  WRITE(lu,61) y%NAME,y%TYPE
203  WRITE(lu,62) a%NAME,a%TYPE
204  WRITE(lu,63)
205 60 FORMAT(1x,'MATRBL (BIEF) : NAME OF X : ',a6,' TYPE : ',1i6)
206 61 FORMAT(1x,' NAME OF Y : ',a6,' TYPE : ',1i6)
207 62 FORMAT(1x,' NAME OF A : ',a6,' TYPE : ',1i6)
208 63 FORMAT(1x,' NOT IMPLEMENTED')
209  CALL plante(1)
210  stop
211 !
212  ENDIF
213 !
214 !-----------------------------------------------------------------------
215 !
216 ! COMPLEMENTS THE VECTOR (PARALLEL MODE)
217 !
218  IF(ncsize.GT.1) THEN
219  IF (modass .EQ. 1) THEN
220  CALL parcom(x,2,mesh)
221  ELSEIF (modass .EQ. 3) THEN
222  CALL parcom_comp(x,x%ADR(1)%P%E,2,mesh)
223  ENDIF
224  ENDIF
225  IF (modass .EQ.3)THEN
226  x%ADR(1)%P%R=x%ADR(1)%P%R+x%ADR(1)%P%E
227  x%ADR(1)%P%E=0.d0
228  ENDIF
229 !
230 !-----------------------------------------------------------------------
231 !
232  RETURN
233  END
subroutine parcom_comp(X, ERRX, ICOM, MESH)
Definition: parcom_comp.f:7
subroutine parcom(X, ICOM, MESH)
Definition: parcom.f:7
subroutine matrbl(OP, X, A, Y, C, MESH)
Definition: matrbl.f:7
subroutine matvec(OP, X, A, Y, C, MESH, LEGO)
Definition: matvec.f:7
Definition: bief.f:3