Sources Esope | Cast3M

arpver
C ARPVER    SOURCE    CB215821  20/11/25    13:18:25     10792          
      SUBROUTINE ARPVER (IPRTRA,TYPRO,I,QUAD,SYM,EPSI,INVER,
     &                                      IPVECR,IPVECI,VALP)
 
C **********************************************************************
C
C                           A R P V E R
C
C  FONCTION:
C  ---------
C
C      CALCUL DE NORMES ET DU RESIDU D'UN MODE
C
C
C  REMARQUES:
C  ---------
C
C      POUR UN MODE (x,LAMBDA) SONT CALCULES :
C
C      * SA NORME PAR RAPPORT AU B-PRODUIT SCALAIRE
C          ||X||= X*BX
C
C      * SON RESIDU
C        - EPS = Kx-LAMBDA*Mx OU Kx-LAMBDA*KSIGx  DANS LE CAS LINEAIRE
C        - EPS = Kx+LABMDA*Cx+LAMBDA**2*Mx DANS LE CAS QUADRATIQUE
C
C      * DES NORMES DU RESIDU
C        - NORME EUCLIDIENNE ||EPS|| = EPS*EPS
C          (PRODUIT SCALAIRE PAR LE CONJUGUE)
C        - NORME INFINIE ||EPS|| = MAX( ABS(RE(EPS)) + ABS(IM(EPS)) )
C
C
C  PARAMETRES:  (E)=ENTREE   (S)=SORTIE
C  -----------
C
C      IPRTRA  ENTIER     (E)   POINTEUR DES OPERATEURS DE TRAVAIL
C
C      TYPRO   ENTIER     (E)    TYPE DE PROBLEME
C
C      I       ENTIER     (E)    NUMERO DE MODE
C
C      QUAD    LOGIQUE    (E)    PROBLEME QUADRATIQUE OU NON
C
C      SYM     LOGIQUE    (E)    PROBLEME SYMETRIQUE OU NON
C
C      EPSI    REEL DP    (E)    ZERO DE TOLERANCE
C
C      INVER   LOGIQUE    (E)   .TRUE. -> PRODUIT SCALAIRE X'KX
C                               .FALSE. -> PRODUIT SCALAIRE X'MX
C
C      IPVECR  ENTIER     (E)    POINTEUR DE LA PARTIE RELLE DU MODE
C
C      IPVECI  ENTIER     (E)    POINTEUR DE LA PARTIE IMAGINAIRE
C                                DU MODE (OPTIONNEL)
C
C      VALP  COMPLEX DP (E)   VALEUR PROPRE ASSOCIEE AU MODE
C
C
C  SOUS-PROGRAMMES APPELES:
C  ------------------------
C
C    MUCPRI, MUCHPO, ADCHPO, OPCHP1, DTCHPO
C    MOTS1, MAXIM1, CORRSP, XTX1, XTMX, XTY1,
C
C
C  AUTEUR, DATE DE CREATION:
C  -------------------------
C
C      PASCAL BOUDA     18 SEPTEMBRE 2015
C
C  LANGAGE:
C  --------
C
C      FORTRAN 77 & 90
C
C **********************************************************************
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
 
-INC PPARAM
-INC CCOPTIO
c -INC TARTRAK
-INC TARWORK
 
      INTEGER IPRTRA
      INTEGER TYPRO
      INTEGER I
      LOGICAL QUAD
      LOGICAL SYM
      REAL*8 EPSI
      LOGICAL INVER
      INTEGER IPVECR
      INTEGER IPVECI
      COMPLEX*16 VALP
 
      INTEGER IPRIGI,IPMASS,IPSCAL
      INTEGER IPLMOT
      CHARACTER*(LOCOMP) MOTCLE
 
      COMPLEX*16 RAY
      COMPLEX*16 UN
      INTEGER IPLMOX, IPLMOY
      INTEGER IPMX, IPMY
      REAL*8 XMX, XMY, YMX, YMY
      REAL*8 RAYR
      LOGICAL MAL
 
      INTEGER IPCTRA(30)
      REAL*8 XTX,RXTX,IXTX,XINF
C      REAL*8 CXINF,CXTCX,ICXTCX,IKXTKX,MXINF,MXTMX,RCXTCX,RKXTKX,RMXTMX
C      REAL*8 IMXTMX,KXTKX,KXINF
      REAL*8 EUC, INF
 
 
 
C ***********************************************************************
C                                 Cas symetrique
C ***********************************************************************
 
      MAL=.FALSE.
      RAYR=0.
      RAY=CMPLX(0.,0.)
 
      MRITRA=IPRTRA
      SEGACT MRITRA
 
      IPRIGI=RIGI(1)
      IPMASS=RIGI(2)
      IPAMOR=RIGI(3)
 
 
      IF (SYM) THEN
 
C ****************
C Calcul de ||x||*
C ****************
 
        IF (.NOT. INVER) THEN
          CALL XTMX (IPVECR,IPMASS,RAYR)
        ELSE
          CALL XTMX (IPVECR,IPRIGI,RAYR)
        ENDIF
 
        RAY=CMPLX(RAYR,0.)
 
C critere de norme "admissible"
        IF (ABS(REAL(RAY)-1.) .GT. (EPSI)) THEN
          MAL=.TRUE.
        ENDIF
 
 
        IF (IIMPI .GT. 0) THEN
 
C **********************************************
C Calcul de Ax-LAMBDA*Bx et de ||Ax-LAMBDA*Bx||*
C **********************************************
 
C Initialisation des chpoints de travail
          DO j=1,30
            IPCTRA(j)=0
          ENDDO
 
C ********Calcul des produits matrices-vecteur***************************
 
          CALL MUCPRI (IPVECR,IPRIGI,IPCTRA(1))
          CALL MUCPRI (IPVECR,IPMASS,IPCTRA(2))
 
C *********Combinaisons lineaires***************************************
 
 
C IPCTRA(10) est le chpoint residu
          CALL MUCHPO(IPCTRA(2),REAL(VALP),IPCTRA(5),1)
          CALL ADCHPO(IPCTRA(1),IPCTRA(5),IPCTRA(10),1.D0,-1.D0)
 
c *Norme euclidienne de Kx
c           CALL XTX1(IPCTRA(1),KXTKX)
c
c *Norme infinie de Kx
C           IOPERA = 14
C           IARGU = 0
C           I1    = 0
C           X1    = 0.D0
C           CALL OPCHP1(IPCTRA(1),IOPERA,IARGU,I1,X1,IPCTRA(6),IRET)
c
c           CALL MOTS1 (IPLMOT,MOTCLE)
c           CALL MAXIM1 (IPCTRA(6),IPLMOT,MOTCLE,0,KXINF)
c
c
c *Norme euclidienne de lambda*Mx
c           CALL XTX1(IPCTRA(5),MXTMX)
c
c *Norme infinie de lambda*Mx
C           IOPERA = 14
C           IARGU = 0
C           I1    = 0
C           X1    = 0.D0
C           CALL OPCHP1(IPCTRA(5),IOPERA,IARGU,I1,X1,IPCTRA(7),IRET)
c
c           CALL MOTS1 (IPLMOT,MOTCLE)
c           CALL MAXIM1 (IPCTRA(7),IPLMOT,MOTCLE,0,MXINF)
 
 
C Norme euclidienne du residu
          CALL XTX1(IPCTRA(10),XTX)
 
C Norme infinie du residu
           IOPERA = 14
           IARGU = 0
           I1    = 0
           X1    = 0.
           CALL OPCHP1(IPCTRA(10),IOPERA,IARGU,I1,X1,IPCTRA(4),IRET)
 
          CALL MOTS1  (IPLMOT,MOTCLE)
          CALL MAXIM1 (IPCTRA(4),IPLMOT,MOTCLE,0,XINF)
 
C Calcul des normes
          EUC=XTX
          INF=XINF
 
        ENDIF
 
 
C ***********************************************************************
C                             Cas non symetrique
C ***********************************************************************
      ELSE
 
 
        IF (TYPRO .EQ. 3) THEN
 
C ****************
C Calcul de ||x||*
C ****************
          IF (INVER) THEN
            IPSCAL=IPRIGI
          ELSE
            IPSCAL=IPMASS
          ENDIF
 
          CALL MUCPRI (IPVECR,IPSCAL,IPMX)
          CALL MUCPRI (IPVECI,IPSCAL,IPMY)
 
C Formation d'un chpoint nul si mode reel
          IF (IPVECI .EQ. 0) THEN
            CALL MUCHPO(IPVECR,0.D0,IPVECI,1)
          ENDIF
 
          CALL CORRSP (IPSCAL,IPVECR,IPMX,IPLMOX,IPLMOY)
          CALL XTY1 (IPVECR,IPMX,IPLMOX,IPLMOY,XMX)
          CALL CORRSP (IPSCAL,IPVECR,IPMY,IPLMOX,IPLMOY)
          CALL XTY1 (IPVECR,IPMY,IPLMOX,IPLMOY,XMY)
          CALL CORRSP (IPSCAL,IPVECI,IPMX,IPLMOX,IPLMOY)
          CALL XTY1 (IPVECI,IPMX,IPLMOX,IPLMOY,YMX)
          CALL CORRSP (IPSCAL,IPVECI,IPMY,IPLMOX,IPLMOY)
          CALL XTY1 (IPVECI,IPMY,IPLMOX,IPLMOY,YMY)
 
          CALL DTCHPO(IPMX)
          CALL DTCHPO(IPMY)
 
          RAY=CMPLX(XMX+YMY,XMY-YMX)
 
C Dans le cas quadratique, le calcul se fait par blocs
 
 
          IF (QUAD) THEN
            UN=CMPLX(1.D0,0.D0)
            RAY=(UN+ABS(VALP)**2)*RAY
          ENDIF
 
C critere de norme "admissible"
          IF ( ABS(REAL(RAY)-1.) .GT. EPSI .OR.
     &                   ABS(AIMAG(RAY)) .GT. EPSI ) THEN
            MAL=.TRUE.
          ENDIF
 
        ENDIF
 
 
        IF (IIMPI .GT. 1) THEN
 
C **********************************************
C Calcul de Ax-LAMBDA*Bx et de ||Ax-LAMBDA*Bx||*
C **********************************************
 
C Initialisation des chpoints de travail
          DO j=1,30
            IPCTRA(j)=0
          ENDDO
 
C ********Calcul des produits matrices-vecteur***************************
          CALL MUCPRI (IPVECR,IPRIGI,IPCTRA(1))
          CALL MUCPRI (IPVECR,IPMASS,IPCTRA(2))
          CALL MUCPRI (IPVECI,IPRIGI,IPCTRA(3))
          CALL MUCPRI (IPVECI,IPMASS,IPCTRA(4))
 
c *Norme euclidienne de Kx
c           CALL XTX1(IPCTRA(1),RKXTKX)
c           CALL XTX1(IPCTRA(3),IKXTKX)
c           KXTKX=RKXTKX+IKXTKX
c
c *Norme infinie de Kx
C           IOPERA = 14
C           IARGU = 0
C           I1    = 0
C           X1    = 0.D0
C           CALL OPCHP1(IPCTRA(1),IOPERA,IARGU,I1,X1,IPCTRA(18),IRET)
C           CALL OPCHP1(IPCTRA(3),IOPERA,IARGU,I1,X1,IPCTRA(19),IRET)
c           CALL ADCHPO(IPCTRA(18),IPCTRA(19),IPCTRA(20),1.D0,1.D0)
c
c           CALL MOTS1 (IPLMOT,MOTCLE)
c           CALL MAXIM1 (IPCTRA(20),IPLMOT,MOTCLE,0,KXINF)
 
 
          IF (QUAD) THEN
 
            CALL MUCPRI (IPVECR,IPAMOR,IPCTRA(5))
            CALL MUCPRI (IPVECI,IPAMOR,IPCTRA(6))
 
          ENDIF
 
C *********Combinaisons lineaires***************************************
 
C IPCTRA(10) est la partie reelle du chpoint residu
C IPCTRA(14) est la partie imaginaire du chpoint residu
 
          IF (.NOT.QUAD) THEN
 
C Partie reelle
            CALL ADCHPO(IPCTRA(2),IPCTRA(4),IPCTRA(7),
     &                                         REAL(VALP),-AIMAG(VALP))
            CALL ADCHPO(IPCTRA(1),IPCTRA(7),IPCTRA(10),1.D0,-1.D0)
C Partie imaginaire
            CALL ADCHPO(IPCTRA(2),IPCTRA(4),IPCTRA(8),
     &                                           AIMAG(VALP),REAL(VALP))
            CALL ADCHPO(IPCTRA(1),IPCTRA(8),IPCTRA(14),1.D0,-1.D0)
 
c *Norme euclidienne de lambda*Mx
c             CALL XTX1(IPCTRA(7),RMXTMX)
c             CALL XTX1(IPCTRA(8),IMXTMX)
c             MXTMX=RMXTMX+IMXTMX
c
c *Norme infinie de lambda*Mx
C             IOPERA = 14
C             IARGU = 0
C             I1    = 0
C             X1    = 0.D0
C             CALL OPCHP1(IPCTRA(7),IOPERA,IARGU,I1,X1,IPCTRA(27),IRET)
C             CALL OPCHP1(IPCTRA(8),IOPERA,IARGU,I1,X1,IPCTRA(28),IRET)
c             CALL ADCHPO (IPCTRA(27),IPCTRA(28),IPCTRA(29),1.D0,1.D0)
c
c             CALL MOTS1 (IPLMOT,MOTCLE)
c             CALL MAXIM1 (IPCTRA(29),IPLMOT,MOTCLE,0,MXINF)
 
          ELSE
 
C Partie reelle
            CALL ADCHPO(IPCTRA(5),IPCTRA(6),IPCTRA(7),
     &                                         REAL(VALP),-AIMAG(VALP))
            CALL ADCHPO(IPCTRA(2),IPCTRA(4),IPCTRA(8),
     &           REAL(VALP)**2-AIMAG(VALP)**2,-2*REAL(VALP)*AIMAG(VALP))
            CALL ADCHPO(IPCTRA(1),IPCTRA(7),IPCTRA(9),1.D0,1.D0)
            CALL ADCHPO(IPCTRA(8),IPCTRA(9),IPCTRA(10),1.D0,1.D0)
 
C Partie imaginaire
            CALL ADCHPO(IPCTRA(5),IPCTRA(6),IPCTRA(11),
     &                                         REAL(VALP),AIMAG(VALP))
            CALL ADCHPO(IPCTRA(4),IPCTRA(2),IPCTRA(12),
     &           REAL(VALP)**2-AIMAG(VALP)**2,2*REAL(VALP)*AIMAG(VALP))
            CALL ADCHPO(IPCTRA(3),IPCTRA(11),IPCTRA(13),1.D0,1.D0)
            CALL ADCHPO(IPCTRA(12),IPCTRA(13),IPCTRA(14),1.D0,1.D0)
 
c *Norme euclidienne de lambda*Cx
c             CALL XTX1(IPCTRA(7),RCXTCX)
c             CALL XTX1(IPCTRA(9),ICXTCX)
c             CXTCX=RCXTCX+ICXTCX
c
c *Norme infinie de lambda*Cx
C             IOPERA = 14
C             IARGU = 0
C             I1    = 0
C             X1    = 0.D0
C             CALL OPCHP1(IPCTRA(7),IOPERA,IARGU,I1,X1,IPCTRA(21),IRET)
C             CALL OPCHP1(IPCTRA(9),IOPERA,IARGU,I1,X1,IPCTRA(22),IRET)
c             CALL ADCHPO (IPCTRA(21),IPCTRA(22),IPCTRA(23),1.D0,1.D0)
c
c             CALL MOTS1 (IPLMOT,MOTCLE)
c             CALL MAXIM1 (IPCTRA(23),IPLMOT,MOTCLE,0,CXINF)
c
c *Norme euclidienne de lambda^2*Mx
c             CALL XTX1(IPCTRA(8),RMXTMX)
c             CALL XTX1(IPCTRA(10),IMXTMX)
c             MXTMX=RMXTMX+IMXTMX
c
c *Norme infinie de lambda^2*Mx
C             IOPERA = 14
C             IARGU = 0
C             I1    = 0
C             X1    = 0.D0
C             CALL OPCHP1(IPCTRA(8),IOPERA,IARGU,I1,X1,IPCTRA(24),IRET)
C             CALL OPCHP1(IPCTRA(10),IOPERA,IARGU,I1,X1,IPCTRA(25),IRET)
c             CALL ADCHPO (IPCTRA(24),IPCTRA(25),IPCTRA(26),1.D0,1.D0)
c
c             CALL MOTS1 (IPLMOT,MOTCLE)
c             CALL MAXIM1 (IPCTRA(26),IPLMOT,MOTCLE,0,MXINF)
 
          ENDIF
 
C Norme euclidienne du residu
          CALL XTX1(IPCTRA(10),RXTX)
          CALL XTX1(IPCTRA(14),IXTX)
          XTX=RXTX+IXTX
 
 
C Norme infinie du residu
          IOPERA = 14
          IARGU = 0
          I1    = 0
          X1    = 0.
          CALL OPCHP1(IPCTRA(10),IOPERA,IARGU,I1,X1,IPCTRA(15),IRET)
          CALL OPCHP1(IPCTRA(14),IOPERA,IARGU,I1,X1,IPCTRA(16),IRET)
 
          CALL ADCHPO(IPCTRA(15),IPCTRA(16),IPCTRA(17),1.D0,1.D0)
 
          CALL MOTS1 (IPLMOT,MOTCLE)
          CALL MAXIM1 (IPCTRA(17),IPLMOT,MOTCLE,0,XINF)
 
C Calcul des normes
            EUC=XTX
            INF=XINF
 
        ENDIF
 
      ENDIF
 
 
 
C ***********************************************************************
C                      AFFICHAGE DES SORTIES
C ***********************************************************************
 
      IF (MAL) THEN
        WRITE(IOIMP,*) 'LE MODE' , i , 'EST MAL NORMALISE ET DONC'
        WRITE(IOIMP,*) 'PEUT ETRE ERRONE. TOLERANCE ADMISSIBLE' , EPSI
        WRITE(IOIMP,*) 'VALEUR DE ||MODE|| :' , i , REAL(RAY)
      ENDIF
 
      IF (IIMPI .GT. 2) THEN
          WRITE(IOIMP,*) ' '
          WRITE(IOIMP,*) '**************************************'
          WRITE(IOIMP,*) '*ETUDE DU MODE',i,'*'
          WRITE(IOIMP,*) '**************************************'
          WRITE(IOIMP,*) ' '
          WRITE(IOIMP,*) ' '
          WRITE(IOIMP,*) 'Norme du mode :'
          WRITE(IOIMP,*) '---------------'
          WRITE(IOIMP,*) ' '
          WRITE(IOIMP,*) REAL(RAY)
          WRITE(IOIMP,*) ' '
          WRITE(IOIMP,*) ' '
 
          IF (SYM) THEN
            WRITE(IOIMP,*) 'Chpoint residu :'
            WRITE(IOIMP,*) '----------------'
            CALL ECCHPO (IPCTRA(10),1)
          ELSE
            WRITE(IOIMP,*) 'Partie reelle du chpoint residu :'
            WRITE(IOIMP,*) '---------------------------------'
            CALL ECCHPO (IPCTRA(10),1)
            WRITE(IOIMP,*) 'Partie imaginaire du chpoint residu : '
            WRITE(IOIMP,*) '--------------------------------------'
            CALL ECCHPO (IPCTRA(14),1)
          ENDIF
 
          WRITE(IOIMP,*) 'Norme euclidienne :',EUC
          WRITE(IOIMP,*) '-------------------'
          WRITE(IOIMP,*) ' '
          WRITE(IOIMP,*) 'Norme infinie :',INF
          WRITE(IOIMP,*) '---------------'
          WRITE(IOIMP,*) ' '
 
 
      ENDIF
 
 
C Destruction des chpoints de travail
      IF (IIMPI .GT. 2) THEN
        DO j=1,30
          IF (IPCTRA(j) .NE. 0) THEN
            CALL DTCHPO (IPCTRA(j))
          ENDIF
        ENDDO
      ENDIF
 
      SEGDES MRITRA
 
      END