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 :
© Cast3M 2003 - All rights reserved.
Disclaimer