Télécharger qzvalp.eso

Retour à la liste

Numérotation des lignes :

  1. C QZVALP SOURCE CB215821 16/04/21 21:18:09 8920
  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 CCOPTIO
  41. -INC SMRIGID
  42. -INC SMLREEL
  43. C
  44. REAL*8 SWAP
  45. INTEGER I , J , EN , NA , NN , ISW
  46. C
  47. LOGICAL CALV
  48. C-----------------------------------------------------------------------
  49. C
  50. POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI
  51. cbp POINTEUR ALFR.XMATRI, ALFI.XMATRI, BETA.XMATRI
  52. POINTEUR ALFR.MLREEL, ALFI.MLREEL, BETA.MLREEL
  53. C
  54. c SEGACT , MATA*mod
  55. c SEGACT , MATB*mod
  56. c SEGACT , MATZ*mod
  57. NM=MATA.RE(/1)
  58. N=MATA.RE(/2)
  59. C
  60. JG = N
  61. SEGINI,ALFR,ALFI,BETA
  62. C-----------------------------------------------------------------------
  63. C========================= DEBUT DES CALCULS ===========================
  64. C-----------------------------------------------------------------------
  65. C
  66. EPSB = MATB.RE(N,1,1)
  67. ISW = 1
  68. C
  69. C-----------------------------------------------------------------------
  70. C-------------------------- BOUCLE PRINCIPALE --------------------------
  71. C-----------------------------------------------------------------------
  72. C
  73. DO 510 NN = 1, N
  74. EN = N + 1 - NN
  75. NA = EN - 1
  76. IF ( ISW . EQ . 2 ) GOTO 505
  77. IF ( EN . EQ . 1 ) GOTO 410
  78. SWAP = MATA.RE(EN,NA,1)
  79. IF ( SWAP . NE . 0.0D0 ) GOTO 420
  80. C
  81. C-----------------------------------------------------------------------
  82. C---------------------- UNE VALEUR PROPRE REELLE -----------------------
  83. C-----------------------------------------------------------------------
  84. C
  85. 410 continue
  86. ALFR.PROG(EN) = MATA.RE(EN,EN,1)
  87. IF ( MATB.RE(EN,EN,1) . LT . 0.0D0 )
  88. & ALFR.PROG(EN) = -ALFR.PROG(EN)
  89. BETA.PROG(EN) = ABS(MATB.RE(EN,EN,1))
  90. ALFI.PROG(EN) = 0.0D0
  91. GOTO 510
  92. C
  93. 420 continue
  94. SWAP = ABS(MATB.RE(NA,NA,1))
  95. IF ( SWAP . LE . EPSB ) GOTO 455
  96. SWAP = ABS(MATB.RE(EN,EN,1))
  97. IF ( SWAP . GT . EPSB ) GOTO 430
  98. A1 = MATA.RE(EN,EN,1)
  99. A2 = MATA.RE(EN,NA,1)
  100. BN = 0.0D0
  101. GOTO 435
  102. 430 continue
  103. AN = ABS(MATA.RE(NA,NA,1)) + ABS(MATA.RE(NA,EN,1))
  104. & + ABS(MATA.RE(EN,NA,1)) + ABS(MATA.RE(EN,EN,1))
  105. BN = ABS(MATB.RE(NA,NA,1)) + ABS(MATB.RE(NA,EN,1))
  106. & + ABS(MATB.RE(EN,EN,1))
  107. A11 = MATA.RE(NA,NA,1) / AN
  108. A12 = MATA.RE(NA,EN,1) / AN
  109. A21 = MATA.RE(EN,NA,1) / AN
  110. A22 = MATA.RE(EN,EN,1) / AN
  111. B11 = MATB.RE(NA,NA,1) / BN
  112. B12 = MATB.RE(NA,EN,1) / BN
  113. B22 = MATB.RE(EN,EN,1) / BN
  114. E = A11 / B11
  115. EI = A22 / B22
  116. S = A21 / (B11 * B22)
  117. T = (A22 - E * B22) / B22
  118. IF ( ABS(E) . LE . ABS(EI) ) GOTO 431
  119. E = EI
  120. T = (A11 - E * B11) / B11
  121. 431 C = 0.5D0 * (T - S * B12)
  122. D = C * C + S * (A12 - E * B12)
  123. IF ( D . LT . 0.0D0 ) GOTO 480
  124. C
  125. C-----------------------------------------------------------------------
  126. C-------------------- DEUX VALEURS PROPRES REELLES ---------------------
  127. C-----------------------------------------------------------------------
  128. C
  129. E = E + (C + DSIGN(SQRT(D),C))
  130. A11 = A11 - E * B11
  131. A12 = A12 - E * B12
  132. A22 = A22 - E * B22
  133. IF ( ABS(A11) + ABS(A12) . LT .
  134. & ABS(A21) + ABS(A22) ) GOTO 432
  135. A1 = A12
  136. A2 = A11
  137. GOTO 435
  138. 432 A1 = A22
  139. A2 = A21
  140. C
  141. C-----------------------------------------------------------------------
  142. C------------------------- CHOISIR Z REEL ----------------------------
  143. C-----------------------------------------------------------------------
  144. C
  145. 435 S = ABS(A1) + ABS(A2)
  146. U1 = A1 / S
  147. U2 = A2 / S
  148. R = DSIGN(SQRT(U1*U1+U2*U2),U1)
  149. V1 = -(U1 + R) / R
  150. V2 = -U2 / R
  151. U2 = V2 / V1
  152. C
  153. DO 440 I = 1, EN
  154. T = MATA.RE(I,EN,1) + U2 * MATA.RE(I,NA,1)
  155. MATA.RE(I,EN,1) = MATA.RE(I,EN,1) + T * V1
  156. MATA.RE(I,NA,1) = MATA.RE(I,NA,1) + T * V2
  157. T = MATB.RE(I,EN,1) + U2 * MATB.RE(I,NA,1)
  158. MATB.RE(I,EN,1) = MATB.RE(I,EN,1) + T * V1
  159. MATB.RE(I,NA,1) = MATB.RE(I,NA,1) + T * V2
  160. 440 CONTINUE
  161. C
  162. IF ( .NOT. CALV ) GOTO 450
  163. C
  164. DO 445 I = 1, N
  165. T = MATZ.RE(I,EN,1) + U2 * MATZ.RE(I,NA,1)
  166. MATZ.RE(I,EN,1) = MATZ.RE(I,EN,1) + T * V1
  167. MATZ.RE(I,NA,1) = MATZ.RE(I,NA,1) + T * V2
  168. 445 CONTINUE
  169. C
  170. 450 IF ( BN . EQ . 0.0D0 ) GOTO 475
  171. IF ( AN . LT . ABS(E) * BN ) GOTO 455
  172. A1 = MATB.RE(NA,NA,1)
  173. A2 = MATB.RE(EN,NA,1)
  174. GOTO 460
  175. 455 continue
  176. A1 = MATA.RE(NA,NA,1)
  177. A2 = MATA.RE(EN,NA,1)
  178. C
  179. C-----------------------------------------------------------------------
  180. C------------------------- CHOISIR Q REEL ----------------------------
  181. C-----------------------------------------------------------------------
  182. C
  183. 460 S = ABS(A1) + ABS(A2)
  184. IF ( S . EQ . 0.0D0 ) GOTO 475
  185. U1 = A1 / S
  186. U2 = A2 / S
  187. R = DSIGN(SQRT(U1*U1+U2*U2),U1)
  188. V1 = -(U1 + R) / R
  189. V2 = -U2 / R
  190. U2 = V2 / V1
  191. C
  192. DO 470 J = NA, N
  193. T = MATA.RE(NA,J,1) + U2 * MATA.RE(EN,J,1)
  194. MATA.RE(NA,J,1) = MATA.RE(NA,J,1) + T * V1
  195. MATA.RE(EN,J,1) = MATA.RE(EN,J,1) + T * V2
  196. T = MATB.RE(NA,J,1) + U2 * MATB.RE(EN,J,1)
  197. MATB.RE(NA,J,1) = MATB.RE(NA,J,1) + T * V1
  198. MATB.RE(EN,J,1) = MATB.RE(EN,J,1) + T * V2
  199. 470 CONTINUE
  200. C
  201. 475 continue
  202. MATA.RE(EN,NA,1) = 0.0D0
  203. ALFR.PROG(NA) = MATA.RE(NA,NA,1)
  204. ALFR.PROG(EN) = MATA.RE(EN,EN,1)
  205. MATB.RE(EN,NA,1) = 0.0D0
  206. IF ( MATB.RE(NA,NA,1) . LT . 0.0D0 )
  207. & ALFR.PROG(NA) = -ALFR.PROG(NA)
  208. IF ( MATB.RE(EN,EN,1) . LT . 0.0D0 )
  209. & ALFR.PROG(EN) = -ALFR.PROG(EN)
  210. BETA.PROG(NA) = ABS(MATB.RE(NA,NA,1))
  211. BETA.PROG(EN) = ABS(MATB.RE(EN,EN,1))
  212. ALFI.PROG(EN) = 0.0D0
  213. ALFI.PROG(NA) = 0.0D0
  214. GOTO 505
  215. C
  216. C-----------------------------------------------------------------------
  217. C------------------- DEUX VALEURS PROPRES COMPLEXES --------------------
  218. C-----------------------------------------------------------------------
  219. C
  220. 480 E = E + C
  221. EI = SQRT(-D)
  222. A11R = A11 - E * B11
  223. A11I = EI * B11
  224. A12R = A12 - E * B12
  225. A12I = EI * B12
  226. A22R = A22 - E * B22
  227. A22I = EI * B22
  228. IF ( ABS(A11R) + ABS(A11I) + ABS(A12R) + ABS(A12I) . LT .
  229. & ABS(A21) + ABS(A22R) + ABS(A22I) ) GOTO 482
  230. A1 = A12R
  231. A1I = A12I
  232. A2 = -A11R
  233. A2I = -A11I
  234. GOTO 485
  235. 482 A1 = A22R
  236. A1I = A22I
  237. A2 = -A21
  238. A2I = 0.0D0
  239. C
  240. C-----------------------------------------------------------------------
  241. C----------------------- CHOISIR Z COMPLEXE --------------------------
  242. C-----------------------------------------------------------------------
  243. C
  244. 485 CZ = SQRT(A1*A1+A1I*A1I)
  245. IF ( CZ . EQ . 0.0D0 ) GOTO 487
  246. SZR = (A1 * A2 + A1I * A2I) / CZ
  247. SZI = (A1 * A2I - A1I * A2) / CZ
  248. R = SQRT(CZ*CZ+SZR*SZR+SZI*SZI)
  249. CZ = CZ / R
  250. SZR = SZR / R
  251. SZI = SZI / R
  252. GOTO 490
  253. 487 SZR = 1.0D0
  254. SZI = 0.0D0
  255. 490 IF ( AN . LT . (ABS(E) + EI) * BN ) GOTO 492
  256. A1 = CZ * B11 + SZR * B12
  257. A1I = SZI * B12
  258. A2 = SZR * B22
  259. A2I = SZI * B22
  260. GOTO 495
  261. 492 A1 = CZ * A11 + SZR * A12
  262. A1I = SZI * A12
  263. A2 = CZ * A21 + SZR * A22
  264. A2I = SZI * A22
  265. C
  266. C-----------------------------------------------------------------------
  267. C----------------------- CHOISIR Q COMPLEXE --------------------------
  268. C-----------------------------------------------------------------------
  269. C
  270. 495 CQ = SQRT(A1*A1+A1I*A1I)
  271. IF ( CQ . EQ . 0.0D0 ) GOTO 497
  272. SQR = (A1 * A2 + A1I * A2I) / CQ
  273. SQI = (A1 * A2I - A1I * A2) / CQ
  274. R = SQRT(CQ*CQ+SQR*SQR+SQI*SQI)
  275. CQ = CQ / R
  276. SQR = SQR / R
  277. SQI = SQI / R
  278. GOTO 500
  279. 497 SQR = 1.0D0
  280. SQI = 0.0D0
  281. C
  282. C-----------------------------------------------------------------------
  283. C------------------- CALCUL DES ELEMENTS DIAGONAUX --------------------
  284. C-----------------------------------------------------------------------
  285. C
  286. 500 SSR = SQR * SZR + SQI * SZI
  287. SSI = SQR * SZI - SQI * SZR
  288. I = 1
  289. TR = CQ * CZ * A11 + CQ * SZR * A12 + SQR * CZ * A21
  290. & + SSR * A22
  291. TI = CQ * SZI * A12 - SQI * CZ * A21 + SSI * A22
  292. DR = CQ * CZ * B11 + CQ * SZR * B12 + SSR * B22
  293. DI = CQ * SZI * B12 + SSI * B22
  294. GOTO 503
  295. 502 I = 2
  296. TR = SSR * A11 - SQR * CZ * A12 - CQ * SZR * A21
  297. & + CQ * CZ * A22
  298. TI = -SSI * A11 - SQI * CZ * A12 + CQ * SZI * A21
  299. DR = SSR * B11 - SQR * CZ * B12 + CQ * CZ * B22
  300. DI = -SSI * B11 - SQI * CZ * B12
  301. 503 T = TI * DR - TR * DI
  302. J = NA
  303. IF ( T . LT . 0.0D0 ) J = EN
  304. R = SQRT(DR*DR+DI*DI)
  305. BETA.PROG(J) = BN * R
  306. ALFR.PROG(J) = AN * (TR * DR + TI * DI) / R
  307. ALFI.PROG(J) = AN * T / R
  308. IF ( I . EQ . 1 ) GOTO 502
  309. 505 ISW = 3 - ISW
  310. 510 CONTINUE
  311. MATB.RE(N,1,1) = EPSB
  312. C
  313. c SEGDES , MATA,MATB,MATZ
  314. c SEGDES , ALFR,ALFI,BETA
  315. c
  316. C=======================================================================
  317. C========================== FIN DES CALCULS ============================
  318. C=======================================================================
  319. RETURN
  320. END
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.  
  328.  
  329.  
  330.  
  331.  
  332.  
  333.  
  334.  
  335.  

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