Sources Esope | Cast3M

rayle2
C RAYLE2    SOURCE    CB215821  20/11/25    13:38:31     10792          
      SUBROUTINE RAYLE2 (IPA,IPB,IV,IW,XRVP,XIVP,IPLMOX,IPLMOY,DELTA
     &      ,IPRX, IPIX, ICVG)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H, O-Z)
 
*****************************************************************************
*                                                                           *
*                               RAYLEI                                      *
*                               ------                                      *
*****************************************************************************
*
*       Fonction:
*       Trouve les valeurs propres d'une matrice de Rayleigh.
*
*       création :              Benoit Prabel    mars 2009
*
*
*****************************************************************************
*
*       PROCEDURE D APPEL: CALL RAYLEI (IPA, IPB, IV, IW, XRVP, XIVP,IPLMOX,
*                                      IPLMOY,DELTA, IPRX, IPIX)
*
 
-INC PPARAM
-INC CCOPTIO
-INC SMLMOTS
-INC SMCHPOI
-INC SMRIGID
-INC CCHAMP
 
        REAL*8 A01(1:2,1:2),B01(1:2,1:2)
        REAL*8 etaR(1:2),etaI(1:2)
        REAL*8 XRVP, XIVP, XRES, alpha0, alpha1, alpha2, DELTA
        CHARACTER*(LOCOMP) MOTCLE
 
        ICVG = 0
 
*-------Projection du probleme sur la base de vecteurs (V, W)
        call MUCPRI (IV, IPA, IPAV)
        call MUCPRI (IW, IPA, IPAW)
        CALL MUCPRI (IV, IPB, IPBV)
        CALL MUCPRI (IW, IPB, IPBW)
        call XTY1 (IV, IPAV, IPLMOX, IPLMOY, A01(1,1) )
        call XTY1 (IV, IPAW, IPLMOX, IPLMOY, A01(1,2) )
        call XTY1 (IW, IPAV, IPLMOX, IPLMOY, A01(2,1) )
        call XTY1 (IW, IPAW, IPLMOX, IPLMOY, A01(2,2) )
        call XTY1 (IV, IPBV, IPLMOX, IPLMOY, B01(1,1) )
        call XTY1 (IV, IPBW, IPLMOX, IPLMOY, B01(1,2) )
        call XTY1 (IW, IPBV, IPLMOX, IPLMOY, B01(2,1) )
        call XTY1 (IW, IPBW, IPLMOX, IPLMOY, B01(2,2) )
*        write(*,*) A01(1,1),A01(1,2),'| -lambda |',B01(1,1),B01(1,2)
*        write(*,*) A01(2,1),A01(2,2),'|         |',B01(2,1),B01(2,2)
 
 
*-------%il faut annuler le determinant de la matrice |A01 - lambda B01|
*       % det(A01 - lambda B01) = alpha0 + alpha1 lambda + alpha2 lambda^2
        alpha0 = (A01(1,1)*A01(2,2)) - (A01(2,1)*A01(1,2))
        alpha1 = (A01(2,1)*B01(1,2)) + (B01(2,1)*A01(1,2))
     &         - (A01(1,1)*B01(2,2)) - (B01(1,1)*A01(2,2))
        alpha2 = (B01(1,1)*B01(2,2)) - (B01(2,1)*B01(1,2))
        DELTA = (alpha1**2) - (4.*alpha0*alpha2)
*        write(*,*) 'Rayleigh trouve Delta=',DELTA
 
*-------% selon DELTA, on en deduit le type de  val propres:
*       %- 1 val propre Reelle double
        if( (abs(DELTA)) .lt. 1.D-10) then
            goto 200
        else
*        %- 2 val propres Reelles simples => il faut continuer a iterer
         if( DELTA .gt. 0.) then
            goto 300
*        %- 2 val propres Complexes conjuguées => on resout
         elseif( DELTA .lt. 0.) then
            goto 100
         else
*            %normalement, on n'arrive pas la...
            return
         endif
        endif
 
 100    CONTINUE
*-------%on se place dans le cas de 2 vp Complexes conjuguées XRVP +/- i*XIVP
        ICVG=1
        XRVP = -1.*alpha1 / (2.*alpha2)
        XIVP = ((abs(DELTA))**0.5) / (2.*alpha2)
*       %reste a calculer eta (=vecteur solution du pb projeté),
*       %puis les vect propres associés
*       % on fixe la 1ere valeur à 1 + 0*i
        etaR(1)=1.
        etaI(1)=0.
*        % on crée
        a11R = A01(1,1) - XRVP*B01(1,1)
        a11I =          - XIVP*B01(1,1)
        a12R = A01(1,2) - XRVP*B01(1,2)
        a12I =          - XIVP*B01(1,2)
        xdenR = (a12R**2) + (a12I**2)
        etaR(2) = -1.*(a11R*a12R + a11I*a12I) / xdenR
        etaI(2) = -1.*(a11I*a12R - a11R*a12I) / xdenR
*       % calcul du vect propre
        call ADCHPO (IV, IW, IPRX, etaR(1),  etaR(2))
        call ADCHPO (IV, IW, IPIX, etaI(1),  etaI(2))
*       % normalisation Complexe => nouvelle routine NORM1C
*       appel a MOTS1 copié de itinv1.eso pour exclure LX
*        CALL MOTS1 (IPLMOT,MOTCLE)
        CALL MOTS2 (IPLMOT,MOTCLE)
        CALL NORM1C (IPRX,IPIX,IPLMOT,MOTCLE,IPRX1,IPIX1)
        IPRX=IPRX1
        IPIX=IPIX1
*       % c'est fini !
        goto 900
 
 200    CONTINUE
*-------%on se place dans le cas de 1 val propre Reelle double XRVP
        ICVG=2
        XRVP = -1.*alpha1 / (2.*alpha2)
        XIVP = 0.D0
*       % on met les 2 vecteurs orthonormée Réels dans IPRX et IPIX !!!
*       % => il faudra donc bien testé DELTA dans le programme
*       % qui appelle Rayleigh pour différencier du cas complexe !!!
*        CALL MOTS1 (IPLMOT,MOTCLE)
        CALL MOTS2 (IPLMOT,MOTCLE)
        call NORMA1 (IV,IPLMOT,MOTCLE,IPRX)
        call NORMA1 (IW,IPLMOT,MOTCLE,IPIX)
        goto 900
*
 300    CONTINUE
*-------%on se place dans le cas de 2 val propres Reelles simples distinctes
        XRVP = ((DELTA**0.5) - alpha1) / (2.*alpha2)
        XRVP2 = -1.*( (DELTA**0.5) + alpha1 ) /  (2.*alpha2)
*        write(*,*) XRVP,XRVP2
        if( abs(1. - abs(XRVP/XRVP2)) .gt. 1.D-5) then
*          % => il faut donc continuer a iterer pour les séparer (ou les confondre) + nettement
           ICVG=0
           XRVP = 0.D0
           XIVP = 0.D0
           IPRX = 0
           IPIX = 0
        else
*          %on a 2 val p reelle simple de meme module -> on resout
*           -> non! on prefere n'en donner qu'une...
           ICVG=3
*           XIVP = XRVP2
           if( (XRVP .lt. 0.) .and. (XRVP2 .gt. 0.))    XRVP=XRVP2
           XIVP = 0.D0
*          % on fixe la 1ere valeur du 1er vecteur propre à  1 + 0*i
            etaR(1) =  1.
            a11R = A01(1,1) - XRVP*B01(1,1)
            a12R = A01(1,2) - XRVP*B01(1,2)
            if( abs(a12R) .lt. 1.E-10 ) then
               etaR(1) =  0.
               etaR(2) =  1.
            else
               etaR(2) = -1.*(a11R/a12R)
            endif
            call ADCHPO (IV, IW, IPRX, etaR(1),  etaR(2))
*            CALL MOTS1  (IPLMOT,MOTCLE)
            CALL MOTS2  (IPLMOT,MOTCLE)
            CALL NORMA1 (IPRX,IPLMOT,MOTCLE,IPRX)
*          % on fixe la 1ere valeur du 2eme vecteur propre à  1 + 0*i
*            etaR(1) =  1.
*            a11R = A01(1,1) - XRVP2*B01(1,1)
*            a12R = A01(1,2) - XRVP2*B01(1,2)
*            if( abs(a12R) .lt. 1.E-10 ) then
*               etaR(1) =  0.
*               etaR(2) =  1.
*            else
*               etaR(2) = -1.*(a11R/a12R)
*            endif
*            call ADCHPO (IV, IW, IPRX2, etaR(1),  etaR(2))
*            CALL MOTS1  (IPLMOT,MOTCLE)
*            CALL NORMA1 (IPRX2,IPLMOT,MOTCLE,IPRX2)
*            IPIX = IPRX2
            IPIX = 0
        endif
 
 900    CONTINUE
 
        END