Télécharger qzhess.eso

Retour à la liste

Numérotation des lignes :

  1. C QZHESS SOURCE BP208322 12/10/03 21:15:28 7517
  2. SUBROUTINE QZHESS (MATA, MATB, CALV, MATZ)
  3. c SUBROUTINE QZHESS (XMATA, XMATB, CALV, XMATZ, NM,N)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. *
  7. *************************************************************************
  8. * *
  9. * MISE SOUS FORME DE HESSENBERG SUPERIEURE DE MATA *
  10. * ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE MATB *
  11. * ( PREMIER PAS DE LA METHODE "QZ" POUR LE CALCUL MODAL COMPLEXE ) *
  12. * _____________________________________________________________________ *
  13. * *
  14. * *
  15. * DATE : le 23 Mars 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 generale (XMATRI) *
  26. * - MATB : matrice reelle generale (XMATRI) *
  27. * - CALV : indique si l'on calcule ou non les vecteurs propres *
  28. * *
  29. * EN SORTIE : *
  30. * - MATA : matrice sous forme Hessenberg superieure *
  31. * - MATB : matrice sous forme triangulaire superieure *
  32. * - MATZ : produit des transformations effectuees (si CALV=.TRUE.)*
  33. * *
  34. *************************************************************************
  35. *
  36. -INC CCOPTIO
  37. -INC SMRIGID
  38. *
  39. INTEGER I , J , K , L , N , LB , L1 , NM , NK1 , NM1 , NM2
  40. LOGICAL CALV
  41. *-----------------------------------------------------------------------
  42. POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI
  43. * REAL*8 XMATA(NM,N),XMATB(NM,N),XMATZ(NM,N)
  44. *
  45. c SEGACT , MATA*mod
  46. c SEGACT , MATB*mod
  47. NM=MATA.RE(/1)
  48. N=MATA.RE(/2)
  49. *
  50. *=======================================================================
  51. *========================= DEBUT DES CALCULS ===========================
  52. *=======================================================================
  53. *
  54. *-----------------------------------------------------------------------
  55. *----------------------- INITIALISATION DE Z -------------------------
  56. *-----------------------------------------------------------------------
  57. *
  58. NLIGRD=NM
  59. NLIGRP=N
  60. nelrig=1
  61. SEGINI , MATZ
  62. *
  63. IF ( .NOT. CALV ) GOTO 9
  64. DO 3 J = 1, N
  65. DO 2 I = 1, N
  66. MATZ.RE(I,J,1) = 0.0D0
  67. 2 CONTINUE
  68. MATZ.RE(J,J,1) = 1.0D0
  69. 3 CONTINUE
  70. *
  71. 9 continue
  72. *
  73. *-----------------------------------------------------------------------
  74. *--------- REDUCTION DE B A LA FORME TRIANGULAIRE SUPERIEURE ---------
  75. *-----------------------------------------------------------------------
  76. *
  77. 10 IF ( N .LE. 1 ) GOTO 170
  78. *
  79. *
  80. NM1 = N - 1
  81. *
  82. DO 100 L = 1, NM1
  83. L1 = L + 1
  84.  
  85. c S = somme des termes de B_il i=l+1..n
  86. S = 0.0D0
  87. DO 20 I = L1, N
  88. S = S + ABS(MATB.RE(I,L,1))
  89. 20 CONTINUE
  90. c S=0 : la ligne L verifie deja la forme triangulaire infer de B
  91. IF ( S .EQ. 0.0D0 ) GOTO 100
  92. c S = somme des termes de B_il i=l..n : ajout du terme diagonal
  93. S = S + ABS(MATB.RE(L,L,1))
  94.  
  95. c methode de Householder pour annuler les termes non-nul
  96. c
  97. c calcul de r = sign(B_ll)*sqrt(sum B_il^2)
  98. c (rem: r est souvent appelé \alpha dans la litterature)
  99. R = 0.0D0
  100. DO 25 I = L, N
  101. MATB.RE(I,L,1) = MATB.RE(I,L,1) / S
  102. R = R + MATB.RE(I,L,1)**2
  103. 25 CONTINUE
  104. R = SIGN(SQRT(R),MATB.RE(L,L,1))
  105. MATB.RE(L,L,1) = MATB.RE(L,L,1) + R
  106. RHO = R * MATB.RE(L,L,1)
  107. c triangularisation effective de B
  108. DO 50 J = L1, N
  109. T = 0.0D0
  110. DO 30 I = L, N
  111. T = T + MATB.RE(I,L,1) * MATB.RE(I,J,1)
  112. 30 CONTINUE
  113. T = -T / RHO
  114. DO 40 I = L, N
  115. MATB.RE(I,J,1) = MATB.RE(I,J,1) + T * MATB.RE(I,L,1)
  116. 40 CONTINUE
  117. 50 CONTINUE
  118. c la meme transformation est appliquee a A
  119. DO 80 J = 1, N
  120. T = 0.0D0
  121. DO 60 I = L, N
  122. T = T + MATB.RE(I,L,1) * MATA.RE(I,J,1)
  123. 60 CONTINUE
  124. T = -T / RHO
  125. DO 70 I = L, N
  126. MATA.RE(I,J,1) = MATA.RE(I,J,1) + T * MATB.RE(I,L,1)
  127. 70 CONTINUE
  128. *
  129. 80 CONTINUE
  130. c diagonale de B
  131. MATB.RE(L,L,1) = -S * R
  132. c mise a 0 des termes post-diagonaux
  133. DO 90 I = L1, N
  134. MATB.RE(I,L,1) = 0.0D0
  135. 90 CONTINUE
  136. *
  137. 100 CONTINUE
  138. *
  139. *
  140. C-----------------------------------------------------------------------
  141. C-------- REDUCTION DE A A LA FORME DE HESSENBERG SUPERIEURE ---------
  142. C-----------------------------------------------------------------------
  143. c while keeping B triangular...
  144. *
  145. IF ( N .EQ. 2 ) GOTO 170
  146. C
  147. NM2 = N - 2
  148. *
  149. DO 160 K = 1, NM2
  150. NK1 = NM1 - K
  151. *
  152. DO 150 LB = 1, NK1
  153. L = N - LB
  154. L1 = L + 1
  155. *
  156. S = ABS(MATA.RE(L,K,1)) + ABS(MATA.RE(L1,K,1))
  157. IF ( S .EQ. 0.0D0 ) GOTO 150
  158. U1 = MATA.RE(L,K,1) / S
  159. U2 = MATA.RE(L1,K,1) / S
  160. R = SIGN(SQRT(U1*U1+U2*U2),U1)
  161. V1 = -(U1 + R) / R
  162. V2 = -U2 / R
  163. U2 = V2 / V1
  164. *
  165. DO 110 J = K, N
  166. T = MATA.RE(L,J,1) + U2 * MATA.RE(L1,J,1)
  167. MATA.RE(L,J,1) = MATA.RE(L,J,1) + T * V1
  168. MATA.RE(L1,J,1) = MATA.RE(L1,J,1) + T * V2
  169. 110 CONTINUE
  170. MATA.RE(L1,K,1) = 0.0D0
  171. *
  172. DO 120 J = L, N
  173. T = MATB.RE(L,J,1) + U2 * MATB.RE(L1,J,1)
  174. MATB.RE(L,J,1) = MATB.RE(L,J,1) + T * V1
  175. MATB.RE(L1,J,1) = MATB.RE(L1,J,1) + T * V2
  176. 120 CONTINUE
  177. *
  178. S = ABS(MATB.RE(L1,L1,1)) + ABS(MATB.RE(L1,L,1))
  179. IF ( S .EQ. 0.0D0 ) GOTO 150
  180. U1 = MATB.RE(L1,L1,1) / S
  181. U2 = MATB.RE(L1,L,1) / S
  182. R = SIGN(SQRT(U1*U1+U2*U2),U1)
  183. V1 = -(U1 + R) / R
  184. V2 = -U2 / R
  185. U2 = V2 / V1
  186. *
  187. DO 130 I = 1, L1
  188. T = MATB.RE(I,L1,1) + U2 * MATB.RE(I,L,1)
  189. MATB.RE(I,L1,1) = MATB.RE(I,L1,1) + T * V1
  190. MATB.RE(I,L,1) = MATB.RE(I,L,1) + T * V2
  191. 130 CONTINUE
  192. MATB.RE(L1,L,1) = 0.0D0
  193. *
  194. DO 140 I = 1, N
  195. T = MATA.RE(I,L1,1) + U2 * MATA.RE(I,L,1)
  196. MATA.RE(I,L1,1) = MATA.RE(I,L1,1) + T * V1
  197. MATA.RE(I,L,1) = MATA.RE(I,L,1) + T * V2
  198. 140 CONTINUE
  199. *
  200. IF ( .NOT. CALV ) GOTO 150
  201. *
  202. DO 145 I = 1, N
  203. T = MATZ.RE(I,L1,1) + U2 * MATZ.RE(I,L,1)
  204. MATZ.RE(I,L1,1) = MATZ.RE(I,L1,1) + T * V1
  205. MATZ.RE(I,L,1) = MATZ.RE(I,L,1) + T * V2
  206. 145 CONTINUE
  207. *
  208. 150 CONTINUE
  209. *
  210. 160 CONTINUE
  211. *
  212. 170 CONTINUE
  213. *
  214. *========================== FIN DES CALCULS ============================
  215. c SEGDES , MATA,MATB,MATZ
  216. *
  217. RETURN
  218. END
  219.  
  220.  
  221.  
  222.  
  223.  
  224.  
  225.  
  226.  
  227.  
  228.  
  229.  
  230.  

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