The TELEMAC-MASCARET system  trunk
matvec.f
Go to the documentation of this file.
1 ! *****************
2  SUBROUTINE matvec
3 ! *****************
4 !
5  &( op , x , a , y , c , mesh , lego )
6 !
7 !***********************************************************************
8 ! BIEF V7P2
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.
17 !+
18 !+ THESE OPERATIONS ARE DIFFERENTS DEPENDING ON THE DIAGONAL TYPE
19 !+ AND THE OFF-DIAGONAL TERMS TYPE.
20 !+
21 !+ IMPLEMENTED OPERATIONS :
22 !+
23 !+ OP = 'X=AY ' : X = AY
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=CAY ' : X = C AY
29 !+ OP = 'X=TAY ' : X = TA Y (TA: TRANSPOSE OF A)
30 !+ OP = 'X=-TAY ' : X = - TA Y (TA: TRANSPOSE OF A)
31 !+ OP = 'X=X+TAY ' : X = X + TA Y
32 !+ OP = 'X=X-TAY ' : X = X - TA Y
33 !+ OP = 'X=X+CTAY' : X = X + C TA Y
34 !+ OP = 'X=CTAY ' : X = C TA Y
35 !
36 !note IMPORTANT :
37 !+
38 !+ 1) X, Y AND A CAN BE STRUCTURES OF BLOCKS
39 !+
40 !+ IF X IS A SIMPLE VECTOR, CALLS MATVEC
41 !+
42 !+ IF X IS A BLOCK OF 2 VECTORS, CALLS MA4VEC
43 !+
44 !+ IF X IS A BLOCK OF 3 VECTORS, CALLS MA9VEC
45 !+
46 !+
47 !+
48 !+ 2) X AND Y CAN BE THE SAME AT THE TIME OF THE CALL; IN THIS
49 !+
50 !+ CASE, USES AN INTERMEDIATE WORKING ARRAY: MESH%T
51 !
52 !history J-M HERVOUET (LNHE)
53 !+ 28/12/05
54 !+ V5P6
55 !+ First version.
56 !
57 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
58 !+ 13/07/2010
59 !+ V6P0
60 !+ Translation of French comments within the FORTRAN sources into
61 !+ English comments
62 !
63 !history N.DURAND (HRW), S.E.BOURBAN (HRW)
64 !+ 21/08/2010
65 !+ V6P0
66 !+ Creation of DOXYGEN tags for automated documentation and
67 !+ cross-referencing of the FORTRAN sources
68 !
69 !history R.NHEILI (Univerte de Perpignan, DALI)
70 !+ 24/02/2016
71 !+ V7P3
72 !+ ADD MODASS=3
73 !
74 !history J-M HERVOUET (EDF LAB, LNHE)
75 !+ 22/03/2016
76 !+ V7P2
77 !+ Adding argument STOX to MATVCT.
78 !
79 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
80 !| A |-->| MATRIX
81 !| C |-->| A GIVEN CONSTANT
82 !| LEGO |-->| IF PRESENT AND FALSE, NO ASSEMBLY
83 !| MESH |-->| MESH STRUCTURE
84 !| OP |-->| OPERATION TO BE DONE
85 !| X |<--| RESULTING VECTOR
86 !| Y |-->| A GIVEN VECTOR USED IN OPERATION OP
87 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
88 !
89  USE bief, ex_matvec => matvec
90  USE declarations_telemac, ONLY : modass
91  USE declarations_telemac, ONLY : w_is_full
92  USE declarations_special
93 !
94  IMPLICIT NONE
95 !
96 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
97 !
98  CHARACTER(LEN=8) , INTENT(IN) :: OP
99  TYPE(bief_obj) , INTENT(INOUT) :: X
100  TYPE(bief_obj) , INTENT(IN) :: A,Y
101  DOUBLE PRECISION, INTENT(IN) :: C
102  TYPE(bief_mesh) , INTENT(INOUT) :: MESH
103  LOGICAL , INTENT(IN), OPTIONAL :: LEGO
104 !
105 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
106 !
107  INTEGER IELM1,IELM2,IELMX,IELMY,NELEM,NELMAX,SIZXA
108  INTEGER NPT,NPT1,NPT2,NPOIN,NPMAX,DIMIKM,NPTFR
109 !
110  INTEGER, DIMENSION(:), POINTER :: IKLE
111 !
112  LOGICAL TRANS,LEGO2
113 !
114 !-----------------------------------------------------------------------
115 !
116  IF(PRESENT(lego)) THEN
117  lego2 = lego
118  ELSE
119  lego2 = .true.
120  ENDIF
121 !
122 !-----------------------------------------------------------------------
123 !
124  IF(y%TYPE.NE.2.OR.a%TYPE.NE.3) THEN
125  WRITE(lu,60) x%NAME,x%TYPE
126  WRITE(lu,61) y%NAME,y%TYPE
127  WRITE(lu,62) a%NAME,a%TYPE
128  WRITE(lu,63)
129  CALL plante(1)
130  stop
131  ENDIF
132 !
133 !-----------------------------------------------------------------------
134 !
135 ! OPERATION WITH THE MATRIX TRANSPOSE ?
136 !
137  trans =.false.
138  IF(op(3:3).EQ.'T'.OR.op(4:4).EQ.'T'.OR.op(5:5).EQ.'T'.OR.
139  & op(6:6).EQ.'T') trans = .true.
140 !
141 ! EXTRACTS THE CHARACTERISTICS OF MATRIX M
142 !
143  ielm1 = a%ELMLIN
144  ielm2 = a%ELMCOL
145  npt1 = bief_nbpts(ielm1,mesh)
146  npt2 = bief_nbpts(ielm2,mesh)
147  npt = min(npt1,npt2)
148 !
149  IF(trans) THEN
150  mesh%T%ELM = ielm2
151  mesh%T%DIM1 = bief_nbpts(ielm2,mesh)
152  ielmx=ielm2
153  ELSE
154  mesh%T%ELM = ielm1
155  mesh%T%DIM1 = bief_nbpts(ielm1,mesh)
156  ielmx=ielm1
157  ENDIF
158 ! TRIAL
159  CALL cpstvc(mesh%T,x)
160 ! END TRIAL
161 !
162  ielmy = y%ELM
163 !
164  IF(.NOT.trans.AND.ielm2.NE.ielmy) THEN
165  WRITE(lu,60) x%NAME,x%TYPE
166  WRITE(lu,61) y%NAME,y%TYPE
167  WRITE(lu,62) a%NAME,a%TYPE
168  WRITE(lu,64) ielm1,ielm2,ielmy
169 60 FORMAT(1x,'MATVEC (BIEF) : NAME OF X : ',a6,' TYPE : ',1i6)
170 61 FORMAT(1x,' NAME OF Y : ',a6,' TYPE : ',1i6)
171 62 FORMAT(1x,' NAME OF A : ',a6,' TYPE : ',1i6)
172 63 FORMAT(1x,' NOT IMPLEMENTED')
173 64 FORMAT(1x,'A AND Y INCOMPATIBLE : ',1i6,2x,1i6,' AND ',1i6)
174  CALL plante(1)
175  stop
176  ENDIF
177  IF(trans.AND.ielm1.NE.ielmy) THEN
178  WRITE(lu,60) x%NAME,x%TYPE
179  WRITE(lu,61) y%NAME,y%TYPE
180  WRITE(lu,62) a%NAME,a%TYPE
181  WRITE(lu,164) ielm1,ielm2,ielmy
182 164 FORMAT(1x,'A AND Y INCOMPATIBLE : ',1i6,2x,1i6,' AND ',1i6,/,
183  & 1x,'BECAUSE THE TRANSPOSED OF A IS USED')
184  CALL plante(1)
185  stop
186  ENDIF
187  IF(.NOT.trans.AND.ielm1.NE.ielmx) THEN
188  WRITE(lu,60) x%NAME,x%TYPE
189  WRITE(lu,61) y%NAME,y%TYPE
190  WRITE(lu,62) a%NAME,a%TYPE
191  WRITE(lu,65) ielm1,ielmx
192 65 FORMAT(1x,'A AND X INCOMPATIBLE : ',1i6,2x,1i6,' AND ',1i6)
193  CALL plante(1)
194  stop
195  ENDIF
196  IF(trans.AND.ielm2.NE.ielmx) THEN
197  WRITE(lu,60) x%NAME,x%TYPE
198  WRITE(lu,61) y%NAME,y%TYPE
199  WRITE(lu,62) a%NAME,a%TYPE
200  WRITE(lu,165) ielm1,ielmx
201 165 FORMAT(1x,'A AND X INCOMPATIBLE: ',1i6,2x,1i6,' AND ',/,1x,
202  & 1x,'BECAUSE THE TRANSPOSED OF A IS USED')
203  CALL plante(1)
204  stop
205  ENDIF
206 !
207  IF(dimens(ielm1).EQ.mesh%DIM1) THEN
208 !
209 ! NORMAL MATRIX
210 !
211  nelem = mesh%NELEM
212  nelmax= mesh%NELMAX
213  ikle=>mesh%IKLE%I
214  IF(a%STO.EQ.1) THEN
215  sizxa=nelmax
216  ELSEIF(a%STO.EQ.3) THEN
217  sizxa=a%X%DIM1
218  ELSE
219  WRITE(lu,*) 'UNKNOWN STORAGE IN MATVEC : ',a%STO
220  CALL plante(1)
221  stop
222  ENDIF
223 !
224  ELSE
225 !
226 ! BOUNDARY MATRIX (NEVER WITH EDGE-BASED STORAGE)
227 !
228  nelem = mesh%NPTFR
229  nelmax= mesh%NPTFRX
230  ikle=>mesh%NBOR%I
231  sizxa=nelmax
232 !
233  ENDIF
234 !
235  nptfr= mesh%NPTFR
236  npoin= mesh%NPOIN
237  npmax= mesh%NPMAX
238  dimikm=mesh%IKLEM1%DIM1
239 !
240 !-----------------------------------------------------------------------
241 ! CALLS MATVCT
242 !-----------------------------------------------------------------------
243 !
244  IF(w_is_full.AND.op(3:3).NE.'X') THEN
245 !
246  WRITE(lu,501)
247 501 FORMAT(1x,'MATVEC (BIEF) : A CALL WITH LEGO = .FALSE.',/,
248  & 1x,' MUST BE FOLLOWED BY A CALL WITH',/,
249  & 1x,' OP=''X=X+....''')
250  CALL plante(1)
251  stop
252 !
253  ELSEIF(w_is_full.OR.(x%NAME.NE.y%NAME.AND.a%STO.EQ.3)) THEN
254 !
255  IF (modass .EQ.1) THEN
256  CALL matvct( op,x%R,a%D%R,a%TYPDIA,a%X%R,a%TYPEXT,y%R,
257  & c,ikle,npt,nelem,nelmax,mesh%W%R,
258  & lego2,ielm1,ielm2,ielmx,mesh%LV,a%STO,a%PRO,
259  & mesh%IKLEM1%I,dimikm,mesh%LIMVOI%I,mesh%MXPTVS,
260  & npmax,npoin,
261  & mesh%GLOSEG%I,mesh%GLOSEG%MAXDIM1,sizxa,
262  & bief_nbpel(ielmx,mesh),mesh,a%STOX)
263  ELSEIF (modass .EQ.3) THEN
264  CALL matvct( op,x%R,a%D%R,a%TYPDIA,a%X%R,a%TYPEXT,y%R,
265  & c,ikle,npt,nelem,nelmax,mesh%W%R,
266  & lego2,ielm1,ielm2,ielmx,mesh%LV,a%STO,a%PRO,
267  & mesh%IKLEM1%I,dimikm,mesh%LIMVOI%I,mesh%MXPTVS,
268  & npmax,npoin,
269  & mesh%GLOSEG%I,mesh%GLOSEG%MAXDIM1,sizxa,
270  & bief_nbpel(ielmx,mesh),mesh,a%STOX,
271  & x_err=x%E, y_err=y%E , da_err=a%D%E)
272  ENDIF
273 !
274  ELSE
275 !
276  IF(trans) THEN
277 !
278  CALL matvct( 'X=TAY ',mesh%T%R,a%D%R,a%TYPDIA,a%X%R,
279  & a%TYPEXT,y%R,c,ikle,npt,nelem,nelmax,mesh%W%R,
280  & lego2,ielm1,ielm2,ielmx,mesh%LV,a%STO,a%PRO,
281  & mesh%IKLEM1%I,dimikm,mesh%LIMVOI%I,mesh%MXPTVS,
282  & npmax,npoin,
283  & mesh%GLOSEG%I,mesh%GLOSEG%MAXDIM1,sizxa,
284  & bief_nbpel(ielmx,mesh),mesh,a%STOX)
285 !
286  IF(op(3:8).EQ.'TAY ') THEN
287  CALL os('X=Y ', x=x, y=mesh%T)
288  ELSEIF(op(3:8).EQ.'-TAY ') THEN
289  CALL os('X=-Y ', x=x, y=mesh%T)
290  ELSEIF(op(3:8).EQ.'X+TAY ') THEN
291  CALL os('X=X+Y ', x=x, y=mesh%T)
292  ELSEIF(op(3:8).EQ.'X-TAY ') THEN
293  CALL os('X=X-Y ', x=x, y=mesh%T)
294  ELSEIF(op(3:8).EQ.'X+CTAY') THEN
295  CALL os('X=X+CY ', x=x, y=mesh%T, c=c)
296  ELSE
297  WRITE(lu,3001) op
298 3001 FORMAT(1x,'MATVEC (BIEF) : UNKNOWN OPERATION : ',a8)
299  CALL plante(1)
300  stop
301  ENDIF
302 !
303  ELSE
304 !
305  IF (modass .EQ.1) THEN
306  CALL matvct( 'X=AY ',mesh%T%R,a%D%R,a%TYPDIA,a%X%R,
307  & a%TYPEXT,y%R,c,ikle,npt,nelem,nelmax,mesh%W%R,
308  & lego2,ielm1,ielm2,ielmx,mesh%LV,a%STO,a%PRO,
309  & mesh%IKLEM1%I,dimikm,mesh%LIMVOI%I,mesh%MXPTVS,
310  & npmax,npoin,
311  & mesh%GLOSEG%I,mesh%GLOSEG%MAXDIM1,sizxa,
312  & bief_nbpel(ielmx,mesh),mesh,a%STOX)
313  ELSEIF (modass .EQ.3) THEN
314  mesh%T%E=0.d0
315  CALL matvct( 'X=AY ',mesh%T%R,a%D%R,a%TYPDIA,a%X%R,
316  & a%TYPEXT,y%R,c,ikle,npt,nelem,nelmax,mesh%W%R,
317  & lego2,ielm1,ielm2,ielmx,mesh%LV,a%STO,a%PRO,
318  & mesh%IKLEM1%I,dimikm,mesh%LIMVOI%I,mesh%MXPTVS,
319  & npmax,npoin,
320  & mesh%GLOSEG%I,mesh%GLOSEG%MAXDIM1,sizxa,
321  & bief_nbpel(ielmx,mesh),mesh,a%STOX,
322  & x_err=mesh%T%E, y_err=y%E , da_err=a%D%E)
323  ENDIF
324 !
325  IF(op(3:8).EQ.'AY ') THEN
326  CALL os('X=Y ', x=x, y=mesh%T)
327  ELSEIF(op(3:8).EQ.'X+AY ') THEN
328  CALL os('X=X+Y ', x=x, y=mesh%T)
329  ELSEIF(op(3:8).EQ.'-AY ') THEN
330  CALL os('X=-Y ', x=x, y=mesh%T)
331  ELSEIF(op(3:8).EQ.'X-AY ') THEN
332  CALL os('X=X-Y ', x=x, y=mesh%T)
333  ELSEIF(op(3:8).EQ.'X+CAY ') THEN
334  CALL os('X=X+CY ', x=x, y=mesh%T, c=c)
335  ELSEIF(op(3:8).EQ.'CAY ') THEN
336  CALL os('X=CY ', x=x, y=mesh%T, c=c)
337  ELSE
338  WRITE(lu,3001) op
339  CALL plante(1)
340  stop
341  ENDIF
342 !
343  ENDIF
344 !
345  ENDIF
346 !
347 !-----------------------------------------------------------------------
348 !
349 ! IF LEGO WAS FALSE, MATVEC WILL HAVE TO TAKE IT INTO ACCOUNT,
350 ! BECAUSE A NON-ASSEMBLED VECTOR WILL BE IN MESH%W
351 !
352  IF(.NOT.lego2) THEN
353  w_is_full = .true.
354  ELSE
355  w_is_full = .false.
356  ENDIF
357 !
358 !-----------------------------------------------------------------------
359 !
360  RETURN
361  END
declarations_telemac::w_is_full
logical w_is_full
Definition: declarations_telemac.f:271
matvct
subroutine matvct(OP, X, DA, TYPDIA, XA, TYPEXT, Y, C, IKLE, NPT, NELEM, NELMAX, W, LEGO, IELM1, IELM2, IELMX, LV, S, P, IKLEM1, DIMIKM, LIMVOI, MXPTVS, NPMAX, NPOIN, GLOSEG, SIZGLO, SIZXA, NDP, MESH, STOX, X_ERR, Y_ERR, DA_ERR)
Definition: matvct.f:11
dimens
integer function dimens(IELM)
Definition: dimens.f:7
declarations_special
Definition: declarations_special.F:3
bief_nbpts
integer function bief_nbpts(IELM, MESH)
Definition: bief_nbpts.f:7
bief_nbpel
integer function bief_nbpel(IELM, MESH)
Definition: bief_nbpel.f:7
declarations_special::lu
integer lu
Definition: declarations_special.F:45
declarations_telemac
Definition: declarations_telemac.f:3
declarations_telemac::modass
integer modass
Definition: declarations_telemac.f:265
cpstvc
subroutine cpstvc(X, Y)
Definition: cpstvc.f:7
os
subroutine os(OP, X, Y, Z, C, IOPT, INFINI, ZERO)
Definition: os.f:7
matvec
subroutine matvec(OP, X, A, Y, C, MESH, LEGO)
Definition: matvec.f:7
bief
Definition: bief.f:3
Generated on Tue Jan 5 2021 14:36:44 for The TELEMAC-MASCARET system by doxygen 1.8.13