Télécharger la plaquette Cast3M


fr en

Télécharger qzvalp.eso

Retour à la liste

Numérotation des lignes : 

qzvalp

  1. C QZVALP    SOURCE    CB215821  17/11/30    21:17:05     9639           
  2.        SUBROUTINE QZVALP (MATA, MATB, ALFR, ALFI, BETA, CALV, MATZ)
  3.       IMPLICIT INTEGER(I-N)
  4.       IMPLICIT REAL*8 (A-H,O-Z)
  5. *
  6. *************************************************************************
  7. *                                                                       *
  8. *      FIN DE MISE SOUS FORME QUASI-TRIANGULAIRE SUPERIEURE DE MATA     *
  9. *           ET   TRIANGULARISATION SUPERIEURE SIMULTANEE DE MATB        *
  10. *                    EXTRACTION DES VALEURS PROPRES                     *
  11. *   ( TROISIEME PAS DE LA METHODE "QZ" POUR LE CALCUL MODAL COMPLEXE )  *
  12. * _____________________________________________________________________ *
  13. *                                                                       *
  14. *                                                                       *
  15. *     DATE : le 7 Avril 1995                                            *
  16. *     AUTEURS : J.ANTUNES , Nicolas BENECH                              *
  17. * _____________________________________________________________________ *
  18. *                                                                       *
  19. *    MODULE(S) APPELANT(S) : VIBRAC                                     *
  20. *                                                                       *
  21. *    MODULE(S) APPELE(S) :                                              *
  22. * _____________________________________________________________________ *
  23. *                                                                       *
  24. *    EN ENTREE :                                                        *
  25. *       - MATA : matrice reelle quasi triangulaire superieure (XMATRI)  *
  26. *       - MATB : matrice reelle triangulaire superieure (XMATRI)        *
  27. *       - CALV : logique indiquant si on calcule les vecteurs propres   *
  28. *       - MATZ : produit des transformations effectuees auparavant      *
  29. *                                ( si CALV = .TRUE. )                   *
  30. *                                                                       *
  31. *    EN SORTIE :                                                        *
  32. *       - MATA : matrice quasi-triangulaire superieure                  *
  33. *       - MATB : matrice sous forme triangulaire superieure             *
  34. *       - ALFR, ALFI, BETA : liste des valeurs propres                  *
  35. *            (parties reelle et imaginaire du numerateur + denominateur)*
  36. *       - MATZ : produit des transformations effectuees (si CALV=.TRUE.)*
  37. *                                                                       *
  38. *************************************************************************
  39. *
  40. -INC PPARAM
  41. -INC CCOPTIO
  42. -INC SMRIGID
  43. -INC SMLREEL
  44. C
  45.       REAL*8 SWAP
  46.       INTEGER I , J , EN , NA , NN , ISW
  47. C
  48.       LOGICAL CALV
  49. C-----------------------------------------------------------------------
  50. C
  51.       POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI
  52. cbp      POINTEUR ALFR.XMATRI, ALFI.XMATRI, BETA.XMATRI
  53.       POINTEUR ALFR.MLREEL, ALFI.MLREEL, BETA.MLREEL
  54. C
  55. c       SEGACT , MATA*mod
  56. c       SEGACT , MATB*mod
  57. c       SEGACT , MATZ*mod
  58.       NM=MATA.RE(/1)
  59.       N=MATA.RE(/2)
  60. C
  61.       JG = N
  62.       SEGINI,ALFR,ALFI,BETA
  63. C-----------------------------------------------------------------------
  64. C========================= DEBUT DES CALCULS ===========================
  65. C-----------------------------------------------------------------------
  66. C
  67.       EPSB = MATB.RE(N,1,1)
  68.       ISW = 1
  69. C
  70. C-----------------------------------------------------------------------
  71. C-------------------------- BOUCLE PRINCIPALE --------------------------
  72. C-----------------------------------------------------------------------
  73. C
  74.       DO 510 NN = 1, N
  75.          EN = N + 1 - NN
  76.          NA = EN - 1
  77.          IF ( ISW . EQ . 2 ) GOTO 505
  78.          IF ( EN . EQ . 1 ) GOTO 410
  79.          SWAP = MATA.RE(EN,NA,1)
  80.          IF ( SWAP . NE . 0.0D0 ) GOTO 420
  81. C
  82. C-----------------------------------------------------------------------
  83. C---------------------- UNE VALEUR PROPRE REELLE -----------------------
  84. C-----------------------------------------------------------------------
  85. C
  86.   410    continue
  87.          ALFR.PROG(EN) = MATA.RE(EN,EN,1)
  88.          IF ( MATB.RE(EN,EN,1) . LT . 0.0D0 )
  89.      &     ALFR.PROG(EN) = -ALFR.PROG(EN)
  90.          BETA.PROG(EN) = ABS(MATB.RE(EN,EN,1))
  91.          ALFI.PROG(EN) = 0.0D0
  92.          GOTO 510
  93. C
  94.   420    continue
  95.          SWAP = ABS(MATB.RE(NA,NA,1))
  96.          IF ( SWAP . LE . EPSB ) GOTO 455
  97.          SWAP = ABS(MATB.RE(EN,EN,1))
  98.          IF ( SWAP . GT . EPSB ) GOTO 430
  99.          A1 = MATA.RE(EN,EN,1)
  100.          A2 = MATA.RE(EN,NA,1)
  101.          BN = 0.0D0
  102.          GOTO 435
  103.   430    continue
  104.          AN = ABS(MATA.RE(NA,NA,1)) + ABS(MATA.RE(NA,EN,1))
  105.      &      + ABS(MATA.RE(EN,NA,1)) + ABS(MATA.RE(EN,EN,1))
  106.          BN = ABS(MATB.RE(NA,NA,1)) + ABS(MATB.RE(NA,EN,1))
  107.      &      + ABS(MATB.RE(EN,EN,1))
  108.          A11 = MATA.RE(NA,NA,1) / AN
  109.          A12 = MATA.RE(NA,EN,1) / AN
  110.          A21 = MATA.RE(EN,NA,1) / AN
  111.          A22 = MATA.RE(EN,EN,1) / AN
  112.          B11 = MATB.RE(NA,NA,1) / BN
  113.          B12 = MATB.RE(NA,EN,1) / BN
  114.          B22 = MATB.RE(EN,EN,1) / BN
  115.          E = A11 / B11
  116.          EI = A22 / B22
  117.          S = A21 / (B11 * B22)
  118.          T = (A22 - E * B22) / B22
  119.          IF ( ABS(E) . LE . ABS(EI) ) GOTO 431
  120.          E = EI
  121.          T = (A11 - E * B11) / B11
  122.   431    C = 0.5D0 * (T - S * B12)
  123.          D = C * C + S * (A12 - E * B12)
  124.          IF ( D . LT . 0.0D0 ) GOTO 480
  125. C
  126. C-----------------------------------------------------------------------
  127. C-------------------- DEUX VALEURS PROPRES REELLES ---------------------
  128. C-----------------------------------------------------------------------
  129. C
  130.          E = E + (C +SIGN(SQRT(D),C))
  131.          A11 = A11 - E * B11
  132.          A12 = A12 - E * B12
  133.          A22 = A22 - E * B22
  134.          IF ( ABS(A11) + ABS(A12) . LT .
  135.      &        ABS(A21) + ABS(A22) ) GOTO 432
  136.          A1 = A12
  137.          A2 = A11
  138.          GOTO 435
  139.   432    A1 = A22
  140.          A2 = A21
  141. C
  142. C-----------------------------------------------------------------------
  143. C------------------------- CHOISIR  Z  REEL ----------------------------
  144. C-----------------------------------------------------------------------
  145. C
  146.   435    S = ABS(A1) + ABS(A2)
  147.          U1 = A1 / S
  148.          U2 = A2 / S
  149.          R =SIGN(SQRT(U1*U1+U2*U2),U1)
  150.          V1 = -(U1 + R) / R
  151.          V2 = -U2 / R
  152.          U2 = V2 / V1
  153. C
  154.          DO 440 I = 1, EN
  155.             T = MATA.RE(I,EN,1) + U2 * MATA.RE(I,NA,1)
  156.             MATA.RE(I,EN,1) = MATA.RE(I,EN,1) + T * V1
  157.             MATA.RE(I,NA,1) = MATA.RE(I,NA,1) + T * V2
  158.             T = MATB.RE(I,EN,1) + U2 * MATB.RE(I,NA,1)
  159.             MATB.RE(I,EN,1) = MATB.RE(I,EN,1) + T * V1
  160.             MATB.RE(I,NA,1) = MATB.RE(I,NA,1) + T * V2
  161.   440    CONTINUE
  162. C
  163.          IF ( .NOT. CALV ) GOTO 450
  164. C
  165.          DO 445 I = 1, N
  166.             T = MATZ.RE(I,EN,1) + U2 * MATZ.RE(I,NA,1)
  167.             MATZ.RE(I,EN,1) = MATZ.RE(I,EN,1) + T * V1
  168.             MATZ.RE(I,NA,1) = MATZ.RE(I,NA,1) + T * V2
  169.   445    CONTINUE
  170. C
  171.   450    IF ( BN . EQ . 0.0D0 ) GOTO 475
  172.          IF ( AN . LT . ABS(E) * BN ) GOTO 455
  173.          A1 = MATB.RE(NA,NA,1)
  174.          A2 = MATB.RE(EN,NA,1)
  175.          GOTO 460
  176.   455    continue
  177.          A1 = MATA.RE(NA,NA,1)
  178.          A2 = MATA.RE(EN,NA,1)
  179. C
  180. C-----------------------------------------------------------------------
  181. C------------------------- CHOISIR  Q  REEL ----------------------------
  182. C-----------------------------------------------------------------------
  183. C
  184.   460    S = ABS(A1) + ABS(A2)
  185.          IF ( S . EQ . 0.0D0 ) GOTO 475
  186.          U1 = A1 / S
  187.          U2 = A2 / S
  188.          R =SIGN(SQRT(U1*U1+U2*U2),U1)
  189.          V1 = -(U1 + R) / R
  190.          V2 = -U2 / R
  191.          U2 = V2 / V1
  192. C
  193.          DO 470 J = NA, N
  194.             T = MATA.RE(NA,J,1) + U2 * MATA.RE(EN,J,1)
  195.             MATA.RE(NA,J,1) = MATA.RE(NA,J,1) + T * V1
  196.             MATA.RE(EN,J,1) = MATA.RE(EN,J,1) + T * V2
  197.             T = MATB.RE(NA,J,1) + U2 * MATB.RE(EN,J,1)
  198.             MATB.RE(NA,J,1) = MATB.RE(NA,J,1) + T * V1
  199.             MATB.RE(EN,J,1) = MATB.RE(EN,J,1) + T * V2
  200.   470    CONTINUE
  201. C
  202.   475    continue
  203.          MATA.RE(EN,NA,1) = 0.0D0
  204.          ALFR.PROG(NA) = MATA.RE(NA,NA,1)
  205.          ALFR.PROG(EN) = MATA.RE(EN,EN,1)
  206.          MATB.RE(EN,NA,1) = 0.0D0
  207.          IF ( MATB.RE(NA,NA,1) . LT . 0.0D0 )
  208.      &     ALFR.PROG(NA) = -ALFR.PROG(NA)
  209.          IF ( MATB.RE(EN,EN,1) . LT . 0.0D0 )
  210.      &     ALFR.PROG(EN) = -ALFR.PROG(EN)
  211.          BETA.PROG(NA) = ABS(MATB.RE(NA,NA,1))
  212.          BETA.PROG(EN) = ABS(MATB.RE(EN,EN,1))
  213.          ALFI.PROG(EN) = 0.0D0
  214.          ALFI.PROG(NA) = 0.0D0
  215.          GOTO 505
  216. C
  217. C-----------------------------------------------------------------------
  218. C------------------- DEUX VALEURS PROPRES COMPLEXES --------------------
  219. C-----------------------------------------------------------------------
  220. C
  221.   480    E = E + C
  222.          EI = SQRT(-D)
  223.          A11R = A11 - E * B11
  224.          A11I = EI * B11
  225.          A12R = A12 - E * B12
  226.          A12I = EI * B12
  227.          A22R = A22 - E * B22
  228.          A22I = EI * B22
  229.          IF ( ABS(A11R) + ABS(A11I) + ABS(A12R) + ABS(A12I) . LT .
  230.      &        ABS(A21) + ABS(A22R) + ABS(A22I) ) GOTO 482
  231.          A1 = A12R
  232.          A1I = A12I
  233.          A2 = -A11R
  234.          A2I = -A11I
  235.          GOTO 485
  236.   482    A1 = A22R
  237.          A1I = A22I
  238.          A2 = -A21
  239.          A2I = 0.0D0
  240. C
  241. C-----------------------------------------------------------------------
  242. C----------------------- CHOISIR  Z  COMPLEXE --------------------------
  243. C-----------------------------------------------------------------------
  244. C
  245.   485    CZ = SQRT(A1*A1+A1I*A1I)
  246.          IF ( CZ . EQ . 0.0D0 ) GOTO 487
  247.          SZR = (A1 * A2 + A1I * A2I) / CZ
  248.          SZI = (A1 * A2I - A1I * A2) / CZ
  249.          R = SQRT(CZ*CZ+SZR*SZR+SZI*SZI)
  250.          CZ = CZ / R
  251.          SZR = SZR / R
  252.          SZI = SZI / R
  253.          GOTO 490
  254.   487    SZR = 1.0D0
  255.          SZI = 0.0D0
  256.   490    IF ( AN . LT . (ABS(E) + EI) * BN ) GOTO 492
  257.          A1 = CZ * B11 + SZR * B12
  258.          A1I = SZI * B12
  259.          A2 = SZR * B22
  260.          A2I = SZI * B22
  261.          GOTO 495
  262.   492    A1 = CZ * A11 + SZR * A12
  263.          A1I = SZI * A12
  264.          A2 = CZ * A21 + SZR * A22
  265.          A2I = SZI * A22
  266. C
  267. C-----------------------------------------------------------------------
  268. C----------------------- CHOISIR  Q  COMPLEXE --------------------------
  269. C-----------------------------------------------------------------------
  270. C
  271.   495    CQ = SQRT(A1*A1+A1I*A1I)
  272.          IF ( CQ . EQ . 0.0D0 ) GOTO 497
  273.          SQR = (A1 * A2 + A1I * A2I) / CQ
  274.          SQI = (A1 * A2I - A1I * A2) / CQ
  275.          R = SQRT(CQ*CQ+SQR*SQR+SQI*SQI)
  276.          CQ = CQ / R
  277.          SQR = SQR / R
  278.          SQI = SQI / R
  279.          GOTO 500
  280.   497    SQR = 1.0D0
  281.          SQI = 0.0D0
  282. C
  283. C-----------------------------------------------------------------------
  284. C------------------- CALCUL DES ELEMENTS DIAGONAUX  --------------------
  285. C-----------------------------------------------------------------------
  286. C
  287.   500    SSR = SQR * SZR + SQI * SZI
  288.          SSI = SQR * SZI - SQI * SZR
  289.          I = 1
  290.          TR = CQ * CZ * A11 + CQ * SZR * A12 + SQR * CZ * A21
  291.      &      + SSR * A22
  292.          TI = CQ * SZI * A12 - SQI * CZ * A21 + SSI * A22
  293.          DR = CQ * CZ * B11 + CQ * SZR * B12 + SSR * B22
  294.          DI = CQ * SZI * B12 + SSI * B22
  295.          GOTO 503
  296.   502    I = 2
  297.          TR = SSR * A11 - SQR * CZ * A12 - CQ * SZR * A21
  298.      &      + CQ * CZ * A22
  299.          TI = -SSI * A11 - SQI * CZ * A12 + CQ * SZI * A21
  300.          DR = SSR * B11 - SQR * CZ * B12 + CQ * CZ * B22
  301.          DI = -SSI * B11 - SQI * CZ * B12
  302.   503    T = TI * DR - TR * DI
  303.          J = NA
  304.          IF ( T . LT . 0.0D0 ) J = EN
  305.          R = SQRT(DR*DR+DI*DI)
  306.          BETA.PROG(J) = BN * R
  307.          ALFR.PROG(J) = AN * (TR * DR + TI * DI) / R
  308.          ALFI.PROG(J) = AN * T / R
  309.          IF ( I . EQ . 1 ) GOTO 502
  310.   505    ISW = 3 - ISW
  311.   510 CONTINUE
  312.       MATB.RE(N,1,1) = EPSB
  313. C
  314. c       SEGDES , MATA,MATB,MATZ
  315. c       SEGDES , ALFR,ALFI,BETA
  316. c
  317. C=======================================================================
  318. C========================== FIN DES CALCULS ============================
  319. C=======================================================================
  320.       RETURN
  321.       END
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  
  336.  
  337.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales