1 : C RAYLEI SOURCE CHAT 09/10/09 21:22:23 6519 2 : SUBROUTINE RAYLEI (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA 3 : & ,IPRX, IPIX) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8 (A-H, O-Z) 6 : 7 : ***************************************************************************** 8 : * * 9 : * RAYLEI * 10 : * ------ * 11 : ***************************************************************************** 12 : * 13 : * Fonction: 14 : * 15 : * Trouve les valeurs propres d'une matrice de Rayleigh. 16 : * 17 : ***************************************************************************** 18 : * 19 : * 20 : * Equation caracteristique: (1, -(LL11+LL22), (LL11*LL22-LL12*LL21)) 21 : * Discriminant: Delta 22 : * 23 : * 23 : * PROCEDURE D APPEL: CALL RAYLEI (IPA, IPB, IV, IW, XRVP, XIVP,IPLMOX, 24 : * IPLMOY,DELTA, IPRX, IPIX) 25 : * 26 : -INC CCOPTIO 27 : -INC SMLMOTS 28 : -INC SMCHPOI 29 : -INC SMRIGID 30 : -INC CCHAMP 31 : SEGMENT IDEMEN (0) 32 : 33 : REAL*8 LL(1:2,1:2) 34 : REAL*8 A, B, C, D, TO, TI, XRVP, XIVP 35 : * COORDONNEES DE LA MATRICE 36 : MCHPOI = 0 37 : MCHPO2 = 0 37 : CALL MUCPRI (IV, IPB, IPBV) 38 : SEGINI IDEMEN 39 : IDEMEN(**)=IPBV 40 : NOID=0 41 : NOEN=1 41 : CALL LDMT (IPA,IDEMEN,NOID,NOEN,1.D-15) 42 : SEGACT IDEMEN 43 : MCHPOI=IDEMEN(1) 44 : SEGSUP IDEMEN 45 : IF (IERR .NE. 0) RETURN 45 : CALL XTY1 (IV, MCHPOI, IPLMOX, IPLMOX, LL(1,1) ) 46 : 46 : CALL MUCPRI (IW, IPB, IPBW) 47 : SEGINI IDEMEN 48 : IDEMEN(**)=IPBW 49 : NOID=0 50 : NOEN=1 50 : CALL LDMT (IPA,IDEMEN,NOID,NOEN,1.D-15) 51 : SEGACT IDEMEN 52 : MCHPO2=IDEMEN(1) 53 : SEGSUP IDEMEN 54 : IF (IERR .NE. 0) RETURN 54 : CALL XTY1 (IV, MCHPO2, IPLMOX, IPLMOX, LL(1,2) ) 54 : CALL XTY1 (IW, MCHPOI, IPLMOX, IPLMOX, LL(2,1) ) 54 : CALL XTY1 (IW, MCHPO2, IPLMOX, IPLMOX, LL(2,2) ) 55 : 56 : A=LL(1,1) 57 : B=LL(1,2) 58 : C=LL(2,1) 59 : D=LL(2,2) 60 : TO=(A+D)/2. 61 : 62 : * Discriminant 63 : DELTA=(A+D)**2-4.D0*(A*D-B*C) 64 : TI=ABS(DELTA) 65 : IF (TI .GE. 1.D-14) THEN 66 : IF (DELTA .LT. 0.D0) THEN 67 : GOTO 100 68 : ELSE 69 : RETURN 70 : END IF 71 : ELSE 72 : GOTO 200 73 : END IF 74 : 75 : 76 : 77 : 100 CONTINUE 78 : 79 : XRVP=TO 80 : XIVP=SQRT(TI)/2.D0 81 : * 82 : * EXTRACTION DES VECTEURS PROPRES 83 : * 84 : IF ( ABS( B ) .GE. ABS( C )) THEN 85 : X2R = (XRVP - A ) / B 86 : X2C = XIVP / B 86 : CALL ADCHPO (IV, IW, IPRX, 1.D0, X2R) 86 : CALL MUCHPO (IW, X2C, IPIX, 1) 87 : ELSE 88 : X1R = (XRVP - D ) / C 89 : X1C = XIVP / C 89 : CALL ADCHPO (IV, IW, IPRX, X1R, 1.D0) 89 : CALL MUCHPO (IV, X1C, IPIX, 1) 90 : END IF 91 : * 92 : * COEFFICIENT DE RAYLEIGH 93 : * EXTRACTION DES VALEURS PROPRES 94 : * 94 : CALL XTMX(IPRX, IPB, RXBRX) 94 : CALL XTMX(IPIX, IPB, RIXBRIX) 95 : ALPHA = RXBRX + RIXBRIX 95 : CALL XTMX(IPRX, IPA, RXARX) 95 : CALL XTMX(IPIX, IPA, RIXARIX) 96 : BETA = RXARX + RIXARIX 97 : XRVP = BETA / ALPHA 97 : CALL MUCPRI(IPIX, IPA, IPAXI) 97 : CALL XTY1(IPRX, IPAXI, IPLMOX, IPLMOY, XRKXI) 97 : CALL MUCPRI(IPRX, IPA, IPAXR) 97 : CALL XTY1(IPIX, IPAXR, IPLMOX, IPLMOY, XIKXR) 98 : XIVP = (XIKXR - XRKXI) / ALPHA 99 : GOTO 300 100 : 200 CONTINUE 101 : IPRX = IV 101 : CALL MUCPRI(IV, IPA, IPAX2) 101 : CALL XTY1(IV, IPAX2, IPLMOX, IPLMOY, X2AX2) 102 : IF (IERR .NE. 0 ) RETURN 102 : CALL MUCPRI(IV, IPB, IPBX2) 102 : CALL XTY1(IV, IPBX2, IPLMOX, IPLMOY, X2BX2) 103 : IF (IERR .NE. 0 ) RETURN 104 : XRVP = X2AX2 / X2BX2 105 : XIVP = 0.D0 106 : 300 CONTINUE 107 : END 108 : 109 : 110 : 111 : 112 : 113 : 114 :
© Cast3M 2003 - All rights reserved.
Disclaimer