1 : C RAYLE2 SOURCE CHAT 09/10/09 21:22:21 6519 2 : C RAYLEI SOURCE AM 08/02/14 21:28:07 6050 3 : SUBROUTINE RAYLE2 (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA 4 : & ,IPRX, IPIX, ICVG) 5 : IMPLICIT INTEGER(I-N) 6 : IMPLICIT REAL*8 (A-H, O-Z) 7 : 8 : ***************************************************************************** 9 : * * 10 : * RAYLEI * 11 : * ------ * 12 : ***************************************************************************** 13 : * 14 : * Fonction: 15 : * Trouve les valeurs propres d'une matrice de Rayleigh. 16 : * 17 : * création : Benoit Prabel mars 2009 18 : * 19 : * 20 : ***************************************************************************** 21 : * 21 : * PROCEDURE D APPEL: CALL RAYLEI (IPA, IPB, IV, IW, XRVP, XIVP,IPLMOX, 22 : * IPLMOY,DELTA, IPRX, IPIX) 23 : * 24 : -INC CCOPTIO 25 : -INC SMLMOTS 26 : -INC SMCHPOI 27 : -INC SMRIGID 28 : -INC CCHAMP 29 : 30 : REAL*8 A01(1:2,1:2),B01(1:2,1:2) 31 : REAL*8 etaR(1:2),etaI(1:2) 32 : REAL*8 XRVP, XIVP, XRES, alpha0, alpha1, alpha2, DELTA 33 : 34 : ICVG = 0 35 : 36 : *-------Projection du probleme sur la base de vecteurs (V, W) 36 : call MUCPRI (IV, IPA, IPAV) 36 : call MUCPRI (IW, IPA, IPAW) 36 : CALL MUCPRI (IV, IPB, IPBV) 36 : CALL MUCPRI (IW, IPB, IPBW) 36 : call XTY1 (IV, IPAV, IPLMOX, IPLMOY, A01(1,1) ) 36 : call XTY1 (IV, IPAW, IPLMOX, IPLMOY, A01(1,2) ) 36 : call XTY1 (IW, IPAV, IPLMOX, IPLMOY, A01(2,1) ) 36 : call XTY1 (IW, IPAW, IPLMOX, IPLMOY, A01(2,2) ) 36 : call XTY1 (IV, IPBV, IPLMOX, IPLMOY, B01(1,1) ) 36 : call XTY1 (IV, IPBW, IPLMOX, IPLMOY, B01(1,2) ) 36 : call XTY1 (IW, IPBV, IPLMOX, IPLMOY, B01(2,1) ) 36 : call XTY1 (IW, IPBW, IPLMOX, IPLMOY, B01(2,2) ) 37 : * write(*,*) A01(1,1),A01(1,2),'| -lambda |',B01(1,1),B01(1,2) 38 : * write(*,*) A01(2,1),A01(2,2),'| |',B01(2,1),B01(2,2) 39 : 40 : 41 : *-------%il faut annuler le determinant de la matrice |A01 - lambda B01| 42 : * % det(A01 - lambda B01) = alpha0 + alpha1 lambda + alpha2 lambda^2 43 : alpha0 = (A01(1,1)*A01(2,2)) - (A01(2,1)*A01(1,2)) 44 : alpha1 = (A01(2,1)*B01(1,2)) + (B01(2,1)*A01(1,2)) 45 : & - (A01(1,1)*B01(2,2)) - (B01(1,1)*A01(2,2)) 46 : alpha2 = (B01(1,1)*B01(2,2)) - (B01(2,1)*B01(1,2)) 47 : DELTA = (alpha1**2) - (4.*alpha0*alpha2) 48 : * write(*,*) 'Rayleigh trouve Delta=',DELTA 49 : 50 : *-------% selon DELTA, on en deduit le type de val propres: 51 : * %- 1 val propre Reelle double 52 : if( (abs(DELTA)) .lt. 1.D-10) then 53 : goto 200 54 : else 55 : * %- 2 val propres Reelles simples => il faut continuer a iterer 56 : if( DELTA .gt. 0.) then 57 : goto 300 58 : * %- 2 val propres Complexes conjuguées => on resout 59 : elseif( DELTA .lt. 0.) then 60 : goto 100 61 : else 62 : * %normalement, on n'arrive pas la... 63 : return 64 : endif 65 : endif 66 : 67 : 100 CONTINUE 68 : *-------%on se place dans le cas de 2 vp Complexes conjuguées XRVP +/- i*XIVP 69 : ICVG=1 70 : XRVP = -1.*alpha1 / (2.*alpha2) 71 : XIVP = ((abs(DELTA))**0.5) / (2.*alpha2) 72 : * %reste a calculer eta (=vecteur solution du pb projeté), 73 : * %puis les vect propres associés 74 : * % on fixe la 1ere valeur à 1 + 0*i 75 : etaR(1)=1. 76 : etaI(1)=0. 77 : * % on crée 78 : a11R = A01(1,1) - XRVP*B01(1,1) 79 : a11I = - XIVP*B01(1,1) 80 : a12R = A01(1,2) - XRVP*B01(1,2) 81 : a12I = - XIVP*B01(1,2) 82 : xdenR = (a12R**2) + (a12I**2) 83 : etaR(2) = -1.*(a11R*a12R + a11I*a12I) / xdenR 84 : etaI(2) = -1.*(a11I*a12R - a11R*a12I) / xdenR 85 : * % calcul du vect propre 85 : call ADCHPO (IV, IW, IPRX, etaR(1), etaR(2)) 85 : call ADCHPO (IV, IW, IPIX, etaI(1), etaI(2)) 86 : * % normalisation Complexe => nouvelle routine NORM1C 87 : * appel a MOTS1 copié de itinv1.eso pour exclure LX 87 : * CALL MOTS1 (IPLMOT,MOTCLE) 87 : CALL MOTS2 (IPLMOT,MOTCLE) 87 : CALL NORM1C (IPRX,IPIX,IPLMOT,MOTCLE,IPRX,IPIX) 88 : * % c'est fini ! 89 : goto 900 90 : 91 : 200 CONTINUE 92 : *-------%on se place dans le cas de 1 val propre Reelle double XRVP 93 : ICVG=2 94 : XRVP = -1.*alpha1 / (2.*alpha2) 95 : XIVP = 0.D0 96 : * % on met les 2 vecteurs orthonormée Réels dans IPRX et IPIX !!! 97 : * % => il faudra donc bien testé DELTA dans le programme 98 : * % qui appelle Rayleigh pour différencier du cas complexe !!! 98 : * CALL MOTS1 (IPLMOT,MOTCLE) 98 : CALL MOTS2 (IPLMOT,MOTCLE) 98 : call NORMA1 (IV,IPLMOT,MOTCLE,IPRX) 98 : call NORMA1 (IW,IPLMOT,MOTCLE,IPIX) 99 : goto 900 100 : * 101 : 300 CONTINUE 102 : *-------%on se place dans le cas de 2 val propres Reelles simples distinctes 103 : XRVP = ((DELTA**0.5) - alpha1) / (2.*alpha2) 104 : XRVP2 = -1.*( (DELTA**0.5) + alpha1 ) / (2.*alpha2) 105 : * write(*,*) XRVP,XRVP2 106 : if( abs(1. - abs(XRVP/XRVP2)) .gt. 1.D-5) then 107 : * % => il faut donc continuer a iterer pour les séparer (ou les confondre) + nettement 108 : ICVG=0 109 : XRVP = 0.D0 110 : XIVP = 0.D0 111 : IPRX = 0 112 : IPIX = 0 113 : else 114 : * %on a 2 val p reelle simple de meme module -> on resout 115 : * -> non! on prefere n'en donner qu'une... 116 : ICVG=3 117 : * XIVP = XRVP2 118 : if( (XRVP .lt. 0.) .and. (XRVP2 .gt. 0.)) XRVP=XRVP2 119 : XIVP = 0.D0 120 : * % on fixe la 1ere valeur du 1er vecteur propre à 1 + 0*i 121 : etaR(1) = 1. 122 : a11R = A01(1,1) - XRVP*B01(1,1) 123 : a12R = A01(1,2) - XRVP*B01(1,2) 124 : if( abs(a12R) .lt. 1.E-10 ) then 125 : etaR(1) = 0. 126 : etaR(2) = 1. 127 : else 128 : etaR(2) = -1.*(a11R/a12R) 129 : endif 129 : call ADCHPO (IV, IW, IPRX, etaR(1), etaR(2)) 129 : * CALL MOTS1 (IPLMOT,MOTCLE) 129 : CALL MOTS2 (IPLMOT,MOTCLE) 129 : CALL NORMA1 (IPRX,IPLMOT,MOTCLE,IPRX) 130 : * % on fixe la 1ere valeur du 2eme vecteur propre à 1 + 0*i 131 : * etaR(1) = 1. 132 : * a11R = A01(1,1) - XRVP2*B01(1,1) 133 : * a12R = A01(1,2) - XRVP2*B01(1,2) 134 : * if( abs(a12R) .lt. 1.E-10 ) then 135 : * etaR(1) = 0. 136 : * etaR(2) = 1. 137 : * else 138 : * etaR(2) = -1.*(a11R/a12R) 139 : * endif 139 : * call ADCHPO (IV, IW, IPRX2, etaR(1), etaR(2)) 139 : * CALL MOTS1 (IPLMOT,MOTCLE) 139 : * CALL NORMA1 (IPRX2,IPLMOT,MOTCLE,IPRX2) 140 : * IPIX = IPRX2 141 : IPIX = 0 142 : endif 143 : 144 : 900 CONTINUE 145 : 146 : END 147 : 148 : 149 :
© Cast3M 2003 - All rights reserved.
Disclaimer