matvec.f

Go to the documentation of this file.
00001 C:\opentelemac\v7p0\sources\utils\bief\matvec.f
00002 !
00102                      SUBROUTINE MATVEC
00103 !                    *****************
00104 !
00105      &( OP , X , A , Y , C , MESH , LEGO )
00106 !
00107 !***********************************************************************
00108 ! BIEF   V6P1                                   21/08/2010
00109 !***********************************************************************
00110 !
00111 !
00112 !
00113 !
00114 !
00115 !
00116 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00117 !| A              |-->| MATRIX
00118 !| C              |-->| A GIVEN CONSTANT
00119 !| LEGO           |-->| IF PRESENT AND FALSE, NO ASSEMBLY
00120 !| MESH           |-->| MESH STRUCTURE
00121 !| OP             |-->| OPERATION TO BE DONE
00122 !| X              |<--| RESULTING VECTOR
00123 !| Y              |-->| A GIVEN VECTOR USED IN OPERATION OP
00124 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
00125 !
00126       USE BIEF, EX_MATVEC => MATVEC
00127 !
00128       IMPLICIT NONE
00129       INTEGER LNG,LU
00130       COMMON/INFO/LNG,LU
00131 !
00132 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00133 !
00134       CHARACTER*8     , INTENT(IN)           :: OP
00135       TYPE(BIEF_OBJ)  , INTENT(INOUT)        :: X
00136       TYPE(BIEF_OBJ)  , INTENT(IN)           :: A,Y
00137       DOUBLE PRECISION, INTENT(IN)           :: C
00138       TYPE(BIEF_MESH) , INTENT(INOUT)        :: MESH
00139       LOGICAL         , INTENT(IN), OPTIONAL :: LEGO
00140 !
00141 !+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
00142 !
00143       INTEGER IELM1,IELM2,IELMX,IELMY,NELEM,NELMAX,SIZXA
00144       INTEGER NPT,NPT1,NPT2,NPOIN,NPMAX,DIMIKM,NPTFR
00145 !
00146       INTEGER, DIMENSION(:), POINTER :: IKLE
00147 !
00148       LOGICAL TRANS,W_IS_FULL,LEGO2
00149 !
00150       DATA W_IS_FULL/.FALSE./
00151       SAVE W_IS_FULL
00152 !
00153 !-----------------------------------------------------------------------
00154 !
00155       IF(PRESENT(LEGO)) THEN
00156         LEGO2 = LEGO
00157       ELSE
00158         LEGO2 = .TRUE.
00159       ENDIF
00160 !
00161 !-----------------------------------------------------------------------
00162 !
00163       IF(Y%TYPE.NE.2.OR.A%TYPE.NE.3) THEN
00164         IF (LNG.EQ.1) WRITE(LU,50) X%NAME,X%TYPE
00165         IF (LNG.EQ.1) WRITE(LU,51) Y%NAME,Y%TYPE
00166         IF (LNG.EQ.1) WRITE(LU,52) A%NAME,A%TYPE
00167         IF (LNG.EQ.1) WRITE(LU,53)
00168         IF (LNG.EQ.2) WRITE(LU,60) X%NAME,X%TYPE
00169         IF (LNG.EQ.2) WRITE(LU,61) Y%NAME,Y%TYPE
00170         IF (LNG.EQ.2) WRITE(LU,62) A%NAME,A%TYPE
00171         IF (LNG.EQ.2) WRITE(LU,63)
00172         CALL PLANTE(1)
00173         STOP
00174       ENDIF
00175 !
00176 !-----------------------------------------------------------------------
00177 !
00178 !  OPERATION WITH THE MATRIX TRANSPOSE ?
00179 !
00180         TRANS =.FALSE.
00181         IF(OP(3:3).EQ.'T'.OR.OP(4:4).EQ.'T'.OR.OP(5:5).EQ.'T'.OR.
00182      &     OP(6:6).EQ.'T') TRANS = .TRUE.
00183 !
00184 !       EXTRACTS THE CHARACTERISTICS OF MATRIX M
00185 !
00186         IELM1 = A%ELMLIN
00187         IELM2 = A%ELMCOL
00188         NPT1  = BIEF_NBPTS(IELM1,MESH)
00189         NPT2  = BIEF_NBPTS(IELM2,MESH)
00190         NPT   = MIN(NPT1,NPT2)
00191 !
00192         IF(TRANS) THEN
00193           MESH%T%ELM = IELM2
00194           MESH%T%DIM1 = BIEF_NBPTS(IELM2,MESH)
00195           IELMX=IELM2
00196         ELSE
00197           MESH%T%ELM = IELM1
00198           MESH%T%DIM1 = BIEF_NBPTS(IELM1,MESH)
00199           IELMX=IELM1
00200         ENDIF
00201 !       TRIAL
00202         CALL CPSTVC(MESH%T,X)
00203 !       END TRIAL
00204 !
00205         IELMY = Y%ELM
00206 !
00207         IF(.NOT.TRANS.AND.IELM2.NE.IELMY) THEN
00208           IF (LNG.EQ.1) WRITE(LU,50) X%NAME,X%TYPE
00209           IF (LNG.EQ.1) WRITE(LU,51) Y%NAME,Y%TYPE
00210           IF (LNG.EQ.1) WRITE(LU,52) A%NAME,A%TYPE
00211           IF (LNG.EQ.1) WRITE(LU,54) IELM1,IELM2,IELMY
00212           IF (LNG.EQ.2) WRITE(LU,60) X%NAME,X%TYPE
00213           IF (LNG.EQ.2) WRITE(LU,61) Y%NAME,Y%TYPE
00214           IF (LNG.EQ.2) WRITE(LU,62) A%NAME,A%TYPE
00215           IF (LNG.EQ.2) WRITE(LU,64) IELM1,IELM2,IELMY
00216 50        FORMAT(1X,'MATVEC (BIEF) : NOM DE X : ',A6,'  TYPE : ',1I6)
00217 51        FORMAT(1X,'                NOM DE Y : ',A6,'  TYPE : ',1I6)
00218 52        FORMAT(1X,'                NOM DE A : ',A6,'  TYPE : ',1I6)
00219 53        FORMAT(1X,'                CAS NON PREVU')
00220 54        FORMAT(1X,'A ET Y INCOMPATIBLES  : ',1I6,2X,1I6,' ET ',1I6)
00221 60        FORMAT(1X,'MATVEC (BIEF) : NAME OF X : ',A6,'  TYPE : ',1I6)
00222 61        FORMAT(1X,'                NAME OF Y : ',A6,'  TYPE : ',1I6)
00223 62        FORMAT(1X,'                NAME OF A : ',A6,'  TYPE : ',1I6)
00224 63        FORMAT(1X,'                NOT IMPLEMENTED')
00225 64        FORMAT(1X,'A AND Y INCOMPATIBLE  : ',1I6,2X,1I6,' AND ',1I6)
00226           CALL PLANTE(1)
00227           STOP
00228         ENDIF
00229         IF(TRANS.AND.IELM1.NE.IELMY) THEN
00230           IF (LNG.EQ.1) WRITE(LU,50) X%NAME,X%TYPE
00231           IF (LNG.EQ.1) WRITE(LU,51) Y%NAME,Y%TYPE
00232           IF (LNG.EQ.1) WRITE(LU,52) A%NAME,A%TYPE
00233           IF (LNG.EQ.1) WRITE(LU,154) IELM1,IELM2,IELMY
00234           IF (LNG.EQ.2) WRITE(LU,60) X%NAME,X%TYPE
00235           IF (LNG.EQ.2) WRITE(LU,61) Y%NAME,Y%TYPE
00236           IF (LNG.EQ.2) WRITE(LU,62) A%NAME,A%TYPE
00237           IF (LNG.EQ.2) WRITE(LU,164) IELM1,IELM2,IELMY
00238 154       FORMAT(1X,'A ET Y INCOMPATIBLES : ',1I6,2X,1I6,' ET ',1I6,/,
00239      &           1X,'POUR UNE OPERATION SUR LA TRANSPOSEE DE A')
00240 164       FORMAT(1X,'A AND Y INCOMPATIBLE  : ',1I6,2X,1I6,' AND ',1I6,/,
00241      &           1X,'BECAUSE THE TRANSPOSED OF A IS USED')
00242           CALL PLANTE(1)
00243           STOP
00244         ENDIF
00245         IF(.NOT.TRANS.AND.IELM1.NE.IELMX) THEN
00246           IF (LNG.EQ.1) WRITE(LU,50) X%NAME,X%TYPE
00247           IF (LNG.EQ.1) WRITE(LU,51) Y%NAME,Y%TYPE
00248           IF (LNG.EQ.1) WRITE(LU,52) A%NAME,A%TYPE
00249           IF (LNG.EQ.1) WRITE(LU,55) IELM1,IELMX
00250           IF (LNG.EQ.2) WRITE(LU,60) X%NAME,X%TYPE
00251           IF (LNG.EQ.2) WRITE(LU,61) Y%NAME,Y%TYPE
00252           IF (LNG.EQ.2) WRITE(LU,62) A%NAME,A%TYPE
00253           IF (LNG.EQ.2) WRITE(LU,65) IELM1,IELMX
00254 55        FORMAT(1X,'A ET X INCOMPATIBLES : ',1I6,2X,1I6,' ET ',1I6)
00255 65        FORMAT(1X,'A AND X INCOMPATIBLE  : ',1I6,2X,1I6,' AND ',1I6)
00256           CALL PLANTE(1)
00257           STOP
00258         ENDIF
00259         IF(TRANS.AND.IELM2.NE.IELMX) THEN
00260           IF (LNG.EQ.1) WRITE(LU,50) X%NAME,X%TYPE
00261           IF (LNG.EQ.1) WRITE(LU,51) Y%NAME,Y%TYPE
00262           IF (LNG.EQ.1) WRITE(LU,52) A%NAME,A%TYPE
00263           IF (LNG.EQ.1) WRITE(LU,155) IELM1,IELMX
00264           IF (LNG.EQ.2) WRITE(LU,60) X%NAME,X%TYPE
00265           IF (LNG.EQ.2) WRITE(LU,61) Y%NAME,Y%TYPE
00266           IF (LNG.EQ.2) WRITE(LU,62) A%NAME,A%TYPE
00267           IF (LNG.EQ.2) WRITE(LU,165) IELM1,IELMX
00268 155       FORMAT(1X,'A ET X INCOMPATIBLES : ',1I6,2X,1I6,' ET ',1I6,/,
00269      &           1X,'POUR UNE OPERATION SUR LA TRANSPOSEE DE A')
00270 165       FORMAT(1X,'A AND X INCOMPATIBLE: ',1I6,2X,1I6,' AND ',/,1X,
00271      &           1X,'BECAUSE THE TRANSPOSED OF A IS USED')
00272           CALL PLANTE(1)
00273           STOP
00274         ENDIF
00275 !
00276         IF(DIMENS(IELM1).EQ.MESH%DIM) THEN
00277 !
00278 !         NORMAL MATRIX
00279 !
00280           NELEM = MESH%NELEM
00281           NELMAX= MESH%NELMAX
00282           IKLE=>MESH%IKLE%I
00283           IF(A%STO.EQ.1) THEN
00284             SIZXA=NELMAX
00285           ELSEIF(A%STO.EQ.3) THEN
00286             SIZXA=A%X%DIM1
00287           ELSE
00288             IF(LNG.EQ.1) THEN
00289               WRITE(LU,*) 'STOCKAGE INCONNU DANS MATVEC : ',A%STO
00290             ENDIF
00291             IF(LNG.EQ.2) THEN
00292               WRITE(LU,*) 'UNKNOWN STORAGE IN MATVEC : ',A%STO
00293             ENDIF
00294             CALL PLANTE(1)
00295             STOP
00296           ENDIF
00297 !
00298         ELSE
00299 !
00300 !         BOUNDARY MATRIX (NEVER WITH EDGE-BASED STORAGE)
00301 !
00302           NELEM = MESH%NPTFR
00303           NELMAX= MESH%NPTFRX
00304           IKLE=>MESH%NBOR%I
00305           SIZXA=NELMAX
00306 !
00307         ENDIF
00308 !
00309         NPTFR= MESH%NPTFR
00310         NPOIN= MESH%NPOIN
00311         NPMAX= MESH%NPMAX
00312         DIMIKM=MESH%IKLEM1%DIM1
00313 !
00314 !-----------------------------------------------------------------------
00315 !       CALLS MATVCT
00316 !-----------------------------------------------------------------------
00317 !
00318       IF(W_IS_FULL.AND.OP(3:3).NE.'X') THEN
00319 !
00320         IF (LNG.EQ.1) WRITE(LU,500)
00321         IF (LNG.EQ.2) WRITE(LU,501)
00322 500     FORMAT(1X,'MATVEC (BIEF) : UN APPEL AVEC LEGO = .FALSE.',/,
00323      &         1X,'                DOIT ETRE SUIVI D''UN APPEL AVEC',/,
00324      &         1X,'                OP=''X=X+....''')
00325 501     FORMAT(1X,'MATVEC (BIEF) : A CALL WITH LEGO = .FALSE.',/,
00326      &         1X,'                MUST BE FOLLOWED BY A CALL WITH',/,
00327      &         1X,'                OP=''X=X+....''')
00328         CALL PLANTE(1)
00329         STOP
00330 !
00331       ELSEIF(W_IS_FULL.OR.(X%NAME.NE.Y%NAME.AND.A%STO.EQ.3)) THEN
00332 !
00333         CALL MATVCT( OP,X%R,A%D%R,A%TYPDIA,A%X%R,A%TYPEXT,Y%R,
00334      &               C,IKLE,NPT,NELEM,NELMAX,MESH%W%R,
00335      &               LEGO2,IELM1,IELM2,IELMX,MESH%LV,A%STO,A%PRO,
00336      &               MESH%IKLEM1%I,DIMIKM,MESH%LIMVOI%I,MESH%MXPTVS,
00337      &               NPMAX,NPOIN,NPTFR,
00338      &               MESH%GLOSEG%I,MESH%GLOSEG%MAXDIM1,SIZXA,
00339      &               BIEF_NBPEL(IELMX,MESH),MESH)
00340 !
00341       ELSE
00342 !
00343         IF(TRANS) THEN
00344 !
00345           CALL MATVCT( 'X=TAY   ',MESH%T%R,A%D%R,A%TYPDIA,A%X%R,
00346      &                 A%TYPEXT,Y%R,C,IKLE,NPT,NELEM,NELMAX,MESH%W%R,
00347      &                 LEGO2,IELM1,IELM2,IELMX,MESH%LV,A%STO,A%PRO,
00348      &                 MESH%IKLEM1%I,DIMIKM,MESH%LIMVOI%I,MESH%MXPTVS,
00349      &                 NPMAX,NPOIN,NPTFR,
00350      &                 MESH%GLOSEG%I,MESH%GLOSEG%MAXDIM1,SIZXA,
00351      &                 BIEF_NBPEL(IELMX,MESH),MESH)
00352 !
00353           IF(OP(3:8).EQ.'TAY   ') THEN
00354             CALL OS( 'X=Y     ' , X , MESH%T , X , C )
00355           ELSEIF(OP(3:8).EQ.'-TAY  ') THEN
00356             CALL OS( 'X=-Y    ' , X , MESH%T , X , C )
00357           ELSEIF(OP(3:8).EQ.'X+TAY ') THEN
00358             CALL OS( 'X=X+Y   ' , X , MESH%T , X , C )
00359           ELSEIF(OP(3:8).EQ.'X-TAY ') THEN
00360             CALL OS( 'X=X-Y   ' , X , MESH%T , X , C )
00361           ELSEIF(OP(3:8).EQ.'X+CTAY') THEN
00362             CALL OS( 'X=X+CY  ' , X , MESH%T , X , C )
00363 !         ELSEIF(OP(3:8).EQ.'CTAY  ') THEN
00364 !           CALL OS( 'X=CY    ' , X , MESH%T , X , C )
00365           ELSE
00366            IF (LNG.EQ.1) WRITE(LU,3000) OP
00367            IF (LNG.EQ.2) WRITE(LU,3001) OP
00368 3000       FORMAT(1X,'MATVEC (BIEF) : OPERATION INCONNUE : ',A8)
00369 3001       FORMAT(1X,'MATVEC (BIEF) : UNKNOWN OPERATION : ',A8)
00370            CALL PLANTE(1)
00371            STOP
00372           ENDIF
00373 !
00374         ELSE
00375 !
00376           CALL MATVCT( 'X=AY    ',MESH%T%R,A%D%R,A%TYPDIA,A%X%R,
00377      &                 A%TYPEXT,Y%R,C,IKLE,NPT,NELEM,NELMAX,MESH%W%R,
00378      &                 LEGO2,IELM1,IELM2,IELMX,MESH%LV,A%STO,A%PRO,
00379      &                 MESH%IKLEM1%I,DIMIKM,MESH%LIMVOI%I,MESH%MXPTVS,
00380      &                 NPMAX,NPOIN,NPTFR,
00381      &                 MESH%GLOSEG%I,MESH%GLOSEG%MAXDIM1,SIZXA,
00382      &                 BIEF_NBPEL(IELMX,MESH),MESH)
00383 !
00384           IF(OP(3:8).EQ.'AY    ') THEN
00385             CALL OS( 'X=Y     ' , X , MESH%T , X , C )
00386           ELSEIF(OP(3:8).EQ.'X+AY  ') THEN
00387             CALL OS( 'X=X+Y   ' , X , MESH%T , X , C )
00388           ELSEIF(OP(3:8).EQ.'-AY   ') THEN
00389             CALL OS( 'X=-Y    ' , X , MESH%T , X , C )
00390           ELSEIF(OP(3:8).EQ.'X-AY  ') THEN
00391             CALL OS( 'X=X-Y   ' , X , MESH%T , X , C )
00392           ELSEIF(OP(3:8).EQ.'X+CAY ') THEN
00393             CALL OS( 'X=X+CY  ' , X , MESH%T , X , C )
00394           ELSEIF(OP(3:8).EQ.'CAY   ') THEN
00395             CALL OS( 'X=CY    ' , X , MESH%T , X , C )
00396           ELSE
00397            IF (LNG.EQ.1) WRITE(LU,3000) OP
00398            IF (LNG.EQ.2) WRITE(LU,3001) OP
00399            CALL PLANTE(1)
00400            STOP
00401           ENDIF
00402 !
00403         ENDIF
00404 !
00405       ENDIF
00406 !
00407 !-----------------------------------------------------------------------
00408 !
00409 !     IF LEGO WAS FALSE, MATVEC WILL HAVE TO TAKE IT INTO ACCOUNT,
00410 !     BECAUSE A NON-ASSEMBLED VECTOR WILL BE IN MESH%W
00411 !
00412       IF(.NOT.LEGO2) THEN
00413         W_IS_FULL = .TRUE.
00414       ELSE
00415         W_IS_FULL = .FALSE.
00416       ENDIF
00417 !
00418 !-----------------------------------------------------------------------
00419 !
00420       RETURN
00421       END

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