C QZREDU SOURCE PV090527 26/04/30 21:16:04 12529 SUBROUTINE QZREDU (MATA, MATB, EPS1, CALV, MATZ, CERR) c SUBROUTINE QZREDU (XMATA, XMATB, EPS1, CALV, XMATZ,NM,N, CERR) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) * ************************************************************************* * * * MISE SOUS FORME QUASI-TRIANGULAIRE SUPERIEURE * * D'UNE MATRICE DE HESSENBERG SUPERIEURE (MATA) * * ET TRIANGULARISATION SUPERIEURE SIMULTANEE DE MATB * * ( DEUXIEME PAS DE LA METHODE "QZ" POUR LE CALCUL MODAL COMPLEXE ) * * _____________________________________________________________________ * * * * * * DATE : le 27 Mars 1995 * * AUTEURS : J.ANTUNES , Nicolas BENECH * * _____________________________________________________________________ * * * * MODULE(S) APPELANT(S) : VIBRAC * * * * MODULE(S) APPELE(S) : * * _____________________________________________________________________ * * * * EN ENTREE : * * - MATA : matrice reelle Hessenberg superieure (XMATRI) * * - MATB : matrice reelle triangulaire superieure (XMATRI) * * - EPS1 : tolerance pour les elements negligeables * * si EPS1 = 0. ou EPS1 < 0., un element est neglige * * s'il est inferieur a l'erreur d'arrondi multiplie * * par la norme de la matrice * * si EPS1 > 0., un element est neglige s'il est inferieur * * a EPS1 multiplie par la norme de la matrice * * (plus rapide mais moins fiable) * * - CALV : indique si l'on calcule ou non les vecteurs propres * * * * EN SORTIE : * * - MATA : matrice quasi-triangulaire superieure * * - MATB : matrice sous forme triangulaire superieure * * - MATZ : produit des transformations effectuees (si CALV=.TRUE.)* * - CERR : code des erreurs * * CERR = 0 : retour normal * * CERR = J : plus de 30*N iterations pour le calcul * * de la valeur propre J * * * ************************************************************************* * -INC PPARAM -INC CCOPTIO -INC SMRIGID C INTEGER I , J , K , L , EN , K1 , K2 , LD , LL , L1 , NA , & ISH , ITN , ITS , KM1 , LM1 , ENM2 , CERR , LOR1 , & ENORN C LOGICAL CALV , NOTLAS C----------------------------------------------------------------------- POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI c REAL*8 XMATA(NM,N),XMATB(NM,N),XMATZ(NM,N) C c SEGACT , MATA*mod c SEGACT , MATB*mod c SEGACT , MATZ*mod NM=MATA.RE(/1) N=MATA.RE(/2) C----------------------------------------------------------------------- C if(iimpi.ge.6) then write(ioimp,*) 'A=' do iou=1,N write(ioimp,*) (MATA.RE(iou,jou,1),jou=1,N) enddo write(ioimp,*) 'B=' do iou=1,N write(ioimp,*) (MATB.RE(iou,jou,1),jou=1,N) enddo endif C*********************************************************************** C*********************************************************************** C*********************************************************************** C C======================================================================= C========================= DEBUT DES CALCULS =========================== C======================================================================= C CERR = 0 C C----------------------------------------------------------------------- C---------------------- CALCUL DE EPSA ET EPSB ---------------------- C----------------------------------------------------------------------- C ANORM = 0.0D0 BNORM = 0.0D0 C DO 30 I = 1, N ANI = 0.0D0 IF ( I .NE. 1 ) ANI = ABS(MATA.RE(I,I-1,1)) BNI = 0.0D0 C DO 20 J = I, N ANI = ANI + ABS(MATA.RE(I,J,1)) BNI = BNI + ABS(MATB.RE(I,J,1)) 20 CONTINUE C IF ( ANI .GT. ANORM ) ANORM = ANI IF ( BNI .GT. BNORM ) BNORM = BNI 30 CONTINUE C IF ( ANORM .EQ. 0.0D0 ) ANORM = 1.0D0 IF ( BNORM .EQ. 0.0D0 ) BNORM = 1.0D0 EP = EPS1 IF ( EP .GT. 0.0D0 ) GOTO 50 C EP = EPSLON(1.0D0) IF (EP .GT. 0.0D0) GO TO 50 C ********** COMPUTE ROUNDOFF LEVEL IF EPSLON(1.) IS ZERO ********** EP = 1.0D0 40 EP = EP / 2.0 IF (1.0D0 + EP .GT. 1.0D0) GO TO 40 50 EPSA = EP * ANORM EPSB = EP * BNORM C ********** REDUCE A TO QUASI-TRIANGULAR FORM, WHILE C KEEPING B TRIANGULAR ********** LOR1 = 1 ENORN = N EN = N c ITN = 30*N c bp 2012-10-02 : on augmente le nombre d'iterations max ITN = 60*N C C----------------------------------------------------------------------- C------------------------- DEBUT DU PAS QZ --------------------------- C----------------------------------------------------------------------- C 60 IF ( EN .LE. 2 ) GOTO 1001 IF ( .NOT. CALV ) ENORN = EN ITS = 0 NA = EN - 1 ENM2 = NA - 1 70 ISH = 2 C C----------------------------------------------------------------------- C------------------- VERIFICATION DE LA CONVERGENCE -------------------- C----------------------------------------------------------------------- C DO 80 LL = 1, EN LM1 = EN - LL L = LM1 + 1 IF ( L .EQ. 1 ) GOTO 95 IF ( ABS(MATA.RE(L,LM1,1)) .LE. EPSA ) GOTO 90 80 CONTINUE C 90 MATA.RE(L,LM1,1) = 0.0D0 IF ( L .LT. NA ) GOTO 95 C ********** 1-BY-1 OR 2-BY-2 BLOCK ISOLATED ********** EN = LM1 GOTO 60 C ********** CHECK FOR SMALL TOP OF B ********** 95 LD = L 100 L1 = L + 1 C B11 = MATB.RE(L,L,1) IF ( ABS(B11) .GT. EPSB ) GOTO 120 MATB.RE(L,L,1) = 0.0D0 C S = ABS(MATA.RE(L,L,1)) + ABS(MATA.RE(L1,L,1)) U1 = MATA.RE(L,L,1) / S U2 = MATA.RE(L1,L,1) / S C C R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1 + R) / R V2 = -U2 / R U2 = V2 / V1 C DO 110 J = L, ENORN T = MATA.RE(L,J,1) + U2 * MATA.RE(L1,J,1) MATA.RE(L,J,1) = MATA.RE(L,J,1) + T * V1 MATA.RE(L1,J,1) = MATA.RE(L1,J,1) + T * V2 T = MATB.RE(L,J,1) + U2 * MATB.RE(L1,J,1) MATB.RE(L,J,1) = MATB.RE(L,J,1) + T * V1 MATB.RE(L1,J,1) = MATB.RE(L1,J,1) + T * V2 110 CONTINUE C IF ( L .NE. 1 ) MATA.RE(L,LM1,1) = -MATA.RE(L,LM1,1) LM1 = L L = L1 GOTO 90 120 CONTINUE A11 = MATA.RE(L,L,1) / B11 A21 = MATA.RE(L1,L,1) / B11 IF ( ISH .EQ. 1 ) GOTO 141 C C----------------------------------------------------------------------- C---------------------- STRATEGIE DE L ITERATION ----------------------- C----------------------------------------------------------------------- C C ********** ITERATION STRATEGY ********** IF ( ITN .EQ. 0 ) GOTO 1000 c bp 2012-10-02 : on change le moment d'appel au "ad hoc shift" dont on c ne trouve pas de justification mathematique mais on la cale c numeriquement sur un cas test problematique : brasero_tube3D.dgibi c l'ideal serait d'avoir un critere qui permet de declencher de shift c (ou bien de passer a linpack et DHGEQZ) c IF ( ITS .EQ. 10 ) GOTO 155 IF ( ITS .EQ. 10*N) GOTO 155 C C ********** DETERMINE TYPE OF SHIFT ********** B22 = MATB.RE(L1,L1,1) IF ( ABS(B22) .LT. EPSB ) B22 = EPSB B33 = MATB.RE(NA,NA,1) IF ( ABS(B33) .LT. EPSB ) B33 = EPSB B44 = MATB.RE(EN,EN,1) IF ( ABS(B44) .LT. EPSB ) B44 = EPSB A33 = MATA.RE(NA,NA,1) / B33 A34 = MATA.RE(NA,EN,1) / B44 A43 = MATA.RE(EN,NA,1) / B33 A44 = MATA.RE(EN,EN,1) / B44 B34 = MATB.RE(NA,EN,1) / B44 T = 0.5D0 * (A43 * B34 - A33 - A44) R = T * T + A34 * A43 - A33 * A44 IF ( R .LT. 0.0D0 ) GOTO 150 C C ********** DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A ********** ISH = 1 R = SQRT(R) SH = -T + R S = -T - R IF ( ABS(S-A44) .LT. ABS(SH-A44) ) SH = S C C ********** LOOK FOR TWO CONSECUTIVE SMALL C SUB-DIAGONAL ELEMENTS OF A. C FOR L=EN-2 STEP -1 UNTIL LD DO -- ********** DO 130 LL = LD, ENM2 L = ENM2 + LD - LL IF ( L .EQ. LD ) GOTO 140 LM1 = L - 1 L1 = L + 1 T = MATA.RE(L,L,1) IF ( ABS(MATB.RE(L,L,1)) .GT. EPSB ) & T = T - SH * MATB.RE(L,L,1) IF (ABS(MATA.RE(L,LM1,1)).LE.ABS(T/MATA.RE(L1,L,1))*EPSA) & GOTO 100 130 CONTINUE 140 continue C 141 A1 = A11 - SH A2 = A21 IF ( L .NE. LD ) MATA.RE(L,LM1,1) = -MATA.RE(L,LM1,1) GOTO 160 C C ********** DETERMINE DOUBLE SHIFT ZEROTH COLUMN OF A ********** 150 continue A12 = MATA.RE(L,L1,1) / B22 A22 = MATA.RE(L1,L1,1) / B22 B12 = MATB.RE(L,L1,1) / B22 A1 = ((A33 - A11) * (A44 - A11) - A34 * A43 + A43 * B34 * A11) & / A21 + A12 - A11 * B12 A2 = (A22 - A11) - A21 * B12 - (A33 - A11) - (A44 - A11) & + A43 * B34 A3 = MATA.RE(L1+1,L1,1) / B22 GOTO 160 C C ********** AD HOC SHIFT ********** 155 A1 = 0.0D0 A2 = 1.0D0 A3 = 1.1605D0 160 ITS = ITS + 1 ITN = ITN - 1 if(iimpi.ge.6) write(*,*)'iteration ',ITS,' A1,A2,A3=',A1,A2,A3 IF ( .NOT. CALV ) LOR1 = LD C C----------------------------------------------------------------------- C------------------------- BOUCLE PRINCIPALE --------------------------- C----------------------------------------------------------------------- C DO 260 K = L, NA NOTLAS = K .NE. NA .AND. ISH .EQ. 2 K1 = K + 1 K2 = K + 2 KM1 = MAX0(K-1,L) LL = MIN0(EN,K1+ISH) IF ( NOTLAS ) GOTO 190 C C ********** ZERO A(K+1,K-1) ********** IF ( K .EQ. L ) GOTO 170 A1 = MATA.RE(K,KM1,1) A2 = MATA.RE(K1,KM1,1) 170 S = ABS(A1) + ABS(A2) IF ( S .EQ. 0.0D0 ) GOTO 70 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 180 J = KM1, ENORN T = MATA.RE(K,J,1) + U2 * MATA.RE(K1,J,1) MATA.RE(K,J,1) = MATA.RE(K,J,1) + T * V1 MATA.RE(K1,J,1) = MATA.RE(K1,J,1) + T * V2 T = MATB.RE(K,J,1) + U2 * MATB.RE(K1,J,1) MATB.RE(K,J,1) = MATB.RE(K,J,1) + T * V1 MATB.RE(K1,J,1) = MATB.RE(K1,J,1) + T * V2 180 CONTINUE C IF ( K .NE. L ) MATA.RE(K1,KM1,1) = 0.0D0 GOTO 240 C C ********** ZERO A(K+1,K-1) AND A(K+2,K-1) ********** 190 IF ( K .EQ. L ) GOTO 200 A1 = MATA.RE(K,KM1,1) A2 = MATA.RE(K1,KM1,1) A3 = MATA.RE(K2,KM1,1) 200 S = ABS(A1) + ABS(A2) + ABS(A3) IF ( S .EQ. 0.0D0 ) GOTO 260 U1 = A1 / S U2 = A2 / S U3 = A3 / S R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1) V1 = -(U1 + R) / R V2 = -U2 / R V3 = -U3 / R U2 = V2 / V1 U3 = V3 / V1 C DO 210 J = KM1, ENORN T = MATA.RE(K,J,1) + U2 * MATA.RE(K1,J,1) & + U3 * MATA.RE(K2,J,1) MATA.RE(K,J,1) = MATA.RE(K,J,1) + T * V1 MATA.RE(K1,J,1) = MATA.RE(K1,J,1) + T * V2 MATA.RE(K2,J,1) = MATA.RE(K2,J,1) + T * V3 T = MATB.RE(K,J,1) + U2 * MATB.RE(K1,J,1) & + U3 * MATB.RE(K2,J,1) MATB.RE(K,J,1) = MATB.RE(K,J,1) + T * V1 MATB.RE(K1,J,1) = MATB.RE(K1,J,1) + T * V2 MATB.RE(K2,J,1) = MATB.RE(K2,J,1) + T * V3 210 CONTINUE C IF ( K .EQ. L ) GOTO 220 MATA.RE(K1,KM1,1) = 0.0D0 MATA.RE(K2,KM1,1) = 0.0D0 C C ********** ZERO B(K+2,K+1) AND B(K+2,K) ********** 220 continue S = ABS(MATB.RE(K2,K2,1)) + ABS(MATB.RE(K2,K1,1)) & + ABS(MATB.RE(K2,K,1)) IF ( S .EQ. 0.0D0 ) GOTO 240 U1 = MATB.RE(K2,K2,1) / S U2 = MATB.RE(K2,K1,1) / S U3 = MATB.RE(K2,K,1) / S R = SIGN(SQRT(U1*U1+U2*U2+U3*U3),U1) V1 = -(U1 + R) / R V2 = -U2 / R V3 = -U3 / R U2 = V2 / V1 U3 = V3 / V1 C DO 230 I = LOR1, LL T = MATA.RE(I,K2,1) + U2 * MATA.RE(I,K1,1) & + U3 * MATA.RE(I,K,1) MATA.RE(I,K2,1) = MATA.RE(I,K2,1) + T * V1 MATA.RE(I,K1,1) = MATA.RE(I,K1,1) + T * V2 MATA.RE(I,K,1) = MATA.RE(I,K,1) + T * V3 T = MATB.RE(I,K2,1) + U2 * MATB.RE(I,K1,1) & + U3 * MATB.RE(I,K,1) MATB.RE(I,K2,1) = MATB.RE(I,K2,1) + T * V1 MATB.RE(I,K1,1) = MATB.RE(I,K1,1) + T * V2 MATB.RE(I,K,1) = MATB.RE(I,K,1) + T * V3 230 CONTINUE C MATB.RE(K2,K,1) = 0.0D0 MATB.RE(K2,K1,1) = 0.0D0 IF ( .NOT. CALV ) GOTO 240 C DO 235 I = 1, N T = MATZ.RE(I,K2,1) + U2 * MATZ.RE(I,K1,1) & + U3 * MATZ.RE(I,K,1) MATZ.RE(I,K2,1) = MATZ.RE(I,K2,1) + T * V1 MATZ.RE(I,K1,1) = MATZ.RE(I,K1,1) + T * V2 MATZ.RE(I,K,1) = MATZ.RE(I,K,1) + T * V3 235 CONTINUE C C ********** ZERO B(K+1,K) ********** 240 CONTINUE S = ABS(MATB.RE(K1,K1,1)) + ABS(MATB.RE(K1,K,1)) IF ( S .EQ. 0.0D0 ) GOTO 260 U1 = MATB.RE(K1,K1,1) / S U2 = MATB.RE(K1,K,1) / S R = SIGN(SQRT(U1*U1+U2*U2),U1) V1 = -(U1 + R) / R V2 = -U2 / R U2 = V2 / V1 C DO 250 I = LOR1, LL T = MATA.RE(I,K1,1) + U2 * MATA.RE(I,K,1) MATA.RE(I,K1,1) = MATA.RE(I,K1,1) + T * V1 MATA.RE(I,K,1) = MATA.RE(I,K,1) + T * V2 T = MATB.RE(I,K1,1) + U2 * MATB.RE(I,K,1) MATB.RE(I,K1,1) = MATB.RE(I,K1,1) + T * V1 MATB.RE(I,K,1) = MATB.RE(I,K,1) + T * V2 250 CONTINUE C MATB.RE(K1,K,1) = 0.0D0 C IF ( .NOT. CALV ) GOTO 260 C DO 255 I = 1, N T = MATZ.RE(I,K1,1) + U2 * MATZ.RE(I,K,1) MATZ.RE(I,K1,1) = MATZ.RE(I,K1,1) + T * V1 MATZ.RE(I,K,1) = MATZ.RE(I,K,1) + T * V2 255 CONTINUE C 260 CONTINUE C GOTO 70 C C----------------------------------------------------------------------- C---- INDICATION D ERREUR SI PAS DE CONVERGENCE EN 30*N ITERATIONS ----- C----------------------------------------------------------------------- C 1000 CERR = EN C C----------------------------------------------------------------------- C------ ON SAUVE EPSB POUR UTILISATION DANS QZVALP ET QZVECP ------ C----------------------------------------------------------------------- C 1001 IF ( N .LE. 1 ) GOTO 1002 MATB.RE(N,1,1) = EPSB c SEGDES , MATA,MATB,MATZ C C======================================================================= C========================== FIN DES CALCULS ============================ C======================================================================= C if(iimpi.ge.4) then if(CERR.eq.0) then write(ioimp,*)'convergence en',ITS,' iterations ITN=',ITN else write(ioimp,*)'NON convergence en',ITS,' iterations ITN=',ITN write(ioimp,*) 'A=' do iou=1,N write(ioimp,*) (MATA.RE(iou,jou,1),jou=1,N) enddo write(ioimp,*) 'B=' do iou=1,N write(ioimp,*) (MATB.RE(iou,jou,1),jou=1,N) enddo endif endif 1002 RETURN END