Télécharger la plaquette Cast3M


fr en

Télécharger hbmbif.eso

Retour à la liste

Numérotation des lignes : 

hbmbif

  1. C HBMBIF    SOURCE    OF166741  26/05/11    21:15:07     12538          
  2.  
  3. *=======================================================================
  4. * Calcul des points de bifurcation par une methode de Newton-Raphson
  5. *=======================================================================
  6. * TYPC indique le type de calcul a realiser
  7. *      = 'L' pour un point limite
  8. *      = 'B' pour un point de branchement
  9. *      = 'P' pour un doublement de periode
  10. *      = 'N' pour une bifurcation de Neimark-Sacker
  11. *=======================================================================
  12.  
  13.       SUBROUTINE HBMBIF(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP,MTLIAA,
  14.      &                  MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,TYPC,
  15.      &                  KPARNUM,KSORT)
  16.  
  17.       IMPLICIT INTEGER(I-N)
  18.       IMPLICIT REAL*8(A-H,O-Z)
  19.  
  20. -INC PPARAM
  21. -INC CCOPTIO
  22.  
  23. -INC SMLREEL
  24.       POINTEUR fvec.MLREEL, fvec0.MLREEL
  25. -INC TMDYNC
  26.  
  27. *     segment local (a nettoyer)
  28.       SEGMENT MWORK
  29.         INTEGER IPIV33(2*NT+2),INDEIG(2*NT),I4(3*NT+2),I5(2*NT+1)
  30.         REAL*8 AVVP(NT),AVVP2(NT+1),V1(NT)
  31.         REAL*8 EIGRE(NT),EIGIM(NT),QIN(NT),OMEGIN,KAPPA,KAPPA2
  32.         REAL*8 EIGRE2(NT+1),EIGIM2(NT+1)
  33.         REAL*8 VR(NT,NT),VL(NT,NT),WORK(8*NT)
  34.         REAL*8 VR2(NT+1,NT+1),VL2(NT+1,NT+1),WORK2(8*(NT+1))
  35.         REAL*8 JAC1(NT,NT),JAC2(NT,NT)
  36.         REAL*8 RxwPhi(NT),RXXPHI(NT,NT),RwwPhi
  37.         REAL*8 vA1(NT), vA2(NT),vA3(NT),vB1(NT),vB2(NT),vB3(NT)
  38.         REAL*8 RX2(NT,NT),SS,EI(NT)
  39.         REAL*8 JBP(2*NT+2,2*NT+2),solBP(2*NT+2)
  40.         REAL*8 dWB, dXB(NT),dPHIR(NT),dPHII(NT),WORK2w(16*NT)
  41.         REAL*8 res2(NT),Rco(NT,NT),PHIR(NT),PHII(NT),DPsiPsi(NT,NT)
  42.         REAL*8 MREG(NT+1,NT+1),DEL1(NT,NT),DEL2(NT),Mi,Ci,JH(2*NT,2*NT)
  43.         REAL*8 EXRE2w(2*NT),EXIM2w(2*NT),VL2w(2*NT,2*NT),VR2w(2*NT,4*NT)
  44.         REAL*8 dRxDdk1(NT),dRxDdk2(NT)
  45.         REAL*8 Rxmk2D(NT,NT),D10(NT,NT),D1w(NT,NT),D2(NT,NT),VQ(NT)
  46.         REAL*8 RxxPHIR(NT,NT),RxxPHII(NT,NT),RxwPhiR(NT),RxwPhiI(NT)
  47.         REAL*8 JNS(3*NT+2,3*NT+2),solNS(3*NT+2),resV1(NT),resV2(NT)
  48.         REAL*8 JPD(3*NT+1,3*NT+1),solPD(3*NT+1),D1wP(NT,NT)
  49.       ENDSEGMENT
  50.  
  51.       CHARACTER*1 TYPC
  52.       LOGICAL CHECK,CANAL
  53.  
  54.       INTEGER MAXITS,MAXVCP,INFO
  55.  
  56.       REAL*8 ZERO,ONE,TWO,XEPS
  57.       PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0, XEPS=1.D-6)
  58.  
  59. C Tableaux locaux :
  60.       REAL*8 AA(3,3),sol(3),arrRes(3),arrRes2(4),arrRes3(5)
  61.  
  62. C Fonctions BLAS/LAPACK
  63.       REAL*8    DDOT, DNRM2
  64.       EXTERNAL  DDOT, DNRM2
  65.  
  66. *-----------------------------------------------------------------------
  67.       PARAMETER(MAXITS=40)
  68. *       TOLF=1.e-4
  69. *       TOLMIN=1.e-6
  70. *       STPMX=100.
  71. *
  72.       MTQ   = KTQ
  73.       MTKAM = KTKAM
  74.       MTEMP = KTEMP
  75.       PSORT = KSORT
  76.       PARNUM = KPARNUM
  77.       CANAL = JANAL
  78.  
  79.       SEGINI,MWORK
  80.       JG = NT
  81.       SEGINI,fvec,fvec0
  82.  
  83.       nmc = NT
  84.       SEGINI,MdZw
  85. c*      DO J = 1, nmd
  86. c*        DO I = 1, nmd
  87. c*          MdZW.MATRC(I,J)=0.D0
  88. c*        ENDDO
  89. c*      ENDDO
  90.  
  91. *      sauvegarde du point de la courbe de reponse
  92.       OMEGIN = OMEG
  93.       DO I = 1,NT
  94.         QIN(I) = Q1(I)
  95.       ENDDO
  96. *
  97.       CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  98.      &            MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec0,fref0,
  99.      &            .false.)
  100.       CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)
  101.  
  102. c      ----------------------------------------------------------------
  103. c      !!! Localisation NeimarkSacker et PeriodDoubling a debugger !!!
  104. c       if(TYPC.eq.'N'.or.TYPC.eq.'P') then
  105. c          WRITE(IOIMP,*) 'Localisation bifurcation ',TYPC,'non prevue !'
  106. c          CHECK = .true.
  107. c          return
  108. c       endif
  109. c      ----------------------------------------------------------------
  110.  
  111. c  message
  112.       IF (IIMPI.GE.2) WRITE(IOIMP,*) 'Localisation bifurcation ',TYPC
  113.       IF (IIMPI.GE.3) WRITE(IOIMP,990) 0,dnrm2(NT,fvec0.prog,1),fref0
  114.  990  FORMAT('+ hbmbif : iter ',I2,' -> |R|=',E10.3,' ref=',E10.3)
  115.  
  116. *=======================================================================
  117. *                            Initialisation
  118. *=======================================================================
  119.  
  120. * La Jacobienne augmentee a la structure de base suivante:
  121. *   JJ = [    Rx      O      Rw  ]
  122. *        [ (RxPhi)x   Rx   RxwPhi]
  123. *        [   0'     2Phi'    0   ]
  124. * qui est a adapter selon la bifurcation traitee.
  125. *
  126. * Il faut initialiser le vecteur propre Phi.
  127. * Calcul de la Jacobienne Rx = Z(OMEGA) - dFNLdX
  128.       DO J=1,NT
  129.         DO I=1,NT
  130.           RX(I,J) = ZZ(I,J)-JAC(I,J)
  131.         ENDDO
  132.       ENDDO
  133.       CALL COPYMAT(NT,Rx,Rco)
  134.  
  135. *-----------------------------------------------------------------------
  136. *      Bifurcation statique: le vecteur propre annule la Jacobienne
  137. *-----------------------------------------------------------------------
  138.  
  139.       IF (TYPC.EQ.'L') THEN
  140.  
  141.          CALL DGEEV('N','V',NT,Rco,NT,EIGRE,EIGIM,VL,1,VR,NT,
  142.      &                                               WORK,8*NT,INFO)
  143. *        On prend le vecteur associe a la partie reelle la plus petite.
  144.          MINVP=1
  145.          XMINVP=ABS(EIGRE(1))
  146.          DO I = 1,NT
  147.            AVVP(I) = ABS(EIGRE(I))
  148.            IF(AVVP(I).LT.XMINVP) THEN
  149.              MINVP=I
  150.              XMINVP=AVVP(I)
  151.            ENDIF
  152.          ENDDO
  153. c          MINVP = MINLOC(AVVP,NT)
  154. *        Les vecteurs propres sont purement reels dans ce cas.
  155.          DO I = 1,NT
  156.            PHIR(I) = VR(I,MINVP)
  157.            PHII(I) = ZERO
  158.          ENDDO
  159. *        Pas de nouvelle frequence pour les cycles bifurques
  160.          KAPPA = ZERO
  161.          MAXVCP=1
  162.          XMAXVCP=ABS(PHIR(1))
  163.          DO I = 1,NT
  164.            AVVP(I) = ABS(PHIR(I))
  165.            IF(AVVP(I).GT.XMAXVCP) THEN
  166.              MAXVCP=I
  167.              XMAXVCP=AVVP(I)
  168.            ENDIF
  169.          ENDDO
  170. c          MAXVCP = MAXLOC(AVVP,NT)
  171. *        Elements pour la resolution par blocs
  172.          SS = ZERO
  173.          DO I=1,NT
  174.            SS = SS + abs(RX(I,I))
  175.          ENDDO
  176.          SS = SS/NT
  177. *        Penalisation de la Jacobienne (evite mauvais conditionnement)
  178.          CALL COPYMAT(NT,RX,RX2)
  179.          DO I=1,NT
  180.           EI(I) = ZERO
  181.          ENDDO
  182.          RX2(MAXVCP,MAXVCP) = RX2(MAXVCP,MAXVCP) + SS
  183.          EI(MAXVCP) = ONE
  184.       ENDIF
  185.  
  186. *-----------------------------------------------------------------------
  187. *      Branch point ?
  188. *-----------------------------------------------------------------------
  189.       IF (TYPC.EQ.'B') THEN
  190.  
  191.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  192.          IF (BAL.EQ.1) THEN
  193.            DO I=1,NT
  194.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  195.            ENDDO
  196.          END IF
  197. *        Regularisation de la Jacobienne
  198.          CALL HBMDVEC(NT,NHBM,NDDL,Q1,ONE,V1)
  199.          DO I = 1,NT
  200.            MREG(I,NT+1) = Rw(I)
  201.            MREG(NT+1,I) = Rw(I)
  202.          ENDDO
  203.          znormV1 = ONE / dnrm2(NT,V1,1)
  204.          DO I = 1,NT
  205.             V1(I) = V1(I)*znormV1
  206.          ENDDO
  207.          DO J = 1,NT
  208.            DO I = 1,NT
  209.              RX2(I,J) = RX(I,J) + V1(I)*V1(J)
  210.              MREG(I,J) = RX2(I,J)
  211.            ENDDO
  212.          ENDDO
  213.          MREG(NT+1,NT+1) = ZERO
  214.          CALL DGEEV('N','V',NT+1,MREG,NT+1,EIGRE2,EIGIM2,VL2,1,VR2,NT+1,
  215.      &                                              WORK2,8*(NT+1),INFO)
  216. *        On prend le vecteur associe a partie reelle la plus petite.
  217.          MINVP=1
  218.          XMINVP=ABS(EIGRE2(1))
  219.          DO I = 1,NT+1
  220.            AVVP2(I) = ABS(EIGRE2(I))
  221.            IF(AVVP2(I).LT.XMINVP) THEN
  222.              MINVP=I
  223.              XMINVP=AVVP2(I)
  224.            ENDIF
  225.          ENDDO
  226. c          MINVP = MINLOC(AVVP2,NT+1)
  227. *        Les vecteurs propres sont purement reels dans ce cas.
  228.          DO I = 1,NT
  229.             PHIR(I) = VR2(I,MINVP)
  230.             PHII(I) = ZERO
  231.          ENDDO
  232. *        Pas de nouvelle frequence pour les cycles bifurques
  233.          KAPPA = ZERO
  234.       ENDIF
  235.  
  236. *-----------------------------------------------------------------------
  237. *      Bifurcations dynamiques (Period doubling ou Neimark Sacker)
  238. *-----------------------------------------------------------------------
  239. *     Construction de la matrice de Hill, JJ
  240.       IF ((TYPC.EQ.'P').OR.(TYPC.EQ.'N')) THEN
  241.         DO I = 1,NT
  242.           DO J = 1,NT
  243.             D10(I,J)=ZERO
  244.             D1w(I,J)=ZERO
  245.             D2(I,J) =ZERO
  246.           ENDDO
  247.         ENDDO
  248.         DO I=1,2*NHBM+1
  249.           DO J=1,NDDL
  250.             DEL2(NDDL*(I-1)+J) = ONE/XM(J,1)
  251.             D2(NDDL*(I-1)+J,NDDL*(I-1)+J) = XM(J,1)
  252.           ENDDO
  253.         ENDDO
  254.         DO I=1,NDDL
  255.           D10(I,I) = XASM(I,1)
  256.         ENDDO
  257.         DO J = 2,2*NHBM,2
  258.           DO I=1,NDDL
  259.             D10(NDDL*(1+(J-2))+I,NDDL*(1+(J-2))+I) = XASM(I,1)
  260.             D1w(NDDL*(1+(J-2))+I,NDDL*(1+(J-1))+I) = J*XM(I,1)
  261.             D1w(NDDL*(1+(J-1))+I,NDDL*(1+(J-2))+I) = -J*XM(I,1)
  262.             D10(NDDL*(1+(J-1))+I,NDDL*(1+(J-1))+I) = XASM(I,1)
  263.           ENDDO
  264.         ENDDO
  265. *        DEL1 = DEL10 + OMEG*DEL1w
  266.         CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1)
  267. *     JJ = [DEL2\DEL1   DEL2\Rx]
  268. *          [     -I         0  ]
  269. *
  270.          DO  I=1,NT
  271.            DO  J=1,NT
  272. *            JJ_11 = DEL2\DEL1
  273.              JH(I,J) = -DEL1(I,J)*DEL2(J)
  274. *            JJ_12 = DEL2\Rx
  275.              JH(I,NT+J) = -RX(I,J)*DEL2(J)
  276. *            JJ_21 = [I]
  277.              IF (I.EQ.J) THEN
  278.                JH(NT+I,J) = ONE
  279.              ELSE
  280.                JH(NT+I,J) = ZERO
  281.              ENDIF
  282. *            JJ_22 = [0]
  283.              JH(NT+I,NT+J) = ZERO
  284.            ENDDO
  285.          ENDDO
  286.          CALL DGEEV('N','V',2*NT,JH,2*NT,EXRE2w,EXIM2w,VL2w,1,VR2w,2*NT,
  287.      &                                                 WORK2,16*NT,INFO)
  288.          CALL HBMORDO(2*NT,NDDL,OMEG,EXIM2w,INDEIG)
  289. *         WRITE(*,*) '   REAL   ', '   IMAG   ', '   INDX   '
  290. *         DO I = 1,2*NT
  291. *           WRITE(*,*) EXRE2w(INDEIG(I)),EXIM2w(I),INDEIG(I)
  292. *         ENDDO
  293. *        WRITE(*,*) 'INDXS:',(INDEIG(IOU),IOU=1,2*NDDL)
  294. *
  295. *        On choisit le vecteur associe a la valeur propre avec partie
  296. *        imaginaire positive du mode instable
  297.          MINVP=1
  298.          DO I=1,2*NDDL
  299.            IF ((EXRE2w(INDEIG(I)).GT.ZERO).AND.(EXIM2w(I).GT.ZERO)) THEN
  300.               MINVP=INDEIG(I)
  301.               INDK =I
  302. *             WRITE(*,*) 'Mode instable:',MINVP
  303.             ENDIF
  304.          ENDDO
  305.          DO I=1,NT
  306.            PHIR(I) = VR2w(I,MINVP)
  307.            PHII(I) = VR2w(I,MINVP+1)
  308.          ENDDO
  309. *         WRITE(*,*) 'PHI_R  ', '  PHI_I'
  310. *         DO IOU = 1,NT
  311. *           WRITE(*,*) PHIR(IOU),PHII(IOU)
  312. *         ENDDO
  313. *         STOP
  314.          MAXVCP=1
  315.          XMAXVCP=ABS(PHII(1))
  316.          DO I = 1,NT
  317.            AVVP(I) = ABS(PHII(I))
  318.            IF (AVVP(I).GT.XMAXVCP) THEN
  319.              MAXVCP=I
  320.              XMAXVCP=AVVP(I)
  321.            ENDIF
  322.          ENDDO
  323. c          MAXVCP = MAXLOC(AVVP,NT)
  324.          DO I = 1,NT
  325.              VQ(I) = ZERO
  326.          ENDDO
  327.          VQ(MAXVCP) = ONE
  328. *        ---------------------
  329.          IF (TYPC.EQ.'P') THEN
  330. *        Nouvelle frequence fixee a OMEGA/2
  331.            KAPPA = OMEG/TWO
  332.          ENDIF
  333. *        ---------------------
  334.          IF (TYPC.EQ.'N') THEN
  335. *        Nouvelle frequence incommensurable KAPPA
  336.             KAPPA = ABS(EXIM2w(INDK))
  337. *            write(*,*) 'KAPPA^0=',KAPPA
  338.          ENDIF
  339. *         STOP
  340.        ENDIF
  341. *=======================================================================
  342. *                              Iterations
  343. *=======================================================================
  344.         DO its=1,MAXITS
  345. *-----------------------------------------------------------------------
  346. *        CAS L: POINT LIMITE
  347. *-----------------------------------------------------------------------
  348.          IF (TYPC.EQ.'L') THEN
  349. ***      Calcul des derivees
  350. *        RxxPhi
  351.          EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS)
  352. *        Perturbation de la Jacobienne: +EPS
  353.          DO JJ = 1,NT
  354.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  355.          ENDDO
  356.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  357.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  358.      &              .false.)
  359.          CALL COPYMAT(NT,JAC,JAC1)
  360. *        Perturbation de la Jacobienne: -EPS
  361.          DO JJ = 1,NT
  362.            Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ)
  363.          ENDDO
  364.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  365.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  366.      &              .false.)
  367.          CALL COPYMAT(NT,JAC,JAC2)
  368. *        retour a la valeur initiale
  369.          DO JJ = 1,NT
  370.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  371.          ENDDO
  372. *        Differences centrees
  373.          r_z = ONE / (TWO*EPS1)
  374.          DO I = 1,NT
  375.            DO J = 1,NT
  376.              RxxPHI(I,J) = (JAC2(I,J)-JAC1(I,J))*r_z
  377.            ENDDO
  378.          ENDDO
  379. *        RxwPhi
  380.          CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,MdZw)
  381.          DO I=1,NT
  382.            r_z = ZERO
  383.            DO J = 1,NT
  384.              r_z = r_z + MdZw.MATRC(I,J)*PHIR(J)
  385.            ENDDO
  386.            RxwPhi(I) = r_z
  387.          ENDDO
  388. *        Rw
  389.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  390.          IF (BAL.EQ.1) THEN
  391.            r_z = TWO*OMEG
  392.            DO I=1,NT
  393.              Rw(I) = Rw(I) - r_z * FEXA(I)
  394.            ENDDO
  395.          END IF
  396. *        Resolution par blocs
  397. *        A1
  398.          DO I = 1,NT
  399.            vA1(I) = -fvec0.PROG(I)
  400.          ENDDO
  401.          CALL COPYMAT(NT,RX2,Rco)
  402.          CALL DGESV(NT,1,Rco,NT,IPIV,vA1,NT,INFO)
  403. *        A2
  404.          DO I = 1,NT
  405.            vA2(I) = -Rw(I)
  406.          ENDDO
  407.          CALL COPYMAT(NT,RX2,Rco)
  408.          CALL DGESV(NT,1,Rco,NT,IPIV,vA2,NT,INFO)
  409. *        A3
  410.          DO I = 1,NT
  411.            vA3(I) = EI(I)
  412.          ENDDO
  413.          CALL COPYMAT(NT,RX2,Rco)
  414.          CALL DGESV(NT,1,Rco,NT,IPIV,vA3,NT,INFO)
  415. *        B1
  416.          DO I = 1,NT
  417.            r_z = ZERO
  418.            DO J = 1,NT
  419.              r_z = r_z + RxxPHI(I,J)*vA1(J) + RX(I,J)*PHIR(J)
  420.            ENDDO
  421.            vB1(I) = -r_z
  422.          ENDDO
  423.          CALL COPYMAT(NT,RX2,Rco)
  424.          CALL DGESV(NT,1,Rco,NT,IPIV,vB1,NT,INFO)
  425. *        B2
  426.          DO I = 1,NT
  427.            r_z = ZERO
  428.            DO J = 1,NT
  429.              r_z = r_z + RxxPHI(I,J)*vA2(J)
  430.            ENDDO
  431.            vB2(I) = -(r_z + RxwPhi(I))
  432.          ENDDO
  433.          CALL COPYMAT(NT,RX2,Rco)
  434.          CALL DGESV(NT,1,Rco,NT,IPIV,vB2,NT,INFO)
  435. *        B3
  436.          DO I = 1,NT
  437.            r_z = ZERO
  438.            DO J = 1,NT
  439.              r_z = r_z + RxxPHI(I,J)*vA3(J)
  440.            ENDDO
  441.            vB3(I) = -r_z
  442.          ENDDO
  443.          CALL COPYMAT(NT,RX2,Rco)
  444.          CALL DGESV(NT,1,Rco,NT,IPIV,vB3,NT,INFO)
  445. *
  446. *        Construction de la matrice du probleme
  447.          AA(1,1) = vA2(MAXVCP)
  448.          AA(1,2) = SS*vA3(MAXVCP)-ONE
  449.          AA(1,3) = ZERO
  450.          AA(2,1) = vB2(MAXVCP)
  451.          AA(2,2) = SS*vB3(MAXVCP)
  452.          AA(2,3) = SS*vA3(MAXVCP)-ONE
  453.          AA(3,1) = DDOT(NT,PHIR,1,vB2,1)
  454.          AA(3,2) = SS*(DDOT(NT,PHIR,1,vB3,1))
  455.          AA(3,3) = SS*(DDOT(NT,PHIR,1,vA3,1))
  456. *        Construction du cote droit
  457.          sol(1) = -vA1(MAXVCP)
  458.          sol(2) = -vB1(MAXVCP)
  459.          sol(3) = 0.5D0*(ONE-DDOT(NT,PHIR,1,PHIR,1))
  460.      &          -DDOT(NT,PHIR,1,vB1,1)
  461. *        Resolution
  462.          CALL DGESV(3,1,AA,3,IPIV,sol,3,INFO)
  463. *        Corrections
  464.          dWB = sol(1)
  465.          beta1 = sol(2)
  466.          beta2 = sol(3)
  467.          DO I = 1,NT
  468.           dXB(I) = vA1(I) + dWB*vA2(I) +SS*beta1*vA3(I)
  469.           dPHIR(I)= vB1(I)+ dWB*vB2(I) +SS*beta1*vB3(I) +SS*beta2*vA3(I)
  470.           Q1(I) = Q1(I) + dXB(I)
  471.           PHIR(I) = PHIR(I) + dPHIR(I)
  472.          ENDDO
  473.          OMEG = OMEG + dWB
  474.          ENDIF
  475. *-----------------------------------------------------------------------
  476. *        CAS B: POINT DE BRANCHEMENT
  477. *-----------------------------------------------------------------------
  478.          IF (TYPC.EQ.'B') THEN
  479. ***      Calcul des derivees
  480. *        RxxPhi
  481.          EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS)
  482. *        Perturbation de la Jacobienne: +EPS
  483.          DO JJ = 1,NT
  484.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  485.          ENDDO
  486.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  487.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  488.      &            .false.)
  489.          CALL COPYMAT(NT,JAC,JAC1)
  490. *        Perturbation de la Jacobienne: -EPS
  491.          DO JJ = 1,NT
  492.            Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ)
  493.          ENDDO
  494.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  495.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  496.      &            .false.)
  497.          CALL COPYMAT(NT,JAC,JAC2)
  498.          DO JJ = 1,NT
  499.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  500.          ENDDO
  501. c        DPsiPsi
  502.          CALL HBMDPP(NT,NDDL,Q1,PHIR,DPsiPsi)
  503. *        Differences centrees
  504.          DO I = 1,NT
  505.            DO J = 1,NT
  506.            RxxPHI(I,J)=((JAC2(I,J)-JAC1(I,J))/(TWO*EPS1)) + DPsiPsi(I,J)
  507.            ENDDO
  508.          ENDDO
  509. *        RxwPhi
  510.          CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,MdZw)
  511.          DO I=1,NT
  512.            r_z = ZERO
  513.            DO J = 1,NT
  514.              r_z = r_z + MdZw.MATRC(I,J)*PHIR(J)
  515.            ENDDO
  516.            RxwPhi(I) = r_z
  517.          ENDDO
  518. *        Rw
  519.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  520.          IF (BAL.EQ.1) THEN
  521.            DO I=1,NT
  522.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  523.            ENDDO
  524.          END IF
  525. *        Rww*Phi
  526.          CALL HBMRWW(NT,NDDL,Q1,PHIR,XM,RwwPhi)
  527.          IF (BAL.EQ.1) THEN
  528.            DO I=1,NT
  529.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  530.              RwwPhi = RwwPhi - TWO*FEXA(I)*PHIR(I)
  531.            ENDDO
  532.          END IF
  533. *        Construction de la matrice du probleme et du cote droit
  534.          solBP(2*NT+1) = ZERO
  535.          DO I = 1,NT
  536.            solBP(NT+I) = ZERO
  537.            DO J = 1,NT
  538.              JBP(I,J) = RX(I,J)
  539.              JBP(NT+I,J) = RxxPHI(I,J)
  540.              JBP(I,NT+J) = ZERO
  541.              JBP(NT+I,NT+J) = RX2(I,J)
  542.              solBP(NT+I) = solBP(NT+I) + RX2(I,J)*PHIR(J)
  543.            ENDDO
  544.            JBP(2*NT+1,NT+I) = RxwPhi(I)
  545.            JBP(2*NT+1,I) = Rw(I)
  546.            JBP(2*NT+2,I) = ZERO
  547.            JBP(2*NT+2,NT+I) = TWO*PHIR(I)
  548.            JBP(I,2*NT+1) = Rw(I)
  549.            JBP(NT+I,2*NT+1) =  RxwPhi(I)
  550.            JBP(I,2*NT+2) = PHIR(I)
  551.            JBP(NT+I,2*NT+2) = ZERO
  552.            solBP(I) = -fvec0.PROG(I)
  553.            solBP(NT+I) = -solBP(NT+I)
  554.          ENDDO
  555.          solBP(2*NT+1) = -(DDOT(NT,Rw,1,PHIR,1))
  556.          solBP(2*NT+2) = -(DDOT(NT,PHIR,1,PHIR,1)-ONE)
  557.          JBP(2*NT+1,2*NT+1) = RwwPhi
  558.          JBP(2*NT+1,2*NT+2) = ZERO
  559.          JBP(2*NT+2,2*NT+1) = ZERO
  560.          JBP(2*NT+2,2*NT+2) = ZERO
  561. *
  562. *        Resolution
  563.          CALL DGESV(2*NT+2,1,JBP,2*NT+2,IPIV33,solBP,2*NT+2,INFO)
  564. *        Corrections
  565.          dWB = solBP(2*NT+1)
  566.          DO I = 1,NT
  567.           dXB(I) = solBP(I)
  568.           dPHIR(I)= solBP(NT+I)
  569.           Q1(I) = Q1(I) + dXB(I)
  570.           PHIR(I) = PHIR(I) + dPHIR(I)
  571.          ENDDO
  572.          OMEG = OMEG + dWB
  573.          ENDIF
  574. *-----------------------------------------------------------------------
  575. *        CAS P: DOUBLEMENT DE PERIODE
  576. *-----------------------------------------------------------------------
  577.          IF (TYPC.EQ.'P') THEN
  578. ***      Calcul des derivees
  579. *        RxxPhiR
  580.          EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS)
  581. *        Perturbation de la Jacobienne: +EPS1
  582.          DO JJ = 1,NT
  583.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  584.          ENDDO
  585.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  586.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  587.      &            .false.)
  588.          CALL COPYMAT(NT,JAC,JAC1)
  589. *        Perturbation de la Jacobienne: -EPS1
  590.          DO JJ = 1,NT
  591.            Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ)
  592.          ENDDO
  593.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  594.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  595.      &            .false.)
  596.          CALL COPYMAT(NT,JAC,JAC2)
  597.          DO JJ = 1,NT
  598.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  599.          ENDDO
  600. *        Differences centrees
  601.          DO I = 1,NT
  602.            DO J = 1,NT
  603.              RxxPHIR(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS1)
  604.            ENDDO
  605.          ENDDO
  606. *        RxxPhiI
  607.          EPS2 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHII,1) + XEPS)
  608. *        Perturbation de la Jacobienne: +EPS2
  609.          DO JJ = 1,NT
  610.            Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ)
  611.          ENDDO
  612.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  613.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  614.      &            .false.)
  615.          CALL COPYMAT(NT,JAC,JAC1)
  616. *        Perturbation de la Jacobienne: -EPS2
  617.          DO JJ = 1,NT
  618.            Q1(JJ) = Q1(JJ) - TWO*EPS2*PHII(JJ)
  619.          ENDDO
  620.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  621.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  622.      &            .false.)
  623.          CALL COPYMAT(NT,JAC,JAC2)
  624.          DO JJ = 1,NT
  625.            Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ)
  626.          ENDDO
  627. *        Differences centrees
  628.          DO I = 1,NT
  629.            DO J = 1,NT
  630.              RxxPHII(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS2)
  631.            ENDDO
  632.          ENDDO
  633. *        Rw
  634.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  635.          IF (BAL.EQ.1) THEN
  636.            DO I=1,NT
  637.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  638.            ENDDO
  639.          END IF
  640. *        Calcul de: (Rxw-(OMEG/2)*D2)*PhiR-(0.5*D10+OMEG*D1w)*PhiI et
  641. *                   (Rxw-(OMEG/2)*D2)*PhiI+(0.5*D10+OMEG*D1w)*PhiR
  642.          CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,MdZw)
  643.          CALL COLIMA(NT,NT,DEL1,D10,ONE,-0.5D0,D1wP)
  644.          DO I=1,NT
  645.           RxwPhiR(I) = ZERO
  646.           RxwPhiI(I) = ZERO
  647.           DO J = 1,NT
  648.            RxwPhiR(I)=RxwPhiR(I)+(MdZw.MATRC(I,J)-KAPPA*D2(I,J))*PHIR(J)
  649.      &                                               - D1wP(I,J)*PHII(J)
  650.            RxwPhiI(I)=RxwPhiI(I)+(MdZw.MATRC(I,J)-KAPPA*D2(I,J))*PHII(J)
  651.      &                                               + D1wP(I,J)*PHIR(J)
  652.           ENDDO
  653.          ENDDO
  654. *        Rxmk2D = Rx-(KAPPA**2)*D2
  655.          CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D)
  656. *        Construction de la matrice du probleme et du cote droit
  657.          DO I = 1,NT
  658.            solPD(NT+I) = ZERO
  659.            solPD(2*NT+I) = ZERO
  660.            DO J = 1,NT
  661.              JPD(I,J) = RX(I,J)
  662.              JPD(NT+I,J) = RxxPHIR(I,J)
  663.              JPD(2*NT+I,J) = RxxPHII(I,J)
  664.              JPD(I,NT+J) = ZERO
  665.              JPD(I,2*NT+J) = ZERO
  666.              JPD(NT+I,NT+J) = Rxmk2D(I,J)
  667.              JPD(NT+I,2*NT+J) = -KAPPA*DEL1(I,J)
  668.              JPD(2*NT+I,NT+J) = KAPPA*DEL1(I,J)
  669.              JPD(2*NT+I,2*NT+J) = Rxmk2D(I,J)
  670.              solPD(NT+I)=solPD(NT+I)+Rxmk2D(I,J)*PHIR(J)
  671.      &                                          -KAPPA*DEL1(I,J)*PHII(J)
  672.              solPD(2*NT+I)=solPD(2*NT+I)+Rxmk2D(I,J)*PHII(J)
  673.      &                                          +KAPPA*DEL1(I,J)*PHIR(J)
  674.            ENDDO
  675.            JPD(     I,3*NT+1) = Rw(I)
  676.            JPD(  NT+I,3*NT+1) = RxwPhiR(I)
  677.            JPD(2*NT+I,3*NT+1) = RxwPhiI(I)
  678.            JPD(3*NT+1,I)      = ZERO
  679.            JPD(3*NT+1,2*NT+1) = VQ(I)
  680.            JPD(3*NT+1,NT+1)   = ZERO
  681.            solPD(I)           = -fvec0.PROG(I)
  682.            solPD(NT+I)        = -solNS(NT+I)
  683.            solPD(2*NT+I)      = -solNS(2*NT+I)
  684.          ENDDO
  685.          solPD(3*NT+1) = -(DDOT(NT,VQ,1,PHII,1)-0.5D0)
  686.          JNS(3*NT+1,3*NT+1) = ZERO
  687. *        Resolution
  688.          CALL DGESV(3*NT+1,1,JPD,3*NT+1,I5,solPD,3*NT+1,INFO)
  689. *        Corrections
  690.          dWB = solPD(3*NT+1)
  691.          DO I = 1,NT
  692.           dXB(I) = solPD(I)
  693.           dPHIR(I)= solPD(NT+I)
  694.           dPHII(I)= solPD(2*NT+I)
  695.           Q1(I) = Q1(I) + dXB(I)
  696.           PHIR(I) = PHIR(I) + dPHIR(I)
  697.           PHII(I) = PHII(I) + dPHII(I)
  698.          ENDDO
  699.          OMEG = OMEG + dWB
  700.          KAPPA= OMEG/TWO
  701.          END IF
  702. *-----------------------------------------------------------------------
  703. *        CAS N: BIFURCATION DE NEIMARK-SACKER
  704. *-----------------------------------------------------------------------
  705.          IF (TYPC.EQ.'N') THEN
  706. ***      Calcul des derivees
  707. *        RxxPhiR
  708.          EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS)
  709. *        Perturbation de la Jacobienne: +EPS1
  710.          DO JJ = 1,NT
  711.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  712.          ENDDO
  713.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  714.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  715.      &            .false.)
  716.          CALL COPYMAT(NT,JAC,JAC1)
  717. *        Perturbation de la Jacobienne: -EPS1
  718.          DO JJ = 1,NT
  719.            Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ)
  720.          ENDDO
  721.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  722.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  723.      &            .false.)
  724.          CALL COPYMAT(NT,JAC,JAC2)
  725.          DO JJ = 1,NT
  726.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  727.          ENDDO
  728. *        Differences centrees
  729.          DO I = 1,NT
  730.            DO J = 1,NT
  731.              RxxPHIR(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS1)
  732.            ENDDO
  733.          ENDDO
  734. *        RxxPhiI
  735.          EPS2 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHII,1) + XEPS)
  736. *        Perturbation de la Jacobienne: +EPS2
  737.          DO JJ = 1,NT
  738.            Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ)
  739.          ENDDO
  740.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  741.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  742.      &            .false.)
  743.          CALL COPYMAT(NT,JAC,JAC1)
  744. *        Perturbation de la Jacobienne: -EPS2
  745.          DO JJ = 1,NT
  746.            Q1(JJ) = Q1(JJ) - TWO*EPS2*PHII(JJ)
  747.          ENDDO
  748.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  749.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  750.      &            .false.)
  751.          CALL COPYMAT(NT,JAC,JAC2)
  752.          DO JJ = 1,NT
  753.            Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ)
  754.          ENDDO
  755. *        Differences centrees
  756.          DO I = 1,NT
  757.            DO J = 1,NT
  758.              RxxPHII(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS2)
  759.            ENDDO
  760.          ENDDO
  761. *        RxwPhiR-KAPPA*D1w*PhiI, RxwPhiI+KAPPA*D1w*PhiR
  762.          CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,MdZw)
  763.          DO I=1,NT
  764.           RxwPhiR(I) = ZERO
  765.           RxwPhiI(I) = ZERO
  766.           dRxDdk1(I) = ZERO
  767.           dRxDdk2(I) = ZERO
  768.           DO J = 1,NT
  769.            RxwPhiR(I)=RxwPhiR(I)+MdZw.MATRC(I,J)*PHIR(J)
  770.      &                          -KAPPA*D1w(I,J)*PHII(J)
  771.            RxwPhiI(I)=RxwPhiI(I)+MdZw.MATRC(I,J)*PHII(J)
  772.      &                          +KAPPA*D1w(I,J)*PHIR(J)
  773.            dRxDdk1(I)=dRxDdk1(I)-TWO*KAPPA*D2(I,J)*PHIR(J)
  774.      &                                                -DEL1(I,J)*PHII(J)
  775.            dRxDdk2(I)=dRxDdk2(I)-TWO*KAPPA*D2(I,J)*PHII(J)
  776.      &                                                +DEL1(I,J)*PHIR(J)
  777.           ENDDO
  778.          ENDDO
  779. *        Rxmk2D = Rx-(KAPPA**2)*D2
  780.          CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D)
  781. *        Rw
  782.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  783.          IF (BAL.EQ.1) THEN
  784.            DO I=1,NT
  785.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  786.            ENDDO
  787.          END IF
  788. *        Construction de la matrice du probleme et du cote droit
  789.          DO I = 1,NT
  790.            solNS(NT+I) = ZERO
  791.            solNS(2*NT+I) = ZERO
  792.            DO J = 1,NT
  793.              JNS(I,J) = RX(I,J)
  794.              JNS(NT+I,J) = RxxPHIR(I,J)
  795.              JNS(2*NT+I,J) = RxxPHII(I,J)
  796.              JNS(I,NT+J) = ZERO
  797.              JNS(I,2*NT+J) = ZERO
  798.              JNS(NT+I,NT+J) = Rxmk2D(I,J)
  799.              JNS(NT+I,2*NT+J) = -KAPPA*DEL1(I,J)
  800.              JNS(2*NT+I,NT+J) = KAPPA*DEL1(I,J)
  801.              JNS(2*NT+I,2*NT+J) = Rxmk2D(I,J)
  802.              solNS(NT+I)=solNS(NT+I)+Rxmk2D(I,J)*PHIR(J)
  803.      &                                          -KAPPA*DEL1(I,J)*PHII(J)
  804.              solNS(2*NT+I)=solNS(2*NT+I)+Rxmk2D(I,J)*PHII(J)
  805.      &                                          +KAPPA*DEL1(I,J)*PHIR(J)
  806.            ENDDO
  807.            JNS(     I,3*NT+1) = ZERO
  808.            JNS(  NT+I,3*NT+1) = dRxDdk1(I)
  809.            JNS(2*NT+I,3*NT+1) = dRxDdk2(I)
  810.            JNS(     I,3*NT+2) = Rw(I)
  811.            JNS(  NT+I,3*NT+2) = RxwPhiR(I)
  812.            JNS(2*NT+I,3*NT+2) = RxwPhiI(I)
  813.            JNS(3*NT+1,I)      = ZERO
  814.            JNS(3*NT+1,NT+I)   = VQ(I)
  815.            JNS(3*NT+1,2*NT+I) = ZERO
  816.            JNS(3*NT+2,I)      = ZERO
  817.            JNS(3*NT+2,NT+I)   = ZERO
  818.            JNS(3*NT+2,2*NT+I) = VQ(I)
  819.            solNS(I)           = -fvec0.PROG(I)
  820.            solNS(NT+I)        = -solNS(NT+I)
  821.            solNS(2*NT+I)      = -solNS(2*NT+I)
  822.          ENDDO
  823.  
  824.          solNS(3*NT+1) = -(DDOT(NT,VQ,1,PHIR,1)-ONE)
  825.          solNS(3*NT+2) = -(DDOT(NT,VQ,1,PHII,1))
  826.          JNS(3*NT+1,3*NT+1) = ZERO
  827.          JNS(3*NT+1,3*NT+2) = ZERO
  828.          JNS(3*NT+2,3*NT+1) = ZERO
  829.          JNS(3*NT+2,3*NT+2) = ZERO
  830. *        Resolution
  831.          CALL DGESV(3*NT+2,1,JNS,3*NT+2,I4,solNS,3*NT+2,INFO)
  832.          DO IOU=1,3*NT+2
  833.          ENDDO
  834. *        Corrections
  835.          dWB = solNS(3*NT+2)
  836.          dK  = solNS(3*NT+1)
  837.          DO I = 1,NT
  838.           dXB(I) = solNS(I)
  839.           dPHIR(I)= solNS(NT+I)
  840.           dPHII(I)= solNS(2*NT+I)
  841.           Q1(I) = Q1(I) + dXB(I)
  842.           PHIR(I) = PHIR(I) + dPHIR(I)
  843.           PHII(I) = PHII(I) + dPHII(I)
  844.          ENDDO
  845.          OMEG = OMEG + dWB
  846.          KAPPA= KAPPA+ dK
  847. *         write(*,*) 'w +',dWB,'=',OMEG
  848. *         write(*,*) 'kappa +',dK,'=',KAPPA
  849.          END IF
  850. *=======================================================================
  851. *          Evaluation du residu
  852. *=======================================================================
  853.  
  854. *          a) R
  855.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  856.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec0,fref,
  857.      &            .false.)
  858.  
  859. *-----------------------------------------------------------------------
  860.          IF (TYPC.EQ.'L') THEN
  861. *-----------------------------------------------------------------------
  862. *          b) Rxs*Phi
  863.            CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)
  864.            DO I=1,NT
  865.              res2(I) = ZERO
  866.              DO J=1,NT
  867.                RX(I,J) = ZZ(I,J) - JAC(I,J)
  868.                res2(I) = res2(I) + RX(I,J)*PHIR(J)
  869.              ENDDO
  870.            ENDDO
  871. *          c) norm(Phi)= 1.
  872.            res3 = DDOT(NT,PHIR,1,PHIR,1)-ONE
  873. *          Y = ( R  ;  RQ * \phi  ;  \phi*\phi-1 )
  874.            arrRes(1) = dnrm2(NT,fvec0.prog(1),1)
  875.            arrRes(2) = dnrm2(NT,res2,1)
  876.            arrRes(3) = abs(res3)
  877.            res = MAXVAL(arrRes,1)
  878.            IF (IIMPI.GE.3)
  879.      &     WRITE(IOIMP,991) its,(arrRes(iou),iou=1,3),res
  880.  991       FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';',
  881.      &     E10.3,';','} -> res=',E10.3)
  882.          ENDIF
  883.  
  884. *-----------------------------------------------------------------------
  885.          IF (TYPC.EQ.'B') THEN
  886. *-----------------------------------------------------------------------
  887. *          b) Rxs*Phi
  888.            CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)
  889.            CALL HBMDVEC(NT,NHBM,NDDL,Q1,ONE,V1)
  890.            znormV1 = ONE / dnrm2(NT,V1,1)
  891.            DO I = 1,NT
  892.              V1(I) = V1(I)*znormV1
  893.            ENDDO
  894.            DO I = 1,NT
  895.              r_z = ZERO
  896.              DO J = 1,NT
  897.                RX(I,J)  = ZZ(I,J) - JAC(I,J)
  898.                RX2(I,J) = RX(I,J) + V1(I)*V1(J)
  899.                r_z = r_z + RX2(I,J)*PHIR(J)
  900.              ENDDO
  901.              res2(I) = r_z
  902.            ENDDO
  903. *          c) Rw'*Phi
  904.            CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)
  905.            IF (BAL.EQ.1) THEN
  906.              DO I=1,NT
  907.                Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  908.              ENDDO
  909.            END IF
  910.            res3 = DDOT(NT,Rw,1,PHIR,1)
  911. *          d) norm(Phi) = 1
  912.            res4 = DDOT(NT,PHIR,1,PHIR,1)-ONE
  913. *          Y = ( R  ;  RQ * \phi  ;  Rw\phi  ;  \phi*\phi-1 )
  914.            arrRes2(1) = dnrm2(NT,fvec0.prog(1),1)
  915.            arrRes2(2) = dnrm2(NT,res2,1)
  916.            arrRes2(3) = abs(res3)
  917.            arrRes2(4) = abs(res4)
  918.            res = MAXVAL(arrRes2,1)
  919.            IF (IIMPI.GE.3)
  920.      &     WRITE(IOIMP,992) its,(arrRes2(iou),iou=1,4),res
  921.  992       FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';',
  922.      &     E10.3,';',E10.3,'} -> res=',E10.3)
  923.          ENDIF
  924.  
  925. *-----------------------------------------------------------------------
  926.          IF (TYPC.EQ.'P') THEN
  927. *-----------------------------------------------------------------------
  928.            CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)
  929.            DO I = 1,NT
  930.              DO J = 1,NT
  931.                RX(I,J) = ZZ(I,J) - JAC(I,J)
  932.              ENDDO
  933.            ENDDO
  934.            CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D)
  935.            CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1)
  936. *          b) (Rx-((OMEG/2)**2)*D2)*PhiR - KAPPA*D1*PhiI
  937. *          c) (Rx-((OMEG/2)**2)*D2)*PhiI + KAPPA*D1*PhiR
  938.            DO I = 1,NT
  939.              r_z1 = ZERO
  940.              r_z2 = ZERO
  941.              DO J = 1,NT
  942.                r_z1 = r_z1 + Rxmk2D(I,J)*PHIR(J)
  943.      &                     - KAPPA*DEL1(I,J)*PHII(J)
  944.                r_z2 = r_z2 + Rxmk2D(I,J)*PHII(J)
  945.      &                     + KAPPA*DEL1(I,J)*PHIR(J)
  946.              ENDDO
  947.              resV1(I) = r_z1
  948.              resV2(I) = r_z2
  949.            ENDDO
  950. *          d) (VQ'*PHIR) - 1
  951.            res4 = DDOT(NT,VQ,1,PHIR,1)-ONE
  952. *          Y = ( R  ;  [D0-k^2*D2]*\phi_R - k*D1*\phi_I  ;
  953. *                      [D0-k^2*D2]*\phi_I + k*D1*\phi_R  ; Q*\phi_I-1 )
  954.            arrRes2(1) = dnrm2(NT,fvec0.prog(1),1)
  955.            arrRes2(2) = dnrm2(NT,resV1,1)
  956.            arrRes2(3) = dnrm2(NT,resV2,1)
  957.            arrRes2(4) = abs(res4)
  958.            res = MAXVAL(arrRes2,1)
  959.            IF (IIMPI.GE.3)
  960.      &     WRITE(IOIMP,993) its,(arrRes2(iou),iou=1,4),res
  961.  993       FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';',
  962.      &     E10.3,';',E10.3,'} -> res=',E10.3)
  963.          ENDIF
  964.  
  965. *-----------------------------------------------------------------------
  966.          IF (TYPC.EQ.'N') THEN
  967. *-----------------------------------------------------------------------
  968.            CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)
  969.            DO J = 1,NT
  970.              DO I = 1,NT
  971.                RX(I,J) = ZZ(I,J) - JAC(I,J)
  972.              ENDDO
  973.            ENDDO
  974.            CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D)
  975.            CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1)
  976. *          b) (Rx-(KAPPA**2)*D2)*PhiR - KAPPA*D1*PhiI
  977. *          c) (Rx-(KAPPA**2)*D2)*PhiI + KAPPA*D1*PhiR
  978.            DO I = 1,NT
  979.              r_z1 = ZERO
  980.              r_z2 = ZERO
  981.              DO J = 1,NT
  982.                r_z1 = r_z1 + Rxmk2D(I,J)*PHIR(J)
  983.      &                     - KAPPA*DEL1(I,J)*PHII(J)
  984.                r_z2 = r_z2 + Rxmk2D(I,J)*PHII(J)
  985.      &                     + KAPPA*DEL1(I,J)*PHIR(J)
  986.              ENDDO
  987.              resV1(I) = r_z1
  988.              resV2(I) = r_z2
  989.            ENDDO
  990. *          d) (VQ'*PHIR) - 1
  991.            res4 = DDOT(NT,VQ,1,PHIR,1)-ONE
  992. *          e) VQ'*PHII
  993.            res5 = DDOT(NT,VQ,1,PHII,1)
  994. *          Y = ( R  ;  [D0-k^2*D2]*\phi_R - k*D1*\phi_I  ;
  995. *                      [D0-k^2*D2]*\phi_I + k*D1*\phi_R  ;
  996. *                                            vQ*\phi_R-1 ; vQ*\phi_I)
  997.            arrRes3(1) = dnrm2(NT,fvec0.prog(1),1)
  998.            arrRes3(2) = dnrm2(NT,resV1,1)
  999.            arrRes3(3) = dnrm2(NT,resV2,1)
  1000.            arrRes3(4) = abs(res4)
  1001.            arrRes3(5) = abs(res5)
  1002.            res = MAXVAL(arrRes3,1)
  1003.            IF (IIMPI.GE.3)
  1004.      &     WRITE(IOIMP,994) its,(arrRes3(iou),iou=1,5),res
  1005. c      &     WRITE(IOIMP,994) its,KAPPA,(arrRes3(iou),iou=1,4),res
  1006.  994       FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';',
  1007.      &     E10.3,';',E10.3,';',E10.3,'} -> res=',E10.3)
  1008.          ENDIF
  1009. *-----------------------------------------------------------------------
  1010. c            IF (IIMPI.GE.3) WRITE(IOIMP,*) 'Iter ',its,' -> res =',res
  1011.            IF (res.LE.TOLMIN) THEN
  1012.              CHECK=.false.
  1013.              CBIF = CBIF + 1
  1014.              TYPBIF(CBIF) = TYPC
  1015.              DO I=1,NT
  1016.                QBIFU(I,CBIF) = Q1(I)
  1017.                QPSIR(I,CBIF) = PHIR(I)
  1018.                QPSII(I,CBIF) = PHII(I)
  1019.                Q1(NT) = QIN(I)
  1020.              ENDDO
  1021.              WBIFU(CBIF) = OMEG
  1022.              WBIF2(CBIF) = KAPPA
  1023.              OMEG = OMEGIN
  1024. c              RETURN
  1025.              GOTO 666
  1026.            ELSE
  1027.              CHECK = .true.
  1028.            ENDIF
  1029.        ENDDO
  1030.        DO I = 1,NT
  1031.          Q1(I) = QIN(I)
  1032.        ENDDO
  1033.        OMEG = OMEGIN
  1034. c      message
  1035.        IF (CHECK) THEN
  1036. c          Pas de convergence apres %i1 iterations. L'execution continue
  1037.            INTERR(1)=MAXITS
  1038.            CALL ERREUR(151)
  1039. c            RETURN
  1040.              GOTO 666
  1041.        ELSEIF (IIMPI.GE.2) THEN
  1042.            WRITE(IOIMP,*) '--> convergence apres ',its,'iterations'
  1043.        ENDIF
  1044.  
  1045. *   MENAGE
  1046.  666  CONTINUE
  1047.       SEGSUP,MWORK
  1048.       SEGSUP,fvec,fvec0
  1049.  
  1050.       RETURN
  1051.       END
  1052.  
  1053.  
  1054.  

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