Sources Esope | Cast3M

qzhess
C QZHESS    SOURCE    PV090527  26/05/04    21:15:05     12529          
      SUBROUTINE QZHESS (MATA, MATB, CALV, MATZ)
c      SUBROUTINE QZHESS (XMATA, XMATB, CALV, XMATZ, NM,N)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
*
*************************************************************************
*                                                                       *
*               MISE SOUS FORME DE HESSENBERG SUPERIEURE DE MATA        *
*           ET   TRIANGULARISATION SUPERIEURE SIMULTANEE DE MATB        *
*     ( PREMIER PAS DE LA METHODE "QZ" POUR LE CALCUL MODAL COMPLEXE )  *
* _____________________________________________________________________ *
*                                                                       *
*                                                                       *
*     DATE : le 23 Mars 1995                                            *
*     AUTEURS : J.ANTUNES , Nicolas BENECH                              *
* _____________________________________________________________________ *
*                                                                       *
*    MODULE(S) APPELANT(S) : VIBRAC                                     *
*                                                                       *
*    MODULE(S) APPELE(S) :                                              *
* _____________________________________________________________________ *
*                                                                       *
*    EN ENTREE :                                                        *
*       - MATA : matrice reelle generale (XMATRI)                       *
*       - MATB : matrice reelle generale (XMATRI)                       *
*       - CALV : indique si l'on calcule ou non les vecteurs propres    *
*                                                                       *
*    EN SORTIE :                                                        *
*       - MATA : matrice sous forme Hessenberg superieure               *
*       - MATB : matrice sous forme triangulaire superieure             *
*       - MATZ : produit des transformations effectuees (si CALV=.TRUE.)*
*                                                                       *
*************************************************************************
*
 
-INC PPARAM
-INC CCOPTIO
-INC SMRIGID
*
      INTEGER I , J , K , L , N , LB , L1 , NM , NK1 , NM1 , NM2
      LOGICAL CALV
*-----------------------------------------------------------------------
      POINTEUR MATA.XMATRI, MATB.XMATRI, MATZ.XMATRI
*      REAL*8 XMATA(NM,N),XMATB(NM,N),XMATZ(NM,N)
*
c       SEGACT , MATA*mod
c       SEGACT , MATB*mod
      NM=MATA.RE(/1)
      N=MATA.RE(/2)
*
*=======================================================================
*========================= DEBUT DES CALCULS ===========================
*=======================================================================
*
*-----------------------------------------------------------------------
*----------------------- INITIALISATION DE  Z  -------------------------
*-----------------------------------------------------------------------
*
      NLIGRD=NM
      NLIGRP=N
      nelrig=1
      rigrel=0
      SEGINI , MATZ
*
      IF ( .NOT. CALV ) GOTO 9
      DO 3 J = 1, N
         DO 2 I = 1, N
            MATZ.RE(I,J,1) = 0.0D0
    2    CONTINUE
         MATZ.RE(J,J,1) = 1.0D0
    3 CONTINUE
*
    9 continue
*
*-----------------------------------------------------------------------
*--------- REDUCTION DE  B  A LA FORME TRIANGULAIRE SUPERIEURE ---------
*-----------------------------------------------------------------------
*
   10 IF ( N .LE. 1 ) GOTO 170
*
*
      NM1 = N - 1
*
      DO 100 L = 1, NM1
         L1 = L + 1
 
c        S = somme des termes de B_il i=l+1..n
         S = 0.0D0
         DO 20 I = L1, N
            S = S + ABS(MATB.RE(I,L,1))
   20    CONTINUE
c        S=0 : la ligne L verifie deja la forme triangulaire infer de B
         IF ( S .EQ. 0.0D0 ) GOTO 100
c        S = somme des termes de B_il i=l..n : ajout du terme diagonal
         S = S + ABS(MATB.RE(L,L,1))
 
c        methode de Householder pour annuler les termes non-nul
c
c        calcul de r = sign(B_ll)*sqrt(sum B_il^2)
c        (rem: r est souvent appelé \alpha dans la litterature)
         R = 0.0D0
         DO 25 I = L, N
            MATB.RE(I,L,1) = MATB.RE(I,L,1) / S
            R = R + MATB.RE(I,L,1)**2
   25    CONTINUE
         R = SIGN(SQRT(R),MATB.RE(L,L,1))
         MATB.RE(L,L,1) = MATB.RE(L,L,1) + R
         RHO = R * MATB.RE(L,L,1)
c        triangularisation effective de B
         DO 50 J = L1, N
            T = 0.0D0
            DO 30 I = L, N
               T = T + MATB.RE(I,L,1) * MATB.RE(I,J,1)
   30       CONTINUE
            T = -T / RHO
            DO 40 I = L, N
               MATB.RE(I,J,1) = MATB.RE(I,J,1) + T * MATB.RE(I,L,1)
   40       CONTINUE
   50    CONTINUE
c        la meme transformation est appliquee a A
         DO 80 J = 1, N
            T = 0.0D0
            DO 60 I = L, N
               T = T + MATB.RE(I,L,1) * MATA.RE(I,J,1)
   60       CONTINUE
            T = -T / RHO
            DO 70 I = L, N
               MATA.RE(I,J,1) = MATA.RE(I,J,1) + T * MATB.RE(I,L,1)
   70       CONTINUE
*
   80    CONTINUE
c        diagonale de B
         MATB.RE(L,L,1) = -S * R
c        mise a 0 des termes post-diagonaux
         DO 90 I = L1, N
            MATB.RE(I,L,1) = 0.0D0
   90    CONTINUE
*
  100 CONTINUE
*
*
C-----------------------------------------------------------------------
C-------- REDUCTION DE  A  A LA FORME DE HESSENBERG SUPERIEURE ---------
C-----------------------------------------------------------------------
c         while keeping B triangular...
*
      IF ( N .EQ. 2 ) GOTO 170
C
      NM2 = N - 2
*
      DO 160 K = 1, NM2
         NK1 = NM1 - K
*
         DO 150 LB = 1, NK1
            L = N - LB
            L1 = L + 1
*
            S = ABS(MATA.RE(L,K,1)) + ABS(MATA.RE(L1,K,1))
            IF ( S .EQ. 0.0D0 ) GOTO 150
            U1 = MATA.RE(L,K,1) / S
            U2 = MATA.RE(L1,K,1) / S
            R = SIGN(SQRT(U1*U1+U2*U2),U1)
            V1 =  -(U1 + R) / R
            V2 = -U2 / R
            U2 = V2 / V1
*
            DO 110 J = K, N
               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
  110       CONTINUE
            MATA.RE(L1,K,1) = 0.0D0
*
            DO 120 J = L, N
               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
  120       CONTINUE
*
            S = ABS(MATB.RE(L1,L1,1)) + ABS(MATB.RE(L1,L,1))
            IF ( S .EQ. 0.0D0 ) GOTO 150
            U1 = MATB.RE(L1,L1,1) / S
            U2 = MATB.RE(L1,L,1) / S
            R = SIGN(SQRT(U1*U1+U2*U2),U1)
            V1 =  -(U1 + R) / R
            V2 = -U2 / R
            U2 = V2 / V1
*
            DO 130 I = 1, L1
               T = MATB.RE(I,L1,1) + U2 * MATB.RE(I,L,1)
               MATB.RE(I,L1,1) = MATB.RE(I,L1,1) + T * V1
               MATB.RE(I,L,1) = MATB.RE(I,L,1) + T * V2
  130       CONTINUE
            MATB.RE(L1,L,1) = 0.0D0
*
              DO 140 I = 1, N
                T = MATA.RE(I,L1,1) + U2 * MATA.RE(I,L,1)
                MATA.RE(I,L1,1) = MATA.RE(I,L1,1) + T * V1
                MATA.RE(I,L,1) = MATA.RE(I,L,1) + T * V2
  140         CONTINUE
*
            IF ( .NOT. CALV ) GOTO 150
*
              DO 145 I = 1, N
                 T = MATZ.RE(I,L1,1) + U2 * MATZ.RE(I,L,1)
                 MATZ.RE(I,L1,1) = MATZ.RE(I,L1,1) + T * V1
                 MATZ.RE(I,L,1) = MATZ.RE(I,L,1) + T * V2
  145         CONTINUE
*
  150    CONTINUE
*
  160 CONTINUE
*
  170 CONTINUE
*
*========================== FIN DES CALCULS ============================
c       SEGDES , MATA,MATB,MATZ
*
      RETURN
      END