Télécharger qzvecp.eso

Retour à la liste

Numérotation des lignes :

qzvecp
  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.  
  47. -INC PPARAM
  48. -INC CCOPTIO
  49. -INC SMRIGID
  50. -INC SMLREEL
  51. C
  52. REAL*8 SWAP
  53. INTEGER I , J , K , M , EN , II , JJ , NA , NN , ISW ,
  54. & ENM2
  55. C-----------------------------------------------------------------------
  56. POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI
  57. cbp POINTEUR ALFR.XMATRI, ALFI.XMATRI, BETA.XMATRI
  58. POINTEUR ALFR.MLREEL, ALFI.MLREEL, BETA.MLREEL
  59. C
  60. c SEGACT , MATA, ALFR,ALFI,BETA
  61. c SEGACT , MATB*mod
  62. c SEGACT , MATZ*mod
  63. NM=MATA.RE(/1)
  64. N=MATA.RE(/2)
  65. C
  66. C***********************************************************************
  67. C***********************************************************************
  68. C***********************************************************************
  69. C
  70. C=======================================================================
  71. C========================= DEBUT DES CALCULS ===========================
  72. C=======================================================================
  73. C
  74. EPSB = MATB.RE(N,1,1)
  75. ISW = 1
  76. C
  77. C-----------------------------------------------------------------------
  78. C-------------------------- BOUCLE PRINCIPALE --------------------------
  79. C-----------------------------------------------------------------------
  80. C
  81. DO 800 NN = 1, N
  82. EN = N + 1 - NN
  83. NA = EN - 1
  84. IF ( ISW . EQ . 2 ) GOTO 795
  85. SWAP = ALFI.PROG(EN)
  86. IF ( SWAP . NE . 0.0D0 ) GOTO 710
  87. C
  88. C-----------------------------------------------------------------------
  89. C------------------------- VECTEUR PROPRE REEL -------------------------
  90. C-----------------------------------------------------------------------
  91. C
  92. M = EN
  93. MATB.RE(EN,EN,1) = 1.0D0
  94. IF ( NA . EQ . 0 ) GOTO 800
  95. ALFM = ALFR.PROG(M)
  96. BETM = BETA.PROG(M)
  97. C
  98. DO 700 II = 1, NA
  99. I = EN - II
  100. W = BETM * MATA.RE(I,I,1) - ALFM * MATB.RE(I,I,1)
  101. R = 0.0D0
  102. C
  103. DO 610 J = M, EN
  104. 610 R = R + (BETM * MATA.RE(I,J,1)
  105. & - ALFM * MATB.RE(I,J,1)) * MATB.RE(J,EN,1)
  106. C
  107. IF ( I . EQ . 1 . OR . ISW . EQ . 2 ) GOTO 630
  108. IF ( BETM * MATA.RE(I,I-1,1) . EQ . 0.0D0 ) GOTO 630
  109. ZZ = W
  110. S = R
  111. GOTO 690
  112. 630 M = I
  113. IF ( ISW . EQ . 2 ) GOTO 640
  114. C
  115. T = W
  116. IF ( W . EQ . 0.0D0 ) T = EPSB
  117. MATB.RE(I,EN,1) = -R / T
  118. GOTO 700
  119. C
  120. 640 X = BETM * MATA.RE(I,I+1,1) - ALFM * MATB.RE(I,I+1,1)
  121. Y = BETM * MATA.RE(I+1,I,1)
  122. Q = W * ZZ - X * Y
  123. T = (X * S - ZZ * R) / Q
  124. MATB.RE(I,EN,1) = T
  125. IF ( ABS(X) . LE . ABS(ZZ) ) GOTO 650
  126. MATB.RE(I+1,EN,1) = (-R - W * T) / X
  127. GOTO 690
  128. 650 MATB.RE(I+1,EN,1) = (-S - Y * T) / ZZ
  129. 690 ISW = 3 - ISW
  130. 700 CONTINUE
  131. C
  132. GOTO 800
  133. C
  134. C-----------------------------------------------------------------------
  135. C----------------------- VECTEUR PROPRE COMPLEXE -----------------------
  136. C-----------------------------------------------------------------------
  137. C
  138. 710 M = NA
  139. ALMR = ALFR.PROG(M)
  140. ALMI = ALFI.PROG(M)
  141. BETM = BETA.PROG(M)
  142. C
  143. Y = BETM * MATA.RE(EN,NA,1)
  144. MATB.RE(NA,EN,1) = (ALMR * MATB.RE(EN,EN,1)
  145. & - BETM * MATA.RE(EN,EN,1)) / Y
  146. MATB.RE(NA,NA,1) = -ALMI * MATB.RE(EN,EN,1) / Y
  147. MATB.RE(EN,NA,1) = 0.0D0
  148. MATB.RE(EN,EN,1) = 1.0D0
  149. ENM2 = NA - 1
  150. IF ( ENM2 . EQ . 0 ) GOTO 795
  151. C
  152. DO 790 II = 1, ENM2
  153. I = NA - II
  154. RA = 0.0D0
  155. SA = 0.0D0
  156. C
  157. W = BETM * MATA.RE(I,I,1) - ALMR * MATB.RE(I,I,1)
  158. W1 = -ALMI * MATB.RE(I,I,1)
  159. C
  160. DO 760 J = M, EN
  161. X = BETM * MATA.RE(I,J,1) - ALMR * MATB.RE(I,J,1)
  162. X1 = -ALMI * MATB.RE(I,J,1)
  163. RA = RA + X * MATB.RE(J,NA,1) - X1 * MATB.RE(J,EN,1)
  164. SA = SA + X * MATB.RE(J,EN,1) + X1 * MATB.RE(J,NA,1)
  165. 760 CONTINUE
  166. C
  167. IF ( I . EQ . 1 . OR . ISW . EQ . 2 ) GOTO 770
  168. SWAP = MATA.RE(I,I-1)
  169. IF ( BETM * SWAP . EQ . 0.0D0 ) GOTO 770
  170. ZZ = W
  171. Z1 = W1
  172. R = RA
  173. S = SA
  174. ISW = 2
  175. GOTO 790
  176. 770 M = I
  177. IF ( ISW . EQ . 2 ) GOTO 780
  178. C
  179. TR = -RA
  180. TI = -SA
  181. 773 DR = W
  182. DI = W1
  183. C
  184. 775 IF ( ABS(DI) . GT . ABS(DR) ) GOTO 777
  185. RR = DI / DR
  186. D = DR + DI * RR
  187. T1 = (TR + TI * RR) / D
  188. T2 = (TI - TR * RR) / D
  189. GOTO (787,782), ISW
  190. 777 RR = DR / DI
  191. D = DR * RR + DI
  192. T1 = (TR * RR + TI) / D
  193. T2 = (TI * RR - TR) / D
  194. GOTO (787,782), ISW
  195. C
  196. 780 continue
  197. X = BETM * MATA.RE(I,I+1,1) - ALMR * MATB.RE(I,I+1,1)
  198. Y = BETM * MATA.RE(I+1,I,1)
  199. X1 = -ALMI * MATB.RE(I,I+1,1)
  200. TR = Y * RA - W * R + W1 * S
  201. TI = Y * SA - W * S - W1 * R
  202. DR = W * ZZ - W1 * Z1 - X * Y
  203. DI = W * Z1 + W1 * ZZ - X1 * Y
  204. IF ( DR . EQ . 0.0D0 . AND . DI . EQ . 0.0D0 ) DR = EPSB
  205. GOTO 775
  206. 782 continue
  207. MATB.RE(I+1,NA,1) = T1
  208. MATB.RE(I+1,EN,1) = T2
  209. ISW = 1
  210. IF ( ABS(Y) . GT . ABS(W) + ABS(W1) ) GOTO 785
  211. TR = -RA - X * MATB.RE(I+1,NA,1) + X1 * MATB.RE(I+1,EN,1)
  212. TI = -SA - X * MATB.RE(I+1,EN,1) - X1 * MATB.RE(I+1,NA,1)
  213. GOTO 773
  214. 785 continue
  215. T1 = (-R - ZZ * MATB.RE(I+1,NA,1)
  216. & + Z1 * MATB.RE(I+1,EN,1)) / Y
  217. T2 = (-S - ZZ * MATB.RE(I+1,EN,1)
  218. & - Z1 * MATB.RE(I+1,NA,1)) / Y
  219. 787 continue
  220. MATB.RE(I,NA,1) = T1
  221. MATB.RE(I,EN,1) = T2
  222. 790 CONTINUE
  223. C
  224. 795 ISW = 3 - ISW
  225. C
  226. C-----------------------------------------------------------------------
  227. C------------------- FIN DE LA SUBSTITUTION INVERSE --------------------
  228. C-----------------------------------------------------------------------
  229. C
  230. 800 CONTINUE
  231. C
  232. C-----------------------------------------------------------------------
  233. C------- ON TRANSFORME DANS LE SYSTEME DE COORDONNEES D ORIGINE --------
  234. C-----------------------------------------------------------------------
  235. C
  236. DO 880 JJ = 1, N
  237. J = N + 1 - JJ
  238. C
  239. DO 880 I = 1, N
  240. ZZ = 0.0D0
  241. C
  242. DO 860 K = 1, J
  243. 860 ZZ = ZZ + MATZ.RE(I,K,1) * MATB.RE(K,J,1)
  244. C
  245. MATZ.RE(I,J,1) = ZZ
  246. 880 CONTINUE
  247. C
  248. C-----------------------------------------------------------------------
  249. C----------------- NORMALISATION DES VECTEURS PROPRES ------------------
  250. C-----------------------------------------------------------------------
  251. C
  252. DO 950 J = 1, N
  253. D = 0.0D0
  254. IF ( ISW . EQ . 2 ) GOTO 920
  255. SWAP = ALFI.PROG(J)
  256. IF ( SWAP . NE . 0.0D0 ) GOTO 945
  257. C
  258. DO 890 I = 1, N
  259. IF ( ABS(MATZ.RE(I,J,1)) . GT . D )
  260. & D = ABS(MATZ.RE(I,J,1))
  261. 890 CONTINUE
  262. C
  263. DO 900 I = 1, N
  264. 900 MATZ.RE(I,J,1) = MATZ.RE(I,J,1) / D
  265. C
  266. GOTO 950
  267. C
  268. 920 continue
  269. DO 930 I = 1, N
  270. R = ABS(MATZ.RE(I,J-1,1)) + ABS(MATZ.RE(I,J,1))
  271. IF ( R . NE . 0.0D0 )
  272. & R = SQRT((MATZ.RE(I,J-1,1) )**2
  273. & + (MATZ.RE(I,J,1) )**2)
  274. IF ( R . GT . D ) D = R
  275. 930 CONTINUE
  276. C
  277. DO 940 I = 1, N
  278. MATZ.RE(I,J-1,1) = MATZ.RE(I,J-1,1) / D
  279. MATZ.RE(I,J,1) = MATZ.RE(I,J,1) / D
  280. 940 CONTINUE
  281. C
  282. 945 ISW = 3 - ISW
  283. 950 CONTINUE
  284. C
  285. c SEGDES , MATA,MATB,MATZ
  286. c SEGDES , ALFR,ALFI,BETA
  287. SEGSUP , MATA,MATB
  288. c
  289. C=======================================================================
  290. C========================== FIN DES CALCULS ============================
  291. C=======================================================================
  292. C
  293. RETURN
  294. END
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  

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