C ARPVER SOURCE CB215821 25/04/08 21:15:04 12227 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.D0 RAY=CMPLX(0.D0,0.D0) 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.D0) 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.D0 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.D0) .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.D0 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