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