C QZVALP SOURCE CB215821 17/11/30 21:17:05 9639 SUBROUTINE QZVALP (MATA, MATB, ALFR, ALFI, BETA, CALV, MATZ) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * ************************************************************************* * * * FIN DE MISE SOUS FORME QUASI-TRIANGULAIRE SUPERIEURE DE MATA * * ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE MATB * * EXTRACTION DES VALEURS PROPRES * * ( TROISIEME PAS DE LA METHODE "QZ" POUR LE CALCUL MODAL COMPLEXE ) * * _____________________________________________________________________ * * * * * * DATE : le 7 Avril 1995 * * AUTEURS : J.ANTUNES , Nicolas BENECH * * _____________________________________________________________________ * * * * MODULE(S) APPELANT(S) : VIBRAC * * * * MODULE(S) APPELE(S) : * * _____________________________________________________________________ * * * * EN ENTREE : * * - MATA : matrice reelle quasi triangulaire superieure (XMATRI) * * - MATB : matrice reelle triangulaire superieure (XMATRI) * * - CALV : logique indiquant si on calcule les vecteurs propres * * - MATZ : produit des transformations effectuees auparavant * * ( si CALV = .TRUE. ) * * * * EN SORTIE : * * - MATA : matrice quasi-triangulaire superieure * * - MATB : matrice sous forme triangulaire superieure * * - ALFR, ALFI, BETA : liste des valeurs propres * * (parties reelle et imaginaire du numerateur + denominateur)* * - MATZ : produit des transformations effectuees (si CALV=.TRUE.)* * * ************************************************************************* * -INC PPARAM -INC CCOPTIO -INC SMRIGID -INC SMLREEL C REAL*8 SWAP INTEGER I , J , EN , NA , NN , ISW C LOGICAL CALV C----------------------------------------------------------------------- C POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI cbp POINTEUR ALFR.XMATRI, ALFI.XMATRI, BETA.XMATRI POINTEUR ALFR.MLREEL, ALFI.MLREEL, BETA.MLREEL C c SEGACT , MATA*mod c SEGACT , MATB*mod c SEGACT , MATZ*mod NM=MATA.RE(/1) N=MATA.RE(/2) C JG = N SEGINI,ALFR,ALFI,BETA C----------------------------------------------------------------------- C========================= DEBUT DES CALCULS =========================== C----------------------------------------------------------------------- C EPSB = MATB.RE(N,1,1) ISW = 1 C C----------------------------------------------------------------------- C-------------------------- BOUCLE PRINCIPALE -------------------------- C----------------------------------------------------------------------- C DO 510 NN = 1, N EN = N + 1 - NN NA = EN - 1 IF ( ISW . EQ . 2 ) GOTO 505 IF ( EN . EQ . 1 ) GOTO 410 SWAP = MATA.RE(EN,NA,1) IF ( SWAP . NE . 0.0D0 ) GOTO 420 C C----------------------------------------------------------------------- C---------------------- UNE VALEUR PROPRE REELLE ----------------------- C----------------------------------------------------------------------- C 410 continue ALFR.PROG(EN) = MATA.RE(EN,EN,1) IF ( MATB.RE(EN,EN,1) . LT . 0.0D0 ) & ALFR.PROG(EN) = -ALFR.PROG(EN) BETA.PROG(EN) = ABS(MATB.RE(EN,EN,1)) ALFI.PROG(EN) = 0.0D0 GOTO 510 C 420 continue SWAP = ABS(MATB.RE(NA,NA,1)) IF ( SWAP . LE . EPSB ) GOTO 455 SWAP = ABS(MATB.RE(EN,EN,1)) IF ( SWAP . GT . EPSB ) GOTO 430 A1 = MATA.RE(EN,EN,1) A2 = MATA.RE(EN,NA,1) BN = 0.0D0 GOTO 435 430 continue AN = ABS(MATA.RE(NA,NA,1)) + ABS(MATA.RE(NA,EN,1)) & + ABS(MATA.RE(EN,NA,1)) + ABS(MATA.RE(EN,EN,1)) BN = ABS(MATB.RE(NA,NA,1)) + ABS(MATB.RE(NA,EN,1)) & + ABS(MATB.RE(EN,EN,1)) A11 = MATA.RE(NA,NA,1) / AN A12 = MATA.RE(NA,EN,1) / AN A21 = MATA.RE(EN,NA,1) / AN A22 = MATA.RE(EN,EN,1) / AN B11 = MATB.RE(NA,NA,1) / BN B12 = MATB.RE(NA,EN,1) / BN B22 = MATB.RE(EN,EN,1) / BN E = A11 / B11 EI = A22 / B22 S = A21 / (B11 * B22) T = (A22 - E * B22) / B22 IF ( ABS(E) . LE . ABS(EI) ) GOTO 431 E = EI T = (A11 - E * B11) / B11 431 C = 0.5D0 * (T - S * B12) D = C * C + S * (A12 - E * B12) IF ( D . LT . 0.0D0 ) GOTO 480 C C----------------------------------------------------------------------- C-------------------- DEUX VALEURS PROPRES REELLES --------------------- C----------------------------------------------------------------------- C E = E + (C +SIGN(SQRT(D),C)) A11 = A11 - E * B11 A12 = A12 - E * B12 A22 = A22 - E * B22 IF ( ABS(A11) + ABS(A12) . LT . & ABS(A21) + ABS(A22) ) GOTO 432 A1 = A12 A2 = A11 GOTO 435 432 A1 = A22 A2 = A21 C C----------------------------------------------------------------------- C------------------------- CHOISIR Z REEL ---------------------------- C----------------------------------------------------------------------- C 435 S = ABS(A1) + ABS(A2) U1 = A1 / S U2 = A2 / S R =SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1 + R) / R V2 = -U2 / R U2 = V2 / V1 C DO 440 I = 1, EN T = MATA.RE(I,EN,1) + U2 * MATA.RE(I,NA,1) MATA.RE(I,EN,1) = MATA.RE(I,EN,1) + T * V1 MATA.RE(I,NA,1) = MATA.RE(I,NA,1) + T * V2 T = MATB.RE(I,EN,1) + U2 * MATB.RE(I,NA,1) MATB.RE(I,EN,1) = MATB.RE(I,EN,1) + T * V1 MATB.RE(I,NA,1) = MATB.RE(I,NA,1) + T * V2 440 CONTINUE C IF ( .NOT. CALV ) GOTO 450 C DO 445 I = 1, N T = MATZ.RE(I,EN,1) + U2 * MATZ.RE(I,NA,1) MATZ.RE(I,EN,1) = MATZ.RE(I,EN,1) + T * V1 MATZ.RE(I,NA,1) = MATZ.RE(I,NA,1) + T * V2 445 CONTINUE C 450 IF ( BN . EQ . 0.0D0 ) GOTO 475 IF ( AN . LT . ABS(E) * BN ) GOTO 455 A1 = MATB.RE(NA,NA,1) A2 = MATB.RE(EN,NA,1) GOTO 460 455 continue A1 = MATA.RE(NA,NA,1) A2 = MATA.RE(EN,NA,1) C C----------------------------------------------------------------------- C------------------------- CHOISIR Q REEL ---------------------------- C----------------------------------------------------------------------- C 460 S = ABS(A1) + ABS(A2) IF ( S . EQ . 0.0D0 ) GOTO 475 U1 = A1 / S U2 = A2 / S R =SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1 + R) / R V2 = -U2 / R U2 = V2 / V1 C DO 470 J = NA, N T = MATA.RE(NA,J,1) + U2 * MATA.RE(EN,J,1) MATA.RE(NA,J,1) = MATA.RE(NA,J,1) + T * V1 MATA.RE(EN,J,1) = MATA.RE(EN,J,1) + T * V2 T = MATB.RE(NA,J,1) + U2 * MATB.RE(EN,J,1) MATB.RE(NA,J,1) = MATB.RE(NA,J,1) + T * V1 MATB.RE(EN,J,1) = MATB.RE(EN,J,1) + T * V2 470 CONTINUE C 475 continue MATA.RE(EN,NA,1) = 0.0D0 ALFR.PROG(NA) = MATA.RE(NA,NA,1) ALFR.PROG(EN) = MATA.RE(EN,EN,1) MATB.RE(EN,NA,1) = 0.0D0 IF ( MATB.RE(NA,NA,1) . LT . 0.0D0 ) & ALFR.PROG(NA) = -ALFR.PROG(NA) IF ( MATB.RE(EN,EN,1) . LT . 0.0D0 ) & ALFR.PROG(EN) = -ALFR.PROG(EN) BETA.PROG(NA) = ABS(MATB.RE(NA,NA,1)) BETA.PROG(EN) = ABS(MATB.RE(EN,EN,1)) ALFI.PROG(EN) = 0.0D0 ALFI.PROG(NA) = 0.0D0 GOTO 505 C C----------------------------------------------------------------------- C------------------- DEUX VALEURS PROPRES COMPLEXES -------------------- C----------------------------------------------------------------------- C 480 E = E + C EI = SQRT(-D) A11R = A11 - E * B11 A11I = EI * B11 A12R = A12 - E * B12 A12I = EI * B12 A22R = A22 - E * B22 A22I = EI * B22 IF ( ABS(A11R) + ABS(A11I) + ABS(A12R) + ABS(A12I) . LT . & ABS(A21) + ABS(A22R) + ABS(A22I) ) GOTO 482 A1 = A12R A1I = A12I A2 = -A11R A2I = -A11I GOTO 485 482 A1 = A22R A1I = A22I A2 = -A21 A2I = 0.0D0 C C----------------------------------------------------------------------- C----------------------- CHOISIR Z COMPLEXE -------------------------- C----------------------------------------------------------------------- C 485 CZ = SQRT(A1*A1+A1I*A1I) IF ( CZ . EQ . 0.0D0 ) GOTO 487 SZR = (A1 * A2 + A1I * A2I) / CZ SZI = (A1 * A2I - A1I * A2) / CZ R = SQRT(CZ*CZ+SZR*SZR+SZI*SZI) CZ = CZ / R SZR = SZR / R SZI = SZI / R GOTO 490 487 SZR = 1.0D0 SZI = 0.0D0 490 IF ( AN . LT . (ABS(E) + EI) * BN ) GOTO 492 A1 = CZ * B11 + SZR * B12 A1I = SZI * B12 A2 = SZR * B22 A2I = SZI * B22 GOTO 495 492 A1 = CZ * A11 + SZR * A12 A1I = SZI * A12 A2 = CZ * A21 + SZR * A22 A2I = SZI * A22 C C----------------------------------------------------------------------- C----------------------- CHOISIR Q COMPLEXE -------------------------- C----------------------------------------------------------------------- C 495 CQ = SQRT(A1*A1+A1I*A1I) IF ( CQ . EQ . 0.0D0 ) GOTO 497 SQR = (A1 * A2 + A1I * A2I) / CQ SQI = (A1 * A2I - A1I * A2) / CQ R = SQRT(CQ*CQ+SQR*SQR+SQI*SQI) CQ = CQ / R SQR = SQR / R SQI = SQI / R GOTO 500 497 SQR = 1.0D0 SQI = 0.0D0 C C----------------------------------------------------------------------- C------------------- CALCUL DES ELEMENTS DIAGONAUX -------------------- C----------------------------------------------------------------------- C 500 SSR = SQR * SZR + SQI * SZI SSI = SQR * SZI - SQI * SZR I = 1 TR = CQ * CZ * A11 + CQ * SZR * A12 + SQR * CZ * A21 & + SSR * A22 TI = CQ * SZI * A12 - SQI * CZ * A21 + SSI * A22 DR = CQ * CZ * B11 + CQ * SZR * B12 + SSR * B22 DI = CQ * SZI * B12 + SSI * B22 GOTO 503 502 I = 2 TR = SSR * A11 - SQR * CZ * A12 - CQ * SZR * A21 & + CQ * CZ * A22 TI = -SSI * A11 - SQI * CZ * A12 + CQ * SZI * A21 DR = SSR * B11 - SQR * CZ * B12 + CQ * CZ * B22 DI = -SSI * B11 - SQI * CZ * B12 503 T = TI * DR - TR * DI J = NA IF ( T . LT . 0.0D0 ) J = EN R = SQRT(DR*DR+DI*DI) BETA.PROG(J) = BN * R ALFR.PROG(J) = AN * (TR * DR + TI * DI) / R ALFI.PROG(J) = AN * T / R IF ( I . EQ . 1 ) GOTO 502 505 ISW = 3 - ISW 510 CONTINUE MATB.RE(N,1,1) = EPSB C c SEGDES , MATA,MATB,MATZ c SEGDES , ALFR,ALFI,BETA c C======================================================================= C========================== FIN DES CALCULS ============================ C======================================================================= RETURN END