1 : C QZVALP SOURCE BP208322 12/10/03 21:15:31 7517 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 : C 43 : REAL*8 SWAP 44 : INTEGER I , J , EN , NA , NN , ISW 45 : C 46 : LOGICAL CALV 47 : C----------------------------------------------------------------------- 48 : C 49 : POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI 50 : POINTEUR ALFR.XMATRI, ALFI.XMATRI, BETA.XMATRI 51 : C 52 : c SEGACT , MATA*mod 53 : c SEGACT , MATB*mod 54 : c SEGACT , MATZ*mod 55 : NM=MATA.RE(/1) 56 : N=MATA.RE(/2) 57 : C 58 : NLIGRD = 1 59 : NLIGRP = N 60 : nelrig=1 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.RE(1,EN,1) = MATA.RE(EN,EN,1) 87 : IF ( MATB.RE(EN,EN,1) . LT . 0.0D0 ) 88 : & ALFR.RE(1,EN,1) = -ALFR.RE(1,EN,1) 89 : BETA.RE(1,EN,1) = DABS(MATB.RE(EN,EN,1)) 90 : ALFI.RE(1,EN,1) = 0.0D0 91 : GOTO 510 92 : C 93 : 420 continue 94 : SWAP = DABS(MATB.RE(NA,NA,1)) 95 : IF ( SWAP . LE . EPSB ) GOTO 455 96 : SWAP = DABS(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 = DABS(MATA.RE(NA,NA,1)) + DABS(MATA.RE(NA,EN,1)) 104 : & + DABS(MATA.RE(EN,NA,1)) + DABS(MATA.RE(EN,EN,1)) 105 : BN = DABS(MATB.RE(NA,NA,1)) + DABS(MATB.RE(NA,EN,1)) 106 : & + DABS(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 ( DABS(E) . LE . DABS(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(DSQRT(D),C)) 130 : A11 = A11 - E * B11 131 : A12 = A12 - E * B12 132 : A22 = A22 - E * B22 133 : IF ( DABS(A11) + DABS(A12) . LT . 134 : & DABS(A21) + DABS(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 = DABS(A1) + DABS(A2) 146 : U1 = A1 / S 147 : U2 = A2 / S 148 : R = DSIGN(DSQRT(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 . DABS(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 = DABS(A1) + DABS(A2) 184 : IF ( S . EQ . 0.0D0 ) GOTO 475 185 : U1 = A1 / S 186 : U2 = A2 / S 187 : R = DSIGN(DSQRT(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.RE(1,NA,1) = MATA.RE(NA,NA,1) 204 : ALFR.RE(1,EN,1) = 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.RE(1,NA,1) = -ALFR.RE(1,NA,1) 208 : IF ( MATB.RE(EN,EN,1) . LT . 0.0D0 ) 209 : & ALFR.RE(1,EN,1) = -ALFR.RE(1,EN,1) 210 : BETA.RE(1,NA,1) = DABS(MATB.RE(NA,NA,1)) 211 : BETA.RE(1,EN,1) = DABS(MATB.RE(EN,EN,1)) 212 : ALFI.RE(1,EN,1) = 0.0D0 213 : ALFI.RE(1,NA,1) = 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 = DSQRT(-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 ( DABS(A11R) + DABS(A11I) + DABS(A12R) + DABS(A12I) . LT . 229 : & DABS(A21) + DABS(A22R) + DABS(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 = DSQRT(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 = DSQRT(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 . (DABS(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 = DSQRT(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 = DSQRT(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 = DSQRT(DR*DR+DI*DI) 305 : BETA.RE(1,J,1) = BN * R 306 : ALFR.RE(1,J,1) = AN * (TR * DR + TI * DI) / R 307 : ALFI.RE(1,J,1) = 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 :
© Cast3M 2003 - All rights reserved.
Disclaimer