1 : C QZREDU SOURCE BP208322 12/10/03 21:15:29 7517 2 : SUBROUTINE QZREDU (MATA, MATB, EPS1, CALV, MATZ, CERR) 3 : c SUBROUTINE QZREDU (XMATA, XMATB, EPS1, CALV, XMATZ,NM,N, CERR) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : * 7 : ************************************************************************* 8 : * * 9 : * MISE SOUS FORME QUASI-TRIANGULAIRE SUPERIEURE * 10 : * D'UNE MATRICE DE HESSENBERG SUPERIEURE (MATA) * 11 : * ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE MATB * 12 : * ( DEUXIEME PAS DE LA METHODE "QZ" POUR LE CALCUL MODAL COMPLEXE ) * 13 : * _____________________________________________________________________ * 14 : * * 15 : * * 16 : * DATE : le 27 Mars 1995 * 17 : * AUTEURS : J.ANTUNES , Nicolas BENECH * 18 : * _____________________________________________________________________ * 19 : * * 20 : * MODULE(S) APPELANT(S) : VIBRAC * 21 : * * 22 : * MODULE(S) APPELE(S) : * 23 : * _____________________________________________________________________ * 24 : * * 25 : * EN ENTREE : * 26 : * - MATA : matrice reelle Hessenberg superieure (XMATRI) * 27 : * - MATB : matrice reelle triangulaire superieure (XMATRI) * 28 : * - EPS1 : tolerance pour les elements negligeables * 29 : * si EPS1 = 0. ou EPS1 < 0., un element est neglige * 30 : * s'il est inferieur a l'erreur d'arrondi multiplie * 31 : * par la norme de la matrice * 32 : * si EPS1 > 0., un element est neglige s'il est inferieur * 33 : * a EPS1 multiplie par la norme de la matrice * 34 : * (plus rapide mais moins fiable) * 35 : * - CALV : indique si l'on calcule ou non les vecteurs propres * 36 : * * 37 : * EN SORTIE : * 38 : * - MATA : matrice quasi-triangulaire superieure * 39 : * - MATB : matrice sous forme triangulaire superieure * 40 : * - MATZ : produit des transformations effectuees (si CALV=.TRUE.)* 41 : * - CERR : code des erreurs * 42 : * CERR = 0 : retour normal * 43 : * CERR = J : plus de 30*N iterations pour le calcul * 44 : * de la valeur propre J * 45 : * * 46 : ************************************************************************* 47 : * 48 : -INC CCOPTIO 49 : -INC SMRIGID 50 : C 51 : INTEGER I , J , K , L , EN , K1 , K2 , LD , LL , L1 , NA , 52 : & ISH , ITN , ITS , KM1 , LM1 , ENM2 , CERR , LOR1 , 53 : & ENORN 54 : C 55 : LOGICAL CALV , NOTLAS 56 : C----------------------------------------------------------------------- 57 : POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI 58 : c REAL*8 XMATA(NM,N),XMATB(NM,N),XMATZ(NM,N) 59 : C 60 : c SEGACT , MATA*mod 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 : if(iimpi.ge.6) then 68 : write(ioimp,*) 'A=' 69 : do iou=1,N 70 : write(ioimp,*) (MATA.RE(iou,jou,1),jou=1,N) 71 : enddo 72 : write(ioimp,*) 'B=' 73 : do iou=1,N 74 : write(ioimp,*) (MATB.RE(iou,jou,1),jou=1,N) 75 : enddo 76 : endif 77 : C*********************************************************************** 78 : C*********************************************************************** 79 : C*********************************************************************** 80 : C 81 : C======================================================================= 82 : C========================= DEBUT DES CALCULS =========================== 83 : C======================================================================= 84 : C 85 : CERR = 0 86 : C 87 : C----------------------------------------------------------------------- 88 : C---------------------- CALCUL DE EPSA ET EPSB ---------------------- 89 : C----------------------------------------------------------------------- 90 : C 91 : ANORM = 0.0D0 92 : BNORM = 0.0D0 93 : C 94 : DO 30 I = 1, N 95 : ANI = 0.0D0 96 : IF ( I .NE. 1 ) ANI = ABS(MATA.RE(I,I-1,1)) 97 : BNI = 0.0D0 98 : C 99 : DO 20 J = I, N 100 : ANI = ANI + ABS(MATA.RE(I,J,1)) 101 : BNI = BNI + ABS(MATB.RE(I,J,1)) 102 : 20 CONTINUE 103 : C 104 : IF ( ANI .GT. ANORM ) ANORM = ANI 105 : IF ( BNI .GT. BNORM ) BNORM = BNI 106 : 30 CONTINUE 107 : C 108 : IF ( ANORM .EQ. 0.0D0 ) ANORM = 1.0D0 109 : IF ( BNORM .EQ. 0.0D0 ) BNORM = 1.0D0 110 : EP = EPS1 111 : IF ( EP .GT. 0.0D0 ) GOTO 50 112 : C 113 : EP = EPSLON(1.0D0) 114 : IF (EP .GT. 0.0D0) GO TO 50 115 : C ********** COMPUTE ROUNDOFF LEVEL IF EPSLON(1.) IS ZERO ********** 116 : EP = 1.0D0 117 : 40 EP = EP / 2.0 118 : IF (1.0D0 + EP .GT. 1.0D0) GO TO 40 119 : 50 EPSA = EP * ANORM 120 : EPSB = EP * BNORM 121 : C ********** REDUCE A TO QUASI-TRIANGULAR FORM, WHILE 122 : C KEEPING B TRIANGULAR ********** 123 : LOR1 = 1 124 : ENORN = N 125 : EN = N 126 : c ITN = 30*N 127 : c bp 2012-10-02 : on augmente le nombre d'iterations max 128 : ITN = 60*N 129 : C 130 : C----------------------------------------------------------------------- 131 : C------------------------- DEBUT DU PAS QZ --------------------------- 132 : C----------------------------------------------------------------------- 133 : C 134 : 60 IF ( EN .LE. 2 ) GOTO 1001 135 : IF ( .NOT. CALV ) ENORN = EN 136 : ITS = 0 137 : NA = EN - 1 138 : ENM2 = NA - 1 139 : 70 ISH = 2 140 : C 141 : C----------------------------------------------------------------------- 142 : C------------------- VERIFICATION DE LA CONVERGENCE -------------------- 143 : C----------------------------------------------------------------------- 144 : C 145 : DO 80 LL = 1, EN 146 : LM1 = EN - LL 147 : L = LM1 + 1 148 : IF ( L .EQ. 1 ) GOTO 95 149 : IF ( ABS(MATA.RE(L,LM1,1)) .LE. EPSA ) GOTO 90 150 : 80 CONTINUE 151 : C 152 : 90 MATA.RE(L,LM1,1) = 0.0D0 153 : IF ( L .LT. NA ) GOTO 95 154 : C ********** 1-BY-1 OR 2-BY-2 BLOCK ISOLATED ********** 155 : EN = LM1 156 : GOTO 60 157 : C ********** CHECK FOR SMALL TOP OF B ********** 158 : 95 LD = L 159 : 100 L1 = L + 1 160 : C 161 : B11 = MATB.RE(L,L,1) 162 : IF ( ABS(B11) .GT. EPSB ) GOTO 120 163 : MATB.RE(L,L,1) = 0.0D0 164 : C 165 : S = ABS(MATA.RE(L,L,1)) + ABS(MATA.RE(L1,L,1)) 166 : U1 = MATA.RE(L,L,1) / S 167 : U2 = MATA.RE(L1,L,1) / S 168 : C 169 : C 170 : R = SIGN(SQRT(U1*U1+U2*U2),U1) 171 : V1 = -(U1 + R) / R 172 : V2 = -U2 / R 173 : U2 = V2 / V1 174 : C 175 : DO 110 J = L, ENORN 176 : T = MATA.RE(L,J,1) + U2 * MATA.RE(L1,J,1) 177 : MATA.RE(L,J,1) = MATA.RE(L,J,1) + T * V1 178 : MATA.RE(L1,J,1) = MATA.RE(L1,J,1) + T * V2 179 : T = MATB.RE(L,J,1) + U2 * MATB.RE(L1,J,1) 180 : MATB.RE(L,J,1) = MATB.RE(L,J,1) + T * V1 181 : MATB.RE(L1,J,1) = MATB.RE(L1,J,1) + T * V2 182 : 110 CONTINUE 183 : C 184 : IF ( L .NE. 1 ) MATA.RE(L,LM1,1) = -MATA.RE(L,LM1,1) 185 : LM1 = L 186 : L = L1 187 : GOTO 90 188 : 120 CONTINUE 189 : A11 = MATA.RE(L,L,1) / B11 190 : A21 = MATA.RE(L1,L,1) / B11 191 : IF ( ISH .EQ. 1 ) GOTO 141 192 : C 193 : C----------------------------------------------------------------------- 194 : C---------------------- STRATEGIE DE L ITERATION ----------------------- 195 : C----------------------------------------------------------------------- 196 : C 197 : C ********** ITERATION STRATEGY ********** 198 : IF ( ITN .EQ. 0 ) GOTO 1000 199 : c bp 2012-10-02 : on change le moment d'appel au "ad hoc shift" dont on 200 : c ne trouve pas de justification mathematique mais on la cale 201 : c numeriquement sur un cas test problematique : brasero_tube3D.dgibi 202 : c l'ideal serait d'avoir un critere qui permet de declencher de shift 203 : c (ou bien de passer a linpack et DHGEQZ) 204 : c IF ( ITS .EQ. 10 ) GOTO 155 205 : IF ( ITS .EQ. 10*N) GOTO 155 206 : 207 : C 208 : C ********** DETERMINE TYPE OF SHIFT ********** 209 : B22 = MATB.RE(L1,L1,1) 210 : IF ( ABS(B22) .LT. EPSB ) B22 = EPSB 211 : B33 = MATB.RE(NA,NA,1) 212 : IF ( ABS(B33) .LT. EPSB ) B33 = EPSB 213 : B44 = MATB.RE(EN,EN,1) 214 : IF ( ABS(B44) .LT. EPSB ) B44 = EPSB 215 : A33 = MATA.RE(NA,NA,1) / B33 216 : A34 = MATA.RE(NA,EN,1) / B44 217 : A43 = MATA.RE(EN,NA,1) / B33 218 : A44 = MATA.RE(EN,EN,1) / B44 219 : B34 = MATB.RE(NA,EN,1) / B44 220 : T = 0.5D0 * (A43 * B34 - A33 - A44) 221 : R = T * T + A34 * A43 - A33 * A44 222 : IF ( R .LT. 0.0D0 ) GOTO 150 223 : C 224 : C ********** DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A ********** 225 : ISH = 1 226 : R = SQRT(R) 227 : SH = -T + R 228 : S = -T - R 229 : IF ( ABS(S-A44) .LT. ABS(SH-A44) ) SH = S 230 : C 231 : C ********** LOOK FOR TWO CONSECUTIVE SMALL 232 : C SUB-DIAGONAL ELEMENTS OF A. 233 : C FOR L=EN-2 STEP -1 UNTIL LD DO -- ********** 234 : DO 130 LL = LD, ENM2 235 : L = ENM2 + LD - LL 236 : IF ( L .EQ. LD ) GOTO 140 237 : LM1 = L - 1 238 : L1 = L + 1 239 : T = MATA.RE(L,L,1) 240 : IF ( ABS(MATB.RE(L,L,1)) .GT. EPSB ) 241 : & T = T - SH * MATB.RE(L,L,1) 242 : IF (ABS(MATA.RE(L,LM1,1)).LE.ABS(T/MATA.RE(L1,L,1))*EPSA) 243 : & GOTO 100 244 : 130 CONTINUE 245 : 140 continue 246 : C 247 : 141 A1 = A11 - SH 248 : A2 = A21 249 : IF ( L .NE. LD ) MATA.RE(L,LM1,1) = -MATA.RE(L,LM1,1) 250 : GOTO 160 251 : C 252 : C ********** DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A ********** 253 : 150 continue 254 : A12 = MATA.RE(L,L1,1) / B22 255 : A22 = MATA.RE(L1,L1,1) / B22 256 : B12 = MATB.RE(L,L1,1) / B22 257 : A1 = ((A33 - A11) * (A44 - A11) - A34 * A43 + A43 * B34 * A11) 258 : & / A21 + A12 - A11 * B12 259 : A2 = (A22 - A11) - A21 * B12 - (A33 - A11) - (A44 - A11) 260 : & + A43 * B34 261 : A3 = MATA.RE(L1+1,L1,1) / B22 262 : GOTO 160 263 : C 264 : C ********** AD HOC SHIFT ********** 265 : 155 A1 = 0.0D0 266 : A2 = 1.0D0 267 : A3 = 1.1605D0 268 : 160 ITS = ITS + 1 269 : ITN = ITN - 1 270 : if(iimpi.ge.6) write(*,*)'iteration ',ITS,' A1,A2,A3=',A1,A2,A3 271 : IF ( .NOT. CALV ) LOR1 = LD 272 : C 273 : C----------------------------------------------------------------------- 274 : C------------------------- BOUCLE PRINCIPALE --------------------------- 275 : C----------------------------------------------------------------------- 276 : C 277 : DO 260 K = L, NA 278 : NOTLAS = K .NE. NA .AND. ISH .EQ. 2 279 : K1 = K + 1 280 : K2 = K + 2 281 : KM1 = MAX0(K-1,L) 282 : LL = MIN0(EN,K1+ISH) 283 : IF ( NOTLAS ) GOTO 190 284 : C 285 : C ********** ZERO A(K+1,K-1) ********** 286 : IF ( K .EQ. L ) GOTO 170 287 : A1 = MATA.RE(K,KM1,1) 288 : A2 = MATA.RE(K1,KM1,1) 289 : 170 S = ABS(A1) + ABS(A2) 290 : IF ( S .EQ. 0.0D0 ) GOTO 70 291 : U1 = A1 / S 292 : U2 = A2 / S 293 : R = SIGN(SQRT(U1*U1+U2*U2),U1) 294 : V1 = -(U1 + R) / R 295 : V2 = -U2 / R 296 : U2 = V2 / V1 297 : C 298 : DO 180 J = KM1, ENORN 299 : T = MATA.RE(K,J,1) + U2 * MATA.RE(K1,J,1) 300 : MATA.RE(K,J,1) = MATA.RE(K,J,1) + T * V1 301 : MATA.RE(K1,J,1) = MATA.RE(K1,J,1) + T * V2 302 : T = MATB.RE(K,J,1) + U2 * MATB.RE(K1,J,1) 303 : MATB.RE(K,J,1) = MATB.RE(K,J,1) + T * V1 304 : MATB.RE(K1,J,1) = MATB.RE(K1,J,1) + T * V2 305 : 180 CONTINUE 306 : C 307 : IF ( K .NE. L ) MATA.RE(K1,KM1,1) = 0.0D0 308 : GOTO 240 309 : C 310 : C ********** ZERO A(K+1,K-1) AND A(K+2,K-1) ********** 311 : 190 IF ( K .EQ. L ) GOTO 200 312 : A1 = MATA.RE(K,KM1,1) 313 : A2 = MATA.RE(K1,KM1,1) 314 : A3 = MATA.RE(K2,KM1,1) 315 : 200 S = ABS(A1) + ABS(A2) + ABS(A3) 316 : IF ( S .EQ. 0.0D0 ) GOTO 260 317 : U1 = A1 / S 318 : U2 = A2 / S 319 : U3 = A3 / S 320 : R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1) 321 : V1 = -(U1 + R) / R 322 : V2 = -U2 / R 323 : V3 = -U3 / R 324 : U2 = V2 / V1 325 : U3 = V3 / V1 326 : C 327 : DO 210 J = KM1, ENORN 328 : T = MATA.RE(K,J,1) + U2 * MATA.RE(K1,J,1) 329 : & + U3 * MATA.RE(K2,J,1) 330 : MATA.RE(K,J,1) = MATA.RE(K,J,1) + T * V1 331 : MATA.RE(K1,J,1) = MATA.RE(K1,J,1) + T * V2 332 : MATA.RE(K2,J,1) = MATA.RE(K2,J,1) + T * V3 333 : T = MATB.RE(K,J,1) + U2 * MATB.RE(K1,J,1) 334 : & + U3 * MATB.RE(K2,J,1) 335 : MATB.RE(K,J,1) = MATB.RE(K,J,1) + T * V1 336 : MATB.RE(K1,J,1) = MATB.RE(K1,J,1) + T * V2 337 : MATB.RE(K2,J,1) = MATB.RE(K2,J,1) + T * V3 338 : 210 CONTINUE 339 : C 340 : IF ( K .EQ. L ) GOTO 220 341 : MATA.RE(K1,KM1,1) = 0.0D0 342 : MATA.RE(K2,KM1,1) = 0.0D0 343 : C 344 : C ********** ZERO B(K+2,K+1) AND B(K+2,K) ********** 345 : 220 continue 346 : S = ABS(MATB.RE(K2,K2,1)) + ABS(MATB.RE(K2,K1,1)) 347 : & + ABS(MATB.RE(K2,K,1)) 348 : IF ( S .EQ. 0.0D0 ) GOTO 240 349 : U1 = MATB.RE(K2,K2,1) / S 350 : U2 = MATB.RE(K2,K1,1) / S 351 : U3 = MATB.RE(K2,K,1) / S 352 : R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1) 353 : V1 = -(U1 + R) / R 354 : V2 = -U2 / R 355 : V3 = -U3 / R 356 : U2 = V2 / V1 357 : U3 = V3 / V1 358 : C 359 : DO 230 I = LOR1, LL 360 : T = MATA.RE(I,K2,1) + U2 * MATA.RE(I,K1,1) 361 : & + U3 * MATA.RE(I,K,1) 362 : MATA.RE(I,K2,1) = MATA.RE(I,K2,1) + T * V1 363 : MATA.RE(I,K1,1) = MATA.RE(I,K1,1) + T * V2 364 : MATA.RE(I,K,1) = MATA.RE(I,K,1) + T * V3 365 : T = MATB.RE(I,K2,1) + U2 * MATB.RE(I,K1,1) 366 : & + U3 * MATB.RE(I,K,1) 367 : MATB.RE(I,K2,1) = MATB.RE(I,K2,1) + T * V1 368 : MATB.RE(I,K1,1) = MATB.RE(I,K1,1) + T * V2 369 : MATB.RE(I,K,1) = MATB.RE(I,K,1) + T * V3 370 : 230 CONTINUE 371 : C 372 : MATB.RE(K2,K,1) = 0.0D0 373 : MATB.RE(K2,K1,1) = 0.0D0 374 : IF ( .NOT. CALV ) GOTO 240 375 : C 376 : DO 235 I = 1, N 377 : T = MATZ.RE(I,K2,1) + U2 * MATZ.RE(I,K1,1) 378 : & + U3 * MATZ.RE(I,K,1) 379 : MATZ.RE(I,K2,1) = MATZ.RE(I,K2,1) + T * V1 380 : MATZ.RE(I,K1,1) = MATZ.RE(I,K1,1) + T * V2 381 : MATZ.RE(I,K,1) = MATZ.RE(I,K,1) + T * V3 382 : 235 CONTINUE 383 : C 384 : C ********** ZERO B(K+1,K) ********** 385 : 240 CONTINUE 386 : S = ABS(MATB.RE(K1,K1,1)) + ABS(MATB.RE(K1,K,1)) 387 : IF ( S .EQ. 0.0D0 ) GOTO 260 388 : U1 = MATB.RE(K1,K1,1) / S 389 : U2 = MATB.RE(K1,K,1) / S 390 : R = SIGN(SQRT(U1*U1+U2*U2),U1) 391 : V1 = -(U1 + R) / R 392 : V2 = -U2 / R 393 : U2 = V2 / V1 394 : C 395 : DO 250 I = LOR1, LL 396 : T = MATA.RE(I,K1,1) + U2 * MATA.RE(I,K,1) 397 : MATA.RE(I,K1,1) = MATA.RE(I,K1,1) + T * V1 398 : MATA.RE(I,K,1) = MATA.RE(I,K,1) + T * V2 399 : T = MATB.RE(I,K1,1) + U2 * MATB.RE(I,K,1) 400 : MATB.RE(I,K1,1) = MATB.RE(I,K1,1) + T * V1 401 : MATB.RE(I,K,1) = MATB.RE(I,K,1) + T * V2 402 : 250 CONTINUE 403 : C 404 : MATB.RE(K1,K,1) = 0.0D0 405 : C 406 : IF ( .NOT. CALV ) GOTO 260 407 : C 408 : DO 255 I = 1, N 409 : T = MATZ.RE(I,K1,1) + U2 * MATZ.RE(I,K,1) 410 : MATZ.RE(I,K1,1) = MATZ.RE(I,K1,1) + T * V1 411 : MATZ.RE(I,K,1) = MATZ.RE(I,K,1) + T * V2 412 : 255 CONTINUE 413 : C 414 : 260 CONTINUE 415 : C 416 : GOTO 70 417 : C 418 : C----------------------------------------------------------------------- 419 : C---- INDICATION D ERREUR SI PAS DE CONVERGENCE EN 30*N ITERATIONS ----- 420 : C----------------------------------------------------------------------- 421 : C 422 : 1000 CERR = EN 423 : C 424 : C----------------------------------------------------------------------- 425 : C------ ON SAUVE EPSB POUR UTILISATION DANS QZVALP ET QZVECP ------ 426 : C----------------------------------------------------------------------- 427 : C 428 : 1001 IF ( N .LE. 1 ) GOTO 1002 429 : MATB.RE(N,1,1) = EPSB 430 : c SEGDES , MATA,MATB,MATZ 431 : C 432 : C======================================================================= 433 : C========================== FIN DES CALCULS ============================ 434 : C======================================================================= 435 : C 436 : if(iimpi.ge.4) then 437 : if(CERR.eq.0) then 438 : write(ioimp,*)'convergence en',ITS,' iterations ITN=',ITN 439 : else 440 : write(ioimp,*)'NON convergence en',ITS,' iterations ITN=',ITN 441 : write(ioimp,*) 'A=' 442 : do iou=1,N 443 : write(ioimp,*) (MATA.RE(iou,jou,1),jou=1,N) 444 : enddo 445 : write(ioimp,*) 'B=' 446 : do iou=1,N 447 : write(ioimp,*) (MATB.RE(iou,jou,1),jou=1,N) 448 : enddo 449 : endif 450 : endif 451 : 452 : 1002 RETURN 453 : END 454 : 455 : 456 : 457 :
© Cast3M 2003 - All rights reserved.
Disclaimer