Télécharger la plaquette Cast3M


fr en

Télécharger qzredu.eso

Retour à la liste

Numérotation des lignes : 

qzredu

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales