Sources Esope | Cast3M

sespc5
C SESPC5    SOURCE    CB215821  20/11/25    13:39:46     10792          
*
************************************************************************
* SESPC5
************************************************************************
* FONCTION:
*     ON TEST LA CONVERGENCE DES ELEMENTS PROPRES Complexes:
*                1./  TEST entre 2 itérés de LA VALEUR PROPRE complexe.
*                2./  TEST SUR Le Residu = [A - lambda B]*X.
*
* AUTEUR, DATE DE CREATION:
*   Benoit Prabel (Novembre 2008)
*
************************************************************
 
      SUBROUTINE SESPC5 ( ILVA1R,ILVA1I, ILVA2R,ILVA2I, IPLCH2,
     1                    IPA,IPB, BOOL1,BOOL2, NBMOD )
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC PPARAM
-INC CCOPTIO
-INC  SMLREEL
-INC  SMLCHPO
 
******
*     -- CONSTANTES --
***
      PARAMETER ( PRECI1 = 1.0D-6 )
      PARAMETER ( PRECI2 = 1.0D-6 )
 
******
*     -- ARGUMENTS --
***
      POINTEUR ILVA1R.MLREEL,ILVA1I.MLREEL,ILVA2R.MLREEL,ILVA2I.MLREEL
      POINTEUR IPLCH1.MLCHPO, IPLCH2.MLCHPO
      INTEGER  NBMOD, IPA,IPB, IPAX,IPBX
      LOGICAL  BSWAP,BOOL1,BOOL2
 
******
*     -- VARIABLES LOCALES --
***
      INTEGER  IB100
      REAL*8   XTMP, ALPHA1, ALPHA2, RESMAX,RESMAR,RESMAI
      CHARACTER*(LOCOMP) MOTCLE
 
 
*---- ON TRIE -----------------------------------------------------
*     le tri ne concerne que la liste 2 + le chpoint associé,
*     car on suppose la liste 1 deja ordonnée
cbp,2019       segact,ILVA2R,ILVA2I
      segact,ILVA2R*mod
      segact,ILVA2I*mod
      segact,IPLCH2*mod
      ILDIM = ILVA2R.PROG(/1)
 
*     algo type tri a bulle: pas tres efficace mais simple (cf ORDVEC.eso)
 10   CONTINUE
      BSWAP = .FALSE.
 
*     on boucle sur toutes les valeurs de la liste
      DO 100 IB100=1,(ILDIM-1)
 
*       on recupere la valeur test (la ib100 ième)
        XR1 = ILVA2R.PROG(IB100)
        XI1 = ILVA2I.PROG(IB100)
        XM1 = ((XR1**2) + (XI1**2))**0.5
        XR2 = ILVA2R.PROG(IB100 + 1)
        XI2 = ILVA2I.PROG(IB100 + 1)
        XM2 = ((XR2**2) + (XI2**2))**0.5
 
*       si le module de 2 val propre est trop proche, alors on trie selon la partie réelle
        if( (abs (XM1 - XM2)) .lt. (1.0D-5 * XM2)) then
          XM1 = XR1
          XM2 = XR2
        endif
 
*       cas ou la valeur est mal placée
        if(XM1 .gt. XM2) then
cbp,2019            segdes,ILVA2R,ILVA2I
           CALL SWAP( IB100, ILVA2R, IPLCH2 )
cbp,2019            segact,ILVA2I*mod
           XTMP = ILVA2I.PROG( IB100 )
           ILVA2I.PROG( IB100 ) = ILVA2I.PROG( IB100 + 1 )
           ILVA2I.PROG( IB100 + 1 ) = XTMP
           IF ( IERR .NE. 0 ) RETURN
           BSWAP = .TRUE.
cbp,2019           segact,ILVA2R
        endif
 
 100  CONTINUE
 
      IF(BSWAP)  GOTO 10
 
 
 
*---- TESTS DE CONVERGENCE sur les NBMOD 1ers modes ----------------
      CALL MOTS1 (IPLMOT,MOTCLE)
      BOOL1 = .true.
      BOOL2 = .true.
      JVEC = 0
      segact,ILVA1R,ILVA1I
cbp,2019      segact,IPLCH2
 
      DO 200 IB200 =1,NBMOD
 
        JVEC = JVEC + 1
 
*------ CONVERGENCE DE LA VALEUR PROPRE ----------------------------
 
        XR1 = ILVA1R.PROG(JVEC)
        XI1 = ILVA1I.PROG(JVEC)
        XR2 = ILVA2R.PROG(JVEC)
        XI2 = ILVA2I.PROG(JVEC)
 
*       pour l'instant on ne juge pas necessaire de tester la vitesse de cvgce de lambda
*       du moins pas de cette maniere qui utilise 2 itérés parfois très lointains....
        goto 201
 
*        write(*,*) 'test de la val p',JVEC,':',XR1,XI1,' ?=? ',XR2,XI2,
*     1  '  ecart=',(XR2-XR1),(XI2-XI1)
        XREFR = max(PRECI1,(abs(XR2)))
        XREFI = max(PRECI1,(abs(XI2)))
        IF(    ( (abs(XR2 - XR1)) .GT. (PRECI1*XREFR) )
     1    .OR. ( (abs(XI2 - XI1)) .GT. (PRECI1*XREFI) ) ) THEN
            BOOL1 = .FALSE.
            GOTO 900
        ENDIF
 
 201  CONTINUE
*------ CONVERGENCE DU VECTEUR PROPRE (Residu nul) ------------------
 
*       Cas d'un mode Réel
        IF(XI2 .eq. 0.D0) THEN
*           write(*,*) 'on va tester le vecteur reel num',JVEC
           IPCHP2 = IPLCH2.ICHPOI( JVEC )
           CALL MUCPRI ( IPCHP2, IPA, IPAX )
           CALL MUCPRI ( IPCHP2, IPB, IPBX )
           CALL COMBI2 ( IPAX, 1.D0, IPBX, (-1.D0*XR2), IRES )
           CALL DTCHPO (IPAX)
           CALL DTCHPO (IPBX)
           CALL MAXIM1 (IRES,IPLMOT,MOTCLE,0,RESMAX)
           CALL DTCHPO (IRES)
           RESMAX = abs (RESMAX)
*           write(*,*) 'mode reel JVEC,RESMAX=',JVEC,RESMAX
 
 
*       Cas d'un mode Complexe
        ELSE
*           write(*,*) 'on va tester le vecteur complexe num',JVEC
           ICHP2R = IPLCH2.ICHPOI( JVEC )
           ICHP2I = IPLCH2.ICHPOI( JVEC + 1 )
           CALL COMBI2 (ICHP2R,(-1.D0*XR2),ICHP2I,XI2,ICHP3R)
           CALL COMBI2 (ICHP2R,(-1.D0*XI2),ICHP2I,(-1.D0*XR2),ICHP3I)
           CALL MUCPRI ( ICHP2R, IPA, IPAXR )
           CALL MUCPRI ( ICHP3R, IPB, IPBXR )
           CALL DTCHPO (ICHP3R)
           CALL MUCPRI ( ICHP2I, IPA, IPAXI )
           CALL MUCPRI ( ICHP3I, IPB, IPBXI )
           CALL DTCHPO (ICHP3I)
           CALL COMBI2 ( IPAXR, 1.D0, IPBXR, 1.D0, IRESR )
           CALL DTCHPO (IPAXR)
           CALL DTCHPO (IPBXR)
           CALL COMBI2 ( IPAXI, 1.D0, IPBXI, 1.D0, IRESI )
           CALL DTCHPO (IPAXI)
           CALL DTCHPO (IPBXI)
           CALL MAXIM1 (IRESR,IPLMOT,MOTCLE,0,RESMAR)
           CALL DTCHPO (IRESR)
           CALL MAXIM1 (IRESI,IPLMOT,MOTCLE,0,RESMAI)
           CALL DTCHPO (IRESI)
           RESMAX = (abs(RESMAR)) + (abs(RESMAI))
*           write(*,*) 'mode complexe JVEC,RESMAX=',JVEC,RESMAX
           JVEC = JVEC + 1
        ENDIF
 
        IF( IERR .NE. 0 )  RETURN
 
*       Test de convergence sur le residu
        IF(RESMAX .ge. PRECI2) then
           BOOL2 = .FALSE.
           GOTO 900
        ENDIF
 
        if(JVEC .ge. NBMOD)  goto 900
 
 
 200  CONTINUE
 
 900  CONTINUE
      SEGDES ,ILVA1R,ILVA1I, ILVA2R,ILVA2I, IPLCH2
*     on modifie le nbmod proposé si le dernier mode est complexe
*     et que l'on a besoin du nbmod+1^ieme vecteur
      if(JVEC .gt. NBMOD)   NBMOD=JVEC
 
      RETURN
      END