Télécharger qzvecp.eso

Retour à la liste

Numérotation des lignes :

  1. C QZVECP SOURCE BP208322 16/01/21 21:15:21 8791
  2. SUBROUTINE QZVECP (MATA, MATB, ALFR, ALFI, BETA, MATZ)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8 (A-H,O-Z)
  5. *
  6. *************************************************************************
  7. * *
  8. * CALCUL DES VECTEURS PROPRES COMPLEXES *
  9. * ( QUATRIEME PAS DE LA METHODE "QZ" POUR LE CALCUL MODAL COMPLEXE ) *
  10. * _____________________________________________________________________ *
  11. * *
  12. * *
  13. * DATE : le 7 Avril 1995 *
  14. * AUTEURS : J.ANTUNES , Nicolas BENECH *
  15. * _____________________________________________________________________ *
  16. * *
  17. * MODULE(S) APPELANT(S) : VIBRAC *
  18. * *
  19. * MODULE(S) APPELE(S) : *
  20. * _____________________________________________________________________ *
  21. * *
  22. * EN ENTREE : *
  23. * - MATA : matrice reelle quasi triangulaire superieure (XMATRI) *
  24. * - MATB : matrice reelle triangulaire superieure (XMATRI) *
  25. * - ALFR, ALFI, BETA : liste des valeurs propres *
  26. * (parties reelle et imaginaire du numerateur + denominateur)*
  27. * - MATZ : produit des transformations effectuees auparavant *
  28. * *
  29. * EN SORTIE : *
  30. * - MATA : matrice identique a la matrice d'entree *
  31. * - MATB : matrice detruite *
  32. * - ALFR, ALFI, BETA : liste des valeurs propres *
  33. * (parties reelle et imaginaire du numerateur + denominateur)*
  34. * - MATZ : matrice des vecteurs propres, normalises de facon *
  35. * telle que le module de la composante la plus large *
  36. * soit unitaire. *
  37. * si ALFI = 0 : le vecteur est dans la colonne J de MATZ *
  38. * si ALFI <> 0 : la partie reelle est dans la colonne J, *
  39. * la partie imag. dans la colonne J+1. *
  40. * (ce vecteur propre correspond a la valeur propre *
  41. * de partie imagianire positive. Le vecteur propre *
  42. * conjugue correspond a la valeur propre conjuguee) *
  43. * *
  44. *************************************************************************
  45. *
  46. -INC CCOPTIO
  47. -INC SMRIGID
  48. -INC SMLREEL
  49. C
  50. REAL*8 SWAP
  51. INTEGER I , J , K , M , EN , II , JJ , NA , NN , ISW ,
  52. & ENM2
  53. C-----------------------------------------------------------------------
  54. POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI
  55. cbp POINTEUR ALFR.XMATRI, ALFI.XMATRI, BETA.XMATRI
  56. POINTEUR ALFR.MLREEL, ALFI.MLREEL, BETA.MLREEL
  57. C
  58. c SEGACT , MATA, ALFR,ALFI,BETA
  59. c SEGACT , MATB*mod
  60. c SEGACT , MATZ*mod
  61. NM=MATA.RE(/1)
  62. N=MATA.RE(/2)
  63. C
  64. C***********************************************************************
  65. C***********************************************************************
  66. C***********************************************************************
  67. C
  68. C=======================================================================
  69. C========================= DEBUT DES CALCULS ===========================
  70. C=======================================================================
  71. C
  72. EPSB = MATB.RE(N,1,1)
  73. ISW = 1
  74. C
  75. C-----------------------------------------------------------------------
  76. C-------------------------- BOUCLE PRINCIPALE --------------------------
  77. C-----------------------------------------------------------------------
  78. C
  79. DO 800 NN = 1, N
  80. EN = N + 1 - NN
  81. NA = EN - 1
  82. IF ( ISW . EQ . 2 ) GOTO 795
  83. SWAP = ALFI.PROG(EN)
  84. IF ( SWAP . NE . 0.0D0 ) GOTO 710
  85. C
  86. C-----------------------------------------------------------------------
  87. C------------------------- VECTEUR PROPRE REEL -------------------------
  88. C-----------------------------------------------------------------------
  89. C
  90. M = EN
  91. MATB.RE(EN,EN,1) = 1.0D0
  92. IF ( NA . EQ . 0 ) GOTO 800
  93. ALFM = ALFR.PROG(M)
  94. BETM = BETA.PROG(M)
  95. C
  96. DO 700 II = 1, NA
  97. I = EN - II
  98. W = BETM * MATA.RE(I,I,1) - ALFM * MATB.RE(I,I,1)
  99. R = 0.0D0
  100. C
  101. DO 610 J = M, EN
  102. 610 R = R + (BETM * MATA.RE(I,J,1)
  103. & - ALFM * MATB.RE(I,J,1)) * MATB.RE(J,EN,1)
  104. C
  105. IF ( I . EQ . 1 . OR . ISW . EQ . 2 ) GOTO 630
  106. IF ( BETM * MATA.RE(I,I-1,1) . EQ . 0.0D0 ) GOTO 630
  107. ZZ = W
  108. S = R
  109. GOTO 690
  110. 630 M = I
  111. IF ( ISW . EQ . 2 ) GOTO 640
  112. C
  113. T = W
  114. IF ( W . EQ . 0.0D0 ) T = EPSB
  115. MATB.RE(I,EN,1) = -R / T
  116. GOTO 700
  117. C
  118. 640 X = BETM * MATA.RE(I,I+1,1) - ALFM * MATB.RE(I,I+1,1)
  119. Y = BETM * MATA.RE(I+1,I,1)
  120. Q = W * ZZ - X * Y
  121. T = (X * S - ZZ * R) / Q
  122. MATB.RE(I,EN,1) = T
  123. IF ( ABS(X) . LE . ABS(ZZ) ) GOTO 650
  124. MATB.RE(I+1,EN,1) = (-R - W * T) / X
  125. GOTO 690
  126. 650 MATB.RE(I+1,EN,1) = (-S - Y * T) / ZZ
  127. 690 ISW = 3 - ISW
  128. 700 CONTINUE
  129. C
  130. GOTO 800
  131. C
  132. C-----------------------------------------------------------------------
  133. C----------------------- VECTEUR PROPRE COMPLEXE -----------------------
  134. C-----------------------------------------------------------------------
  135. C
  136. 710 M = NA
  137. ALMR = ALFR.PROG(M)
  138. ALMI = ALFI.PROG(M)
  139. BETM = BETA.PROG(M)
  140. C
  141. Y = BETM * MATA.RE(EN,NA,1)
  142. MATB.RE(NA,EN,1) = (ALMR * MATB.RE(EN,EN,1)
  143. & - BETM * MATA.RE(EN,EN,1)) / Y
  144. MATB.RE(NA,NA,1) = -ALMI * MATB.RE(EN,EN,1) / Y
  145. MATB.RE(EN,NA,1) = 0.0D0
  146. MATB.RE(EN,EN,1) = 1.0D0
  147. ENM2 = NA - 1
  148. IF ( ENM2 . EQ . 0 ) GOTO 795
  149. C
  150. DO 790 II = 1, ENM2
  151. I = NA - II
  152. RA = 0.0D0
  153. SA = 0.0D0
  154. C
  155. W = BETM * MATA.RE(I,I,1) - ALMR * MATB.RE(I,I,1)
  156. W1 = -ALMI * MATB.RE(I,I,1)
  157. C
  158. DO 760 J = M, EN
  159. X = BETM * MATA.RE(I,J,1) - ALMR * MATB.RE(I,J,1)
  160. X1 = -ALMI * MATB.RE(I,J,1)
  161. RA = RA + X * MATB.RE(J,NA,1) - X1 * MATB.RE(J,EN,1)
  162. SA = SA + X * MATB.RE(J,EN,1) + X1 * MATB.RE(J,NA,1)
  163. 760 CONTINUE
  164. C
  165. IF ( I . EQ . 1 . OR . ISW . EQ . 2 ) GOTO 770
  166. SWAP = MATA.RE(I,I-1)
  167. IF ( BETM * SWAP . EQ . 0.0D0 ) GOTO 770
  168. ZZ = W
  169. Z1 = W1
  170. R = RA
  171. S = SA
  172. ISW = 2
  173. GOTO 790
  174. 770 M = I
  175. IF ( ISW . EQ . 2 ) GOTO 780
  176. C
  177. TR = -RA
  178. TI = -SA
  179. 773 DR = W
  180. DI = W1
  181. C
  182. 775 IF ( ABS(DI) . GT . ABS(DR) ) GOTO 777
  183. RR = DI / DR
  184. D = DR + DI * RR
  185. T1 = (TR + TI * RR) / D
  186. T2 = (TI - TR * RR) / D
  187. GOTO (787,782), ISW
  188. 777 RR = DR / DI
  189. D = DR * RR + DI
  190. T1 = (TR * RR + TI) / D
  191. T2 = (TI * RR - TR) / D
  192. GOTO (787,782), ISW
  193. C
  194. 780 continue
  195. X = BETM * MATA.RE(I,I+1,1) - ALMR * MATB.RE(I,I+1,1)
  196. Y = BETM * MATA.RE(I+1,I,1)
  197. X1 = -ALMI * MATB.RE(I,I+1,1)
  198. TR = Y * RA - W * R + W1 * S
  199. TI = Y * SA - W * S - W1 * R
  200. DR = W * ZZ - W1 * Z1 - X * Y
  201. DI = W * Z1 + W1 * ZZ - X1 * Y
  202. IF ( DR . EQ . 0.0D0 . AND . DI . EQ . 0.0D0 ) DR = EPSB
  203. GOTO 775
  204. 782 continue
  205. MATB.RE(I+1,NA,1) = T1
  206. MATB.RE(I+1,EN,1) = T2
  207. ISW = 1
  208. IF ( ABS(Y) . GT . ABS(W) + ABS(W1) ) GOTO 785
  209. TR = -RA - X * MATB.RE(I+1,NA,1) + X1 * MATB.RE(I+1,EN,1)
  210. TI = -SA - X * MATB.RE(I+1,EN,1) - X1 * MATB.RE(I+1,NA,1)
  211. GOTO 773
  212. 785 continue
  213. T1 = (-R - ZZ * MATB.RE(I+1,NA,1)
  214. & + Z1 * MATB.RE(I+1,EN,1)) / Y
  215. T2 = (-S - ZZ * MATB.RE(I+1,EN,1)
  216. & - Z1 * MATB.RE(I+1,NA,1)) / Y
  217. 787 continue
  218. MATB.RE(I,NA,1) = T1
  219. MATB.RE(I,EN,1) = T2
  220. 790 CONTINUE
  221. C
  222. 795 ISW = 3 - ISW
  223. C
  224. C-----------------------------------------------------------------------
  225. C------------------- FIN DE LA SUBSTITUTION INVERSE --------------------
  226. C-----------------------------------------------------------------------
  227. C
  228. 800 CONTINUE
  229. C
  230. C-----------------------------------------------------------------------
  231. C------- ON TRANSFORME DANS LE SYSTEME DE COORDONNEES D ORIGINE --------
  232. C-----------------------------------------------------------------------
  233. C
  234. DO 880 JJ = 1, N
  235. J = N + 1 - JJ
  236. C
  237. DO 880 I = 1, N
  238. ZZ = 0.0D0
  239. C
  240. DO 860 K = 1, J
  241. 860 ZZ = ZZ + MATZ.RE(I,K,1) * MATB.RE(K,J,1)
  242. C
  243. MATZ.RE(I,J,1) = ZZ
  244. 880 CONTINUE
  245. C
  246. C-----------------------------------------------------------------------
  247. C----------------- NORMALISATION DES VECTEURS PROPRES ------------------
  248. C-----------------------------------------------------------------------
  249. C
  250. DO 950 J = 1, N
  251. D = 0.0D0
  252. IF ( ISW . EQ . 2 ) GOTO 920
  253. SWAP = ALFI.PROG(J)
  254. IF ( SWAP . NE . 0.0D0 ) GOTO 945
  255. C
  256. DO 890 I = 1, N
  257. IF ( ABS(MATZ.RE(I,J,1)) . GT . D )
  258. & D = ABS(MATZ.RE(I,J,1))
  259. 890 CONTINUE
  260. C
  261. DO 900 I = 1, N
  262. 900 MATZ.RE(I,J,1) = MATZ.RE(I,J,1) / D
  263. C
  264. GOTO 950
  265. C
  266. 920 continue
  267. DO 930 I = 1, N
  268. R = ABS(MATZ.RE(I,J-1,1)) + ABS(MATZ.RE(I,J,1))
  269. IF ( R . NE . 0.0D0 )
  270. & R = SQRT((MATZ.RE(I,J-1,1) )**2
  271. & + (MATZ.RE(I,J,1) )**2)
  272. IF ( R . GT . D ) D = R
  273. 930 CONTINUE
  274. C
  275. DO 940 I = 1, N
  276. MATZ.RE(I,J-1,1) = MATZ.RE(I,J-1,1) / D
  277. MATZ.RE(I,J,1) = MATZ.RE(I,J,1) / D
  278. 940 CONTINUE
  279. C
  280. 945 ISW = 3 - ISW
  281. 950 CONTINUE
  282. C
  283. c SEGDES , MATA,MATB,MATZ
  284. c SEGDES , ALFR,ALFI,BETA
  285. SEGSUP , MATA,MATB
  286. c
  287. C=======================================================================
  288. C========================== FIN DES CALCULS ============================
  289. C=======================================================================
  290. C
  291. RETURN
  292. END
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  

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