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    CB215821  23/01/25    21:15:18     11573          
  2.  
  3.       SUBROUTINE HBMBIF(NT,NHBM,NDDL,NFFT,KTQ,KTKAM,MTPHI,KTEMP,MTLIAA,
  4.      & MTLIAB,MTFEX,MTPAS,LOCLFA,LOCLB1,CHECK,TYPC,KPARNUM,KSORT)
  5. *
  6. *=======================================================================  
  7. * Calcul des points de bifurcation par une methode de Newton-Raphson
  8. *=======================================================================  
  9. * TYPC indique le type de calcul a realiser
  10. *      = 'L' pour un point limite
  11. *      = 'B' pour un point de branchement
  12. *      = 'P' pour un doublement de periode
  13. *      = 'N' pour une bifurcation de Neimark-Sacker
  14. *=======================================================================     
  15.         IMPLICIT INTEGER(I-N)
  16.         IMPLICIT REAL*8(A-H,O-Z)
  17.  
  18.         CHARACTER*1 TYPC
  19.         LOGICAL CHECK,CANAL
  20.         INTEGER MAXITS,MAXVCP,INFO
  21.  
  22. *     segment local (a nettoyer et faire remonter dans l'include TMDYNC.INC + tard)
  23. c       SEGMENT MWORK        
  24.         INTEGER IPIV33(2*NT+2),INDEIG(2*NT),I4(3*NT+2),I5(2*NT+1)
  25.         REAL*8 fvec(NT),AVVP(NT),AVVP2(NT+1),V1(NT)
  26.         REAL*8 EIGRE(NT),EIGIM(NT),QIN(NT),OMEGIN,KAPPA,KAPPA2
  27.         REAL*8 EIGRE2(NT+1),EIGIM2(NT+1)
  28.         REAL*8 VR(NT,NT),VL(NT,NT),WORK(8*NT)
  29.         REAL*8 VR2(NT+1,NT+1),VL2(NT+1,NT+1),WORK2(8*(NT+1))
  30.         REAL*8 JAC1(NT,NT),JAC2(NT,NT),fvec0(NT)
  31.         REAL*8 dZw(NT,NT),RxwPhi(NT),RXXPHI(NT,NT),RwwPhi
  32.         REAL*8 vA1(NT), vA2(NT),vA3(NT),vB1(NT),vB2(NT),vB3(NT)
  33.         REAL*8 RX2(NT,NT),SS,EI(NT)
  34.         REAL*8 AA(3,3),sol(3),arrRes(3),arrRes2(4),arrRes3(5)
  35.         REAL*8 JBP(2*NT+2,2*NT+2),solBP(2*NT+2)
  36.         REAL*8 dWB, dXB(NT),dPHIR(NT),dPHII(NT),WORK2w(16*NT)
  37.         REAL*8 res2(NT),Rco(NT,NT),PHIR(NT),PHII(NT),DPsiPsi(NT,NT)
  38.         REAL*8 MREG(NT+1,NT+1),DEL1(NT,NT),DEL2(NT),Mi,Ci,JH(2*NT,2*NT)
  39.         REAL*8 EXRE2w(2*NT),EXIM2w(2*NT),VL2w(2*NT,2*NT),VR2w(2*NT,4*NT)
  40.         REAL*8 dRxDdk1(NT),dRxDdk2(NT)
  41.         REAL*8 Rxmk2D(NT,NT),D10(NT,NT),D1w(NT,NT),D2(NT,NT),VQ(NT)
  42.         REAL*8 RxxPHIR(NT,NT),RxxPHII(NT,NT),RxwPhiR(NT),RxwPhiI(NT)
  43.         REAL*8 JNS(3*NT+2,3*NT+2),solNS(3*NT+2),resV1(NT),resV2(NT)
  44.         REAL*8 JPD(3*NT+1,3*NT+1),solPD(3*NT+1),D1wP(NT,NT)
  45. c       ENDSEGMENT
  46.  
  47. *-----------------------------------------------------------------------        
  48. * BIF: contient les coefficients de Fourier au point de bifurcation, 
  49. * ainsi que les deux fréquences correspondantes (KAPPA = 0. pour LP,BP)              
  50. * et les parties reelle/imaginaire du vecteur propre associe. 
  51. *-----------------------------------------------------------------------
  52. -INC PPARAM
  53. -INC CCOPTIO
  54. *-INC TMDYNC.INC
  55. ************************** debut TMDYNC.INC ****************************
  56.  
  57. *    TMDYNC : FUTUR INCLUDE POUR LES SEGMENTS DE L'OPERATEUR DYNC
  58. *    TODO   : a extraire dans un include des que stabilise
  59. *
  60. *     Segment des variables generalisees:
  61. *     -----------------------------------
  62.       SEGMENT MTQ
  63.          REAL*8 Q1(NT1)
  64.          REAL*8 OMEG,XPARA
  65.          REAL*8 JAC(NT1,NT1),ZZ(NT1,NT1),RX(NT1,NT1)
  66.          REAL*8 dX(NT1), dw, dv
  67.       ENDSEGMENT
  68. *     Q1   : vecteur des inconnues frequentielles de dimension (2h+1)*n
  69. *            Q1 = {q_0 q_c1 q_s1 ... q_sh}
  70. *            avec q_i vecteur de dimension n ou n=nombre de modes
  71. *     OMEG : frequence fondamentale de l'approximation
  72. *     XPARA: parametre de continuation (par defaut la frequence) 
  73. *            \in [PARINI,PARFIN]
  74. *     RX   : matrice jacobienne = ZZ + dFnl/dX
  75. *     JAC  : jacobienne des efforts non-lineaires  = dFnl/dX
  76. *     ZZ   : matrice dynamique associee aux matrices modales K, M et C 
  77. *            lineaires et constantes
  78. *     {dX,dw,(dv)} : vecteur tangent utilise pour la prediction
  79. *
  80. *
  81. *     Segment contenant les matrices XK, XASM et XM:
  82. *     ---------------------------------------------
  83.       SEGMENT MTKAM
  84.          REAL*8 XK(NA1,NB1K),XASM(NA1,NB1C),XM(NA1,NB1M)
  85.          REAL*8 GAM(NPC1,nl1),IGAM(nl1,NPC1),DL(nl1)
  86. *         REAL*8 GAMFIN(NPC2,nl1)
  87.       ENDSEGMENT
  88. *     XK,XASM et XM : matrices de raideur, amortissement et masse
  89. *     GAM et IGAM : matrices pour la FFT et son inverse
  90. *     GAMFIN      : 
  91. *
  92. *     Segment des deformees modales:
  93. *     ------------------------------
  94. *     (idem DYNE)
  95.       SEGMENT MTPHI
  96.          INTEGER IBASB(NPLB),IPLSB(NPLB),INMSB(NSB),IORSB(NSB)
  97.          INTEGER IAROTA(NSB)
  98.          REAL*8  XPHILB(NSB,NPLSB,NA2,IDIMB)
  99.       ENDSEGMENT
  100. *
  101. *     Segment descriptif des liaisons en base A:
  102. *     ------------------------------------------
  103. *     (idem DYNE)
  104.       SEGMENT MTLIAA
  105.          INTEGER IPALA(NLIAA,NIPALA),IPLIA(NLIAA,NPLAA),JPLIA(NPLA)
  106.          REAL*8  XPALA(NLIAA,NXPALA)
  107.       ENDSEGMENT
  108. *
  109. *     Segment descriptif des liaisons en base B:
  110. *     ------------------------------------------
  111. *     (idem DYNE)
  112.       SEGMENT MTLIAB
  113.          INTEGER IPALB(NLIAB,NIPALB),IPLIB(NLIAB,NPLBB),JPLIB(NPLB)
  114.          REAL*8  XPALB(NLIAB,NXPALB)
  115.          REAL*8  XABSCI(NLIAB,NIP),XORDON(NLIAB,NIP)
  116.       ENDSEGMENT
  117. *
  118. *     Segment representant les chargements exterieurs:
  119. *     -----------------------------------------------
  120.       SEGMENT MTFEX
  121.          REAL*8  FEXA(NT1)
  122.          REAL*8  FEXPSM(NPLB,NPC1,2,IDIMB)
  123.          INTEGER BAL
  124.       ENDSEGMENT
  125. *     FEXA : Vecteur des efforts ext. sous la forme de coefficients de
  126. *            Fourier et exprimes en base A
  127. *     FEXPSM: chargement/deplacement statique lie aux modes negliges
  128. *             (neglige aussi les Fnl). Dans DYNC toujours =0, cree pour
  129. *             compatibilite avec calcul des Fnl.
  130. *     BAL  : indique s'il s'agit d'un chargement de type balourd
  131. *            (cad proportionnel a OMEG**2)
  132. *
  133. *     Segment "local" pour DEVLFA:
  134. *     ----------------------------
  135.       SEGMENT LOCLFA
  136.          REAL*8 FTEST(NA1,4)
  137.       ENDSEGMENT
  138. *      
  139. *     Segment "local" pour DEVLB1:
  140. *     ----------------------------
  141.       SEGMENT LOCLB1
  142.          REAL*8 FTEST2(NPLB,6)
  143.       ENDSEGMENT
  144. *      
  145. *     Segment contenant les variables au cours d un pas de temps:
  146. *     ----------------------------------------------------------
  147.       SEGMENT MTPAS
  148.          REAL*8 FTOTA(NA1,4),FTOTB(NPLB,IDIMB),FTOTBA(NA1)
  149.          REAL*8 XPTB(NPLB,2,IDIMB),FINERT(NA1,4)
  150.          REAL*8 XVALA(NLIAA,4,NTVAR),XVALB(NLIAB,4,NTVAR)
  151.          REAL*8 FEXB(NPLB,2,IDIM),XCHPFB(2,NLIAB,4,NPLB)
  152.          REAL*8 KTOTXA(NA1,NA1),KTOTVA(NA1,NA1)
  153.          REAL*8 KTOTXB(NPLB,IDIMB,IDIMB), KTOTVB(NPLB,IDIMB,IDIMB) 
  154.       ENDSEGMENT
  155. *     FTOTA/B/BA : forces sur base A, B et B projetees sur A
  156. *     XPTB       : deplacement du point d'une liaison en base B
  157. *     XVALA/B    : grandeurs de la liaison en base A/B a stocker
  158. *     FEXB       : forces exterieures en base B (a priori uniquement 
  159. *                  pour les moments appliques aux rotations rigides ?)
  160. *     XCHPFB     : forces de contact en base B (lorsqu'on  considere un
  161. *                  maillage de contact dans certaines liaisons)
  162. *     KTOTXA/XB/VA/VB : Jacobienne par rapport au deplacement/vitesse en
  163. *                       base A/B (= contributions a dFnl/dX)
  164. *
  165. *
  166. *     Segment des points de reference des modes (base A):
  167. *     --------------------------------------------------
  168.       SEGMENT MPREF
  169.          INTEGER IPOREF(NPREF)
  170.       ENDSEGMENT
  171. *
  172. *     Segment des points en base B:
  173. *     -----------------------------
  174.       SEGMENT NCPR(NBPTS)
  175. *     NCRP(#global) = #local dans XPTB (1er indice)      
  176. *
  177. *     Segment des parametres numeriques pour la continuation: 
  178. *     ------------------------------------------------------
  179.       SEGMENT PARNUM
  180.          CHARACTER*4 TYPS
  181.          REAL*8   DS,DSMAX,DSMIN,ANGMIN,ANGMAX,ITERMOY,ISENS,TOLMIN
  182.          REAL*8   PARINI,PARFIN
  183.          INTEGER  ITERMAX,NBPAS
  184.          LOGICAL  JANAL
  185.       ENDSEGMENT
  186. *
  187. *     Segment des resultats:
  188. *     ---------------------
  189.       SEGMENT  PSORT
  190.         REAL*8 QSAVE(NT1,NPAS),WSAVE(NPAS),LSAVE(2,2*NA1,NPAS)
  191.         REAL*8 VSAVE(NPAS)
  192.         LOGICAL ZSAVE(NPAS)
  193.         CHARACTER*2 TYPBIF(NBIFU)
  194.         REAL*8 QBIFU(NT1,NBIFU),WBIFU(NBIFU),WBIF2(NBIFU)
  195.         REAL*8 QPSIR(NT1,NBIFU),QPSII(NT1,NBIFU)
  196.         INTEGER CBIF
  197.       ENDSEGMENT
  198. *     QSAVE(i,j) = Q harmonique i au pas j
  199. *     VSAVE(j)   = parametre de continuation (si non w) au j-eme pas
  200. *     ZSAVE(j)   = stabilite au j-eme pas
  201. *     LSAVE(1,j) : partie reelle de l'exposant de Floquet
  202. *     LSAVE(2,j) : partie imaginaire de l'exposant de Floquet
  203. *     TYPBIF = {LimitPoint, BranchPoint, NeimarkSacker, PeriodDoubling}
  204. *     QBIFU,WBIFU : vecteur Q et w au point de bifurcation
  205. *           WBIF2 : partie imaginaire de l'exposant de Floquet
  206. *     QPSIR,QPSII : vecteur propre au point de bifurcation
  207.  
  208. *     Segment des tableaux de travail:
  209. *     -------------------------------
  210.       SEGMENT  MTEMP
  211.          REAL*8 RW(NT1),A,T0(NT1+1),TP(NT1+1),AMPX,AUX
  212.          REAL*8 T02(NT1+2), TP2(NT1+2)
  213.          INTEGER IPIV(NT1),IPIV2(NT1+1),IPIV3(NT1+2)
  214.          REAL*8 res
  215.          REAL*8 RHS(NT1+1),Ja(NT1+1,NT1+1)
  216.          REAL*8 QOLD(NT1),OMEGOLD
  217.          REAL*8 MATJA(NT1+1,NT1+1),Rw2(NT1)
  218.          REAL*8 Jaa(NT1+2,NT1+2),RHS2(NT1+2),Ra(NT1),VV,VVOLD
  219.       ENDSEGMENT
  220. *     Jacobiennes augmentees 
  221. *     Ja : [ RX  Rw ; dX dw]
  222. *     Jaa: [ RX  Rw Ra; gx 0 0; dX dw da]
  223.  
  224. *      SEGMENT  NNNN
  225. *         REAL*8 IGAM2(nl1,NPC2),DL2(nl1)
  226. *      ENDSEGMENT
  227.  
  228. *************************** fin TMDYNC.INC *****************************
  229.  
  230.        REAL*8 ZERO,ONE,TWO,XEPS
  231.        PARAMETER (ZERO=0.D0, ONE=1.D0, TWO=2.D0, XEPS=1.D-6)
  232.  
  233. C Fonctions BLAS/LAPACK
  234.        REAL*8    DDOT, DNRM2
  235.        EXTERNAL  DDOT, DNRM2
  236.  
  237. *-----------------------------------------------------------------------
  238.        PARAMETER(MAXITS=40)
  239. *       TOLF=1.e-4
  240. *       TOLMIN=1.e-6
  241. *       STPMX=100.  
  242. *
  243.        MTQ   = KTQ
  244.        MTKAM = KTKAM
  245.        MTEMP = KTEMP
  246.        PSORT = KSORT
  247.        PARNUM = KPARNUM
  248.        CANAL = JANAL
  249.  
  250. c        SEGINI,MWORK
  251.  
  252. *      sauvegarde du point de la courbe de reponse
  253.        OMEGIN = OMEG
  254.        DO I = 1,NT 
  255.          QIN(I) = Q1(I)
  256.        ENDDO       
  257. *        
  258.        CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  259.      &            MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec0,fref0,
  260.      &            .false.)
  261.        CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)  
  262.  
  263. c      ----------------------------------------------------------------       
  264. c      !!! Localisation NeimarkSacker et PeriodDoubling a debugger !!!
  265. c       if(TYPC.eq.'N'.or.TYPC.eq.'P') then
  266. c          WRITE(IOIMP,*) 'Localisation bifurcation ',TYPC,'non prevue !'
  267. c          CHECK = .true.
  268. c          return
  269. c       endif
  270. c      ---------------------------------------------------------------- 
  271.  
  272. c      message
  273.        IF (IIMPI.GE.2) WRITE(IOIMP,*) 'Localisation bifurcation ',TYPC
  274.        IF (IIMPI.GE.3) WRITE(IOIMP,990) 0,dnrm2(NT,fvec0,1),fref0
  275.  990   FORMAT('+ hbmbif : iter ',I2,' -> |R|=',E10.3,' ref=',E10.3) 
  276.  
  277.  
  278. *=======================================================================     
  279. *                            Initialisation                              
  280. *======================================================================= 
  281.  
  282. * La Jacobienne augmentee a la structure de base suivante: 
  283. *   JJ = [    Rx      O      Rw  ]
  284. *        [ (RxPhi)x   Rx   RxwPhi]
  285. *        [   0'     2Phi'    0   ]
  286. * qui est a adapter selon la bifurcation traitee.
  287. *
  288. * Il faut initialiser le vecteur propre Phi.
  289. * Calcul de la Jacobienne Rx = Z(OMEGA) - dFNLdX
  290. c        DO I=1,NT
  291.        DO J=1,NT
  292.        DO I=1,NT
  293.            RX(I,J) = ZZ(I,J)-JAC(I,J)
  294.        ENDDO
  295.        ENDDO
  296.        CALL COPYMAT(NT,Rx,Rco)
  297.  
  298.  
  299. *-----------------------------------------------------------------------               
  300. *      Bifurcation statique: le vecteur propre annule la Jacobienne
  301. *----------------------------------------------------------------------- 
  302.  
  303.        IF (TYPC.EQ.'L') THEN
  304.  
  305.          CALL DGEEV('N','V',NT,Rco,NT,EIGRE,EIGIM,VL,1,VR,NT,
  306.      &                                                   WORK,8*NT,INFO)
  307. *        On prend le vecteur associe a la partie reelle la plus petite.     
  308.          MINVP=1
  309.          XMINVP=ABS(EIGRE(1))
  310.          DO I = 1,NT
  311.            AVVP(I) = ABS(EIGRE(I))
  312.            IF(AVVP(I).LT.XMINVP) THEN
  313.              MINVP=I
  314.              XMINVP=AVVP(I)
  315.            ENDIF
  316.          ENDDO           
  317. c          MINVP = MINLOC(AVVP,NT)
  318. *        Les vecteurs propres sont purement reels dans ce cas.         
  319.          DO I = 1,NT
  320.             PHIR(I) = VR(I,MINVP)
  321.             PHII(I) = ZERO
  322.          ENDDO
  323. *        Pas de nouvelle frequence pour les cycles bifurques          
  324.          KAPPA = ZERO
  325.          MAXVCP=1
  326.          XMAXVCP=ABS(PHIR(1))
  327.          DO I = 1,NT
  328.            AVVP(I) = ABS(PHIR(I))
  329.            IF(AVVP(I).GT.XMAXVCP) THEN
  330.              MAXVCP=I
  331.              XMAXVCP=AVVP(I)
  332.            ENDIF
  333.          ENDDO
  334. c          MAXVCP = MAXLOC(AVVP,NT)
  335. *        Elements pour la resolution par blocs          
  336.          SS = ZERO
  337.          DO I=1,NT
  338.            SS = SS + abs(RX(I,I))
  339.          ENDDO         
  340.          SS = SS/NT
  341. *        Penalisation de la Jacobienne (evite mauvais conditionnement)
  342.          CALL COPYMAT(NT,RX,RX2)
  343.          DO I=1,NT
  344.           EI(I) = ZERO
  345.          ENDDO
  346.          RX2(MAXVCP,MAXVCP) = RX2(MAXVCP,MAXVCP) + SS
  347.          EI(MAXVCP) = ONE
  348.        ENDIF
  349.  
  350. *-----------------------------------------------------------------------  
  351. *      Branch point ?
  352. *-----------------------------------------------------------------------               
  353.        IF (TYPC.EQ.'B') THEN
  354.  
  355.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)  
  356.          IF (BAL.EQ.1) THEN
  357.            DO I=1,NT
  358.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  359.            ENDDO
  360.          END IF 
  361. *        Regularisation de la Jacobienne
  362.          CALL HBMDVEC(NT,NHBM,NDDL,Q1,ONE,V1)         
  363.          DO I = 1,NT
  364.            MREG(I,NT+1) = Rw(I)
  365.            MREG(NT+1,I) = Rw(I)
  366.          ENDDO
  367.          znormV1 = dnrm2(NT,V1,1)
  368.          DO I = 1,NT 
  369.             V1(I) = V1(I)/znormV1
  370.          ENDDO
  371.          DO I = 1,NT
  372.            DO J = 1,NT
  373.              RX2(I,J) = RX(I,J) + V1(I)*V1(J)  
  374.              MREG(I,J) = RX2(I,J)
  375.            ENDDO
  376.          ENDDO    
  377.          MREG(NT+1,NT+1) = ZERO
  378.          CALL DGEEV('N','V',NT+1,MREG,NT+1,EIGRE2,EIGIM2,VL2,1,VR2,NT+1,
  379.      &                                              WORK2,8*(NT+1),INFO)         
  380. *        On prend le vecteur associe a partie reelle la plus petite.     
  381.          MINVP=1
  382.          XMINVP=ABS(EIGRE2(1))
  383.          DO I = 1,NT+1
  384.            AVVP2(I) = ABS(EIGRE2(I))
  385.            IF(AVVP2(I).LT.XMINVP) THEN
  386.              MINVP=I
  387.              XMINVP=AVVP2(I)
  388.            ENDIF
  389.          ENDDO           
  390. c          MINVP = MINLOC(AVVP2,NT+1)
  391. *        Les vecteurs propres sont purement reels dans ce cas.         
  392.          DO I = 1,NT
  393.             PHIR(I) = VR2(I,MINVP)
  394.             PHII(I) = ZERO
  395.          ENDDO
  396. *        Pas de nouvelle frequence pour les cycles bifurques          
  397.          KAPPA = ZERO                 
  398.        ENDIF
  399.  
  400. *-----------------------------------------------------------------------               
  401. *      Bifurcations dynamiques (Period doubling ou Neimark Sacker)
  402. *-----------------------------------------------------------------------               
  403. *      Construction de la matrice de Hill, JJ
  404.        IF ((TYPC.EQ.'P').OR.(TYPC.EQ.'N')) THEN
  405.          DO I = 1,NT
  406.             DO J = 1,NT
  407.               D10(I,J)=ZERO
  408.               D1w(I,J)=ZERO    
  409.               D2(I,J) =ZERO
  410.              ENDDO
  411.          ENDDO
  412.          DO I=1,2*NHBM+1
  413.            DO J=1,NDDL
  414.               DEL2(NDDL*(I-1)+J) = ONE/XM(J,1)
  415.               D2(NDDL*(I-1)+J,NDDL*(I-1)+J) = XM(J,1)
  416.            ENDDO
  417.          ENDDO
  418.          DO I=1,NDDL
  419.              D10(I,I) = XASM(I,1)      
  420.          ENDDO       
  421.          DO J = 2,2*NHBM,2
  422.            DO I=1,NDDL
  423.                D10(NDDL*(1+(J-2))+I,NDDL*(1+(J-2))+I) = XASM(I,1)
  424.                D1w(NDDL*(1+(J-2))+I,NDDL*(1+(J-1))+I) = J*XM(I,1)
  425.                D1w(NDDL*(1+(J-1))+I,NDDL*(1+(J-2))+I) = -J*XM(I,1)
  426.                D10(NDDL*(1+(J-1))+I,NDDL*(1+(J-1))+I) = XASM(I,1)
  427.            ENDDO             
  428.          ENDDO
  429. *        DEL1 = DEL10 + OMEG*DEL1w
  430.          CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1)
  431. *     JJ = [DEL2\DEL1   DEL2\Rx]
  432. *          [     -I         0  ] 
  433. * 
  434.          DO  I=1,NT
  435.            DO  J=1,NT
  436. *            JJ_11 = DEL2\DEL1
  437.              JH(I,J) = -DEL1(I,J)*DEL2(J) 
  438. *            JJ_12 = DEL2\Rx 
  439.              JH(I,NT+J) = -RX(I,J)*DEL2(J)           
  440. *            JJ_21 = [I]           
  441.              IF (I.EQ.J) THEN
  442.                JH(NT+I,J) = ONE
  443.              ELSE
  444.                JH(NT+I,J) = ZERO
  445.              ENDIF
  446. *            JJ_22 = [0]
  447.              JH(NT+I,NT+J) = ZERO         
  448.            ENDDO
  449.          ENDDO
  450.          CALL DGEEV('N','V',2*NT,JH,2*NT,EXRE2w,EXIM2w,VL2w,1,VR2w,2*NT,
  451.      &                                                 WORK2,16*NT,INFO)
  452.          CALL HBMORDO(2*NT,NDDL,OMEG,EXIM2w,INDEIG)
  453. *         WRITE(*,*) '   REAL   ', '   IMAG   ', '   INDX   '
  454. *         DO I = 1,2*NT
  455. *           WRITE(*,*) EXRE2w(INDEIG(I)),EXIM2w(I),INDEIG(I)
  456. *         ENDDO         
  457. *        WRITE(*,*) 'INDXS:',(INDEIG(IOU),IOU=1,2*NDDL)
  458. *
  459. *        On choisit le vecteur associe a la valeur propre avec partie
  460. *        imaginaire positive du mode instable
  461.          MINVP=1
  462.          DO I=1,2*NDDL
  463.            IF ((EXRE2w(INDEIG(I)).GT.ZERO).AND.(EXIM2w(I).GT.ZERO)) THEN
  464.               MINVP=INDEIG(I)
  465.               INDK =I
  466. *             WRITE(*,*) 'Mode instable:',MINVP
  467.             ENDIF
  468.          ENDDO
  469.          DO I=1,NT 
  470.            PHIR(I) = VR2w(I,MINVP)
  471.            PHII(I) = VR2w(I,MINVP+1)  
  472.          ENDDO  
  473. *         WRITE(*,*) 'PHI_R  ', '  PHI_I'
  474. *         DO IOU = 1,NT 
  475. *           WRITE(*,*) PHIR(IOU),PHII(IOU)
  476. *         ENDDO
  477. *         STOP
  478.          MAXVCP=1
  479.          XMAXVCP=ABS(PHII(1))
  480.          DO I = 1,NT
  481.            AVVP(I) = ABS(PHII(I))
  482.            IF (AVVP(I).GT.XMAXVCP) THEN
  483.              MAXVCP=I
  484.              XMAXVCP=AVVP(I)
  485.            ENDIF
  486.          ENDDO
  487. c          MAXVCP = MAXLOC(AVVP,NT)
  488.          DO I = 1,NT
  489.              VQ(I) = ZERO
  490.          ENDDO
  491.          VQ(MAXVCP) = ONE
  492. *        ---------------------         
  493.          IF (TYPC.EQ.'P') THEN
  494. *        Nouvelle frequence fixee a OMEGA/2     
  495.            KAPPA = OMEG/TWO
  496.          ENDIF              
  497. *        ---------------------
  498.          IF (TYPC.EQ.'N') THEN
  499. *        Nouvelle frequence incommensurable KAPPA   
  500.             KAPPA = ABS(EXIM2w(INDK))  
  501. *            write(*,*) 'KAPPA^0=',KAPPA
  502.          ENDIF
  503. *         STOP
  504.        ENDIF
  505. *=======================================================================     
  506. *                              Iterations                                   
  507. *=======================================================================     
  508.         DO its=1,MAXITS
  509. *-----------------------------------------------------------------------                   
  510. *        CAS L: POINT LIMITE 
  511. *-----------------------------------------------------------------------                   
  512.          IF (TYPC.EQ.'L') THEN 
  513. ***      Calcul des derivees
  514. *        RxxPhi
  515.          EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS)         
  516. *        Perturbation de la Jacobienne: +EPS         
  517.          DO JJ = 1,NT
  518.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  519.          ENDDO
  520.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  521.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  522.      &            .false.)
  523.          CALL COPYMAT(NT,JAC,JAC1)
  524. *        Perturbation de la Jacobienne: -EPS
  525.          DO JJ = 1,NT
  526.            Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ)
  527.          ENDDO
  528.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  529.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  530.      &            .false.)     
  531.          CALL COPYMAT(NT,JAC,JAC2)
  532. *        retour a la valeur initiale
  533.          DO JJ = 1,NT
  534.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  535.          ENDDO
  536. *        Differences centrees         
  537.          DO I = 1,NT
  538.            DO J = 1,NT
  539.              RxxPHI(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS1) 
  540.            ENDDO
  541.          ENDDO
  542. *        RxwPhi
  543.          CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,dZw)
  544.          DO I=1,NT
  545.            RxwPhi(I) = ZERO
  546.            DO J = 1,NT
  547.              RxwPhi(I) = RxwPhi(I) + dZw(I,J)*PHIR(J)
  548.            ENDDO 
  549.          ENDDO
  550. *        Rw
  551.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)  
  552.          IF (BAL.EQ.1) THEN
  553.            DO I=1,NT
  554.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  555.            ENDDO
  556.          END IF        
  557. *        Resolution par blocs 
  558. *        A1
  559.          DO I = 1,NT
  560.            vA1(I) = -fvec0(I)
  561.          ENDDO
  562.          CALL COPYMAT(NT,RX2,Rco)
  563.          CALL DGESV(NT,1,Rco,NT,IPIV,vA1,NT,INFO)
  564. *        A2
  565.          DO I = 1,NT
  566.            vA2(I) = -Rw(I)
  567.          ENDDO
  568.          CALL COPYMAT(NT,RX2,Rco)
  569.          CALL DGESV(NT,1,Rco,NT,IPIV,vA2,NT,INFO)
  570. *        A3 
  571.          DO I = 1,NT
  572.            vA3(I) = EI(I)
  573.          ENDDO
  574.          CALL COPYMAT(NT,RX2,Rco)
  575.          CALL DGESV(NT,1,Rco,NT,IPIV,vA3,NT,INFO)
  576. *        B1
  577.          DO I = 1,NT
  578.            vB1(I) = ZERO
  579.            DO J = 1,NT
  580.              vB1(I) = vB1(I) + RxxPHI(I,J)*vA1(J) + RX(I,J)*PHIR(J)
  581.            ENDDO
  582.            vB1(I) = -vB1(I) 
  583.          ENDDO
  584.          CALL COPYMAT(NT,RX2,Rco)
  585.          CALL DGESV(NT,1,Rco,NT,IPIV,vB1,NT,INFO)
  586. *        B2
  587.          DO I = 1,NT
  588.            vB2(I) = ZERO
  589.            DO J = 1,NT
  590.              vB2(I) = vB2(I) + RxxPHI(I,J)*vA2(J)
  591.            ENDDO
  592.            vB2(I) = -(vB2(I) + RxwPhi(I))
  593.          ENDDO
  594.          CALL COPYMAT(NT,RX2,Rco)
  595.          CALL DGESV(NT,1,Rco,NT,IPIV,vB2,NT,INFO)
  596. *        B3
  597.          DO I = 1,NT
  598.            vB3(I) = ZERO
  599.            DO J = 1,NT
  600.              vB3(I) = vB3(I) + RxxPHI(I,J)*vA3(J)
  601.            ENDDO
  602.            vB3(I) = -vB3(I)
  603.          ENDDO
  604.          CALL COPYMAT(NT,RX2,Rco)
  605.          CALL DGESV(NT,1,Rco,NT,IPIV,vB3,NT,INFO)         
  606. *        
  607. *        Construction de la matrice du probleme
  608.          AA(1,1) = vA2(MAXVCP)
  609.          AA(1,2) = SS*vA3(MAXVCP)-ONE
  610.          AA(1,3) = ZERO 
  611.          AA(2,1) = vB2(MAXVCP)
  612.          AA(2,2) = SS*vB3(MAXVCP) 
  613.          AA(2,3) = SS*vA3(MAXVCP)-ONE
  614.          AA(3,1) = DDOT(NT,PHIR,1,vB2,1)
  615.          AA(3,2) = SS*(DDOT(NT,PHIR,1,vB3,1))
  616.          AA(3,3) = SS*(DDOT(NT,PHIR,1,vA3,1))
  617. *        Construction du cote droit 
  618.          sol(1) = -vA1(MAXVCP)
  619.          sol(2) = -vB1(MAXVCP)
  620.          sol(3) = 0.5D0*(ONE-DDOT(NT,PHIR,1,PHIR,1))
  621.      &          -DDOT(NT,PHIR,1,vB1,1)
  622. *        Resolution
  623.          CALL DGESV(3,1,AA,3,IPIV,sol,3,INFO) 
  624. *        Corrections 
  625.          dWB = sol(1)
  626.          beta1 = sol(2) 
  627.          beta2 = sol(3)
  628.          DO I = 1,NT
  629.           dXB(I) = vA1(I) + dWB*vA2(I) +SS*beta1*vA3(I)
  630.           dPHIR(I)= vB1(I)+ dWB*vB2(I) +SS*beta1*vB3(I) +SS*beta2*vA3(I)          
  631.           Q1(I) = Q1(I) + dXB(I)
  632.           PHIR(I) = PHIR(I) + dPHIR(I)
  633.          ENDDO
  634.          OMEG = OMEG + dWB        
  635.          ENDIF
  636. *-----------------------------------------------------------------------        
  637. *        CAS B: POINT DE BRANCHEMENT
  638. *-----------------------------------------------------------------------                   
  639.          IF (TYPC.EQ.'B') THEN
  640. ***      Calcul des derivees
  641. *        RxxPhi
  642.          EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS)         
  643. *        Perturbation de la Jacobienne: +EPS         
  644.          DO JJ = 1,NT
  645.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  646.          ENDDO
  647.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  648.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  649.      &            .false.)     
  650.          CALL COPYMAT(NT,JAC,JAC1)
  651. *        Perturbation de la Jacobienne: -EPS
  652.          DO JJ = 1,NT
  653.            Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ)
  654.          ENDDO
  655.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  656.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  657.      &            .false.)     
  658.          CALL COPYMAT(NT,JAC,JAC2)
  659.          DO JJ = 1,NT
  660.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  661.          ENDDO
  662. c        DPsiPsi
  663.          CALL HBMDPP(NT,NDDL,Q1,PHIR,DPsiPsi)
  664. *        Differences centrees         
  665.          DO I = 1,NT
  666.            DO J = 1,NT
  667.            RxxPHI(I,J)=((JAC2(I,J)-JAC1(I,J))/(TWO*EPS1)) + DPsiPsi(I,J)
  668.            ENDDO
  669.          ENDDO
  670. *        RxwPhi
  671.          CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,dZw)
  672.          DO I=1,NT
  673.            RxwPhi(I) = ZERO
  674.            DO J = 1,NT
  675.              RxwPhi(I) = RxwPhi(I) + dZw(I,J)*PHIR(J)
  676.            ENDDO 
  677.          ENDDO
  678. *        Rw
  679.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)  
  680.          IF (BAL.EQ.1) THEN
  681.            DO I=1,NT
  682.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  683.            ENDDO
  684.          END IF      
  685. *        Rww*Phi
  686.          CALL HBMRWW(NT,NDDL,Q1,PHIR,XM,RwwPhi) 
  687.          IF (BAL.EQ.1) THEN
  688.            DO I=1,NT
  689.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  690.              RwwPhi = RwwPhi - TWO*FEXA(I)*PHIR(I)
  691.            ENDDO
  692.          END IF   
  693. *        Construction de la matrice du probleme et du cote droit
  694.          solBP(2*NT+1) = ZERO
  695.          DO I = 1,NT 
  696.            solBP(NT+I) = ZERO
  697.            DO J = 1,NT
  698.              JBP(I,J) = RX(I,J) 
  699.              JBP(NT+I,J) = RxxPHI(I,J) 
  700.              JBP(I,NT+J) = ZERO 
  701.              JBP(NT+I,NT+J) = RX2(I,J) 
  702.              solBP(NT+I) = solBP(NT+I) + RX2(I,J)*PHIR(J)
  703.            ENDDO 
  704.            JBP(2*NT+1,NT+I) = RxwPhi(I) 
  705.            JBP(2*NT+1,I) = Rw(I)
  706.            JBP(2*NT+2,I) = ZERO 
  707.            JBP(2*NT+2,NT+I) = TWO*PHIR(I) 
  708.            JBP(I,2*NT+1) = Rw(I) 
  709.            JBP(NT+I,2*NT+1) =  RxwPhi(I) 
  710.            JBP(I,2*NT+2) = PHIR(I) 
  711.            JBP(NT+I,2*NT+2) = ZERO
  712.            solBP(I) = -fvec0(I)
  713.            solBP(NT+I) = -solBP(NT+I) 
  714.          ENDDO
  715.          solBP(2*NT+1) = -(DDOT(NT,Rw,1,PHIR,1))
  716.          solBP(2*NT+2) = -(DDOT(NT,PHIR,1,PHIR,1)-ONE)
  717.          JBP(2*NT+1,2*NT+1) = RwwPhi
  718.          JBP(2*NT+1,2*NT+2) = ZERO
  719.          JBP(2*NT+2,2*NT+1) = ZERO
  720.          JBP(2*NT+2,2*NT+2) = ZERO
  721. *
  722. *        Resolution
  723.          CALL DGESV(2*NT+2,1,JBP,2*NT+2,IPIV33,solBP,2*NT+2,INFO) 
  724. *        Corrections 
  725.          dWB = solBP(2*NT+1)
  726.          DO I = 1,NT
  727.           dXB(I) = solBP(I)
  728.           dPHIR(I)= solBP(NT+I)          
  729.           Q1(I) = Q1(I) + dXB(I)
  730.           PHIR(I) = PHIR(I) + dPHIR(I)
  731.          ENDDO
  732.          OMEG = OMEG + dWB        
  733.          ENDIF
  734. *-----------------------------------------------------------------------        
  735. *        CAS P: DOUBLEMENT DE PERIODE 
  736. *-----------------------------------------------------------------------                   
  737.          IF (TYPC.EQ.'P') THEN
  738. ***      Calcul des derivees
  739. *        RxxPhiR
  740.          EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS)         
  741. *        Perturbation de la Jacobienne: +EPS1         
  742.          DO JJ = 1,NT
  743.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  744.          ENDDO
  745.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  746.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  747.      &            .false.)     
  748.          CALL COPYMAT(NT,JAC,JAC1)
  749. *        Perturbation de la Jacobienne: -EPS1
  750.          DO JJ = 1,NT
  751.            Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ)
  752.          ENDDO
  753.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  754.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  755.      &            .false.)     
  756.          CALL COPYMAT(NT,JAC,JAC2)
  757.          DO JJ = 1,NT
  758.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  759.          ENDDO
  760. *        Differences centrees         
  761.          DO I = 1,NT
  762.            DO J = 1,NT
  763.              RxxPHIR(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS1) 
  764.            ENDDO
  765.          ENDDO
  766. *        RxxPhiI
  767.          EPS2 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHII,1) + XEPS)         
  768. *        Perturbation de la Jacobienne: +EPS2         
  769.          DO JJ = 1,NT
  770.            Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ)
  771.          ENDDO
  772.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  773.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  774.      &            .false.)     
  775.          CALL COPYMAT(NT,JAC,JAC1)
  776. *        Perturbation de la Jacobienne: -EPS2
  777.          DO JJ = 1,NT
  778.            Q1(JJ) = Q1(JJ) - TWO*EPS2*PHII(JJ)
  779.          ENDDO
  780.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  781.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  782.      &            .false.)     
  783.          CALL COPYMAT(NT,JAC,JAC2)
  784.          DO JJ = 1,NT
  785.            Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ)
  786.          ENDDO
  787. *        Differences centrees         
  788.          DO I = 1,NT
  789.            DO J = 1,NT
  790.              RxxPHII(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS2) 
  791.            ENDDO
  792.          ENDDO
  793. *        Rw
  794.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)  
  795.          IF (BAL.EQ.1) THEN
  796.            DO I=1,NT
  797.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  798.            ENDDO
  799.          END IF  
  800. *        Calcul de: (Rxw-(OMEG/2)*D2)*PhiR-(0.5*D10+OMEG*D1w)*PhiI et
  801. *                   (Rxw-(OMEG/2)*D2)*PhiI+(0.5*D10+OMEG*D1w)*PhiR 
  802.          CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,dZw)
  803.          CALL COLIMA(NT,NT,DEL1,D10,ONE,-0.5D0,D1wP)
  804.          DO I=1,NT
  805.           RxwPhiR(I) = ZERO
  806.           RxwPhiI(I) = ZERO
  807.           DO J = 1,NT
  808.            RxwPhiR(I)=RxwPhiR(I)+(dZw(I,J)-KAPPA*D2(I,J))*PHIR(J)
  809.      &                                               - D1wP(I,J)*PHII(J)
  810.            RxwPhiI(I)=RxwPhiI(I)+(dZw(I,J)-KAPPA*D2(I,J))*PHII(J)
  811.      &                                               + D1wP(I,J)*PHIR(J)
  812.           ENDDO 
  813.          ENDDO
  814. *        Rxmk2D = Rx-(KAPPA**2)*D2
  815.          CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D)
  816. *        Construction de la matrice du probleme et du cote droit
  817.          DO I = 1,NT 
  818.            solPD(NT+I) = ZERO
  819.            solPD(2*NT+I) = ZERO
  820.            DO J = 1,NT
  821.              JPD(I,J) = RX(I,J) 
  822.              JPD(NT+I,J) = RxxPHIR(I,J)
  823.              JPD(2*NT+I,J) = RxxPHII(I,J)             
  824.              JPD(I,NT+J) = ZERO
  825.              JPD(I,2*NT+J) = ZERO             
  826.              JPD(NT+I,NT+J) = Rxmk2D(I,J)
  827.              JPD(NT+I,2*NT+J) = -KAPPA*DEL1(I,J)
  828.              JPD(2*NT+I,NT+J) = KAPPA*DEL1(I,J)
  829.              JPD(2*NT+I,2*NT+J) = Rxmk2D(I,J)    
  830.              solPD(NT+I)=solPD(NT+I)+Rxmk2D(I,J)*PHIR(J)
  831.      &                                          -KAPPA*DEL1(I,J)*PHII(J)             
  832.              solPD(2*NT+I)=solPD(2*NT+I)+Rxmk2D(I,J)*PHII(J)
  833.      &                                          +KAPPA*DEL1(I,J)*PHIR(J)
  834.            ENDDO
  835.            JPD(     I,3*NT+1) = Rw(I) 
  836.            JPD(  NT+I,3*NT+1) = RxwPhiR(I)
  837.            JPD(2*NT+I,3*NT+1) = RxwPhiI(I)
  838.            JPD(3*NT+1,I)      = ZERO
  839.            JPD(3*NT+1,2*NT+1) = VQ(I)
  840.            JPD(3*NT+1,NT+1)   = ZERO           
  841.            solPD(I)           = -fvec0(I)
  842.            solPD(NT+I)        = -solNS(NT+I)
  843.            solPD(2*NT+I)      = -solNS(2*NT+I)
  844.          ENDDO
  845.          solPD(3*NT+1) = -(DDOT(NT,VQ,1,PHII,1)-0.5D0)
  846.          JNS(3*NT+1,3*NT+1) = ZERO
  847. *        Resolution
  848.          CALL DGESV(3*NT+1,1,JPD,3*NT+1,I5,solPD,3*NT+1,INFO) 
  849. *        Corrections 
  850.          dWB = solPD(3*NT+1)
  851.          DO I = 1,NT
  852.           dXB(I) = solPD(I)
  853.           dPHIR(I)= solPD(NT+I)    
  854.           dPHII(I)= solPD(2*NT+I)    
  855.           Q1(I) = Q1(I) + dXB(I)
  856.           PHIR(I) = PHIR(I) + dPHIR(I)
  857.           PHII(I) = PHII(I) + dPHII(I)
  858.          ENDDO
  859.          OMEG = OMEG + dWB     
  860.          KAPPA= OMEG/TWO
  861.          END IF 
  862. *-----------------------------------------------------------------------        
  863. *        CAS N: BIFURCATION DE NEIMARK-SACKER 
  864. *-----------------------------------------------------------------------        
  865.          IF (TYPC.EQ.'N') THEN
  866. ***      Calcul des derivees
  867. *        RxxPhiR
  868.          EPS1 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHIR,1) + XEPS)         
  869. *        Perturbation de la Jacobienne: +EPS1         
  870.          DO JJ = 1,NT
  871.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  872.          ENDDO
  873.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  874.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  875.      &            .false.)     
  876.          CALL COPYMAT(NT,JAC,JAC1)
  877. *        Perturbation de la Jacobienne: -EPS1
  878.          DO JJ = 1,NT
  879.            Q1(JJ) = Q1(JJ) - TWO*EPS1*PHIR(JJ)
  880.          ENDDO
  881.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  882.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  883.      &            .false.)     
  884.          CALL COPYMAT(NT,JAC,JAC2)
  885.          DO JJ = 1,NT
  886.            Q1(JJ) = Q1(JJ) + EPS1*PHIR(JJ)
  887.          ENDDO
  888. *        Differences centrees         
  889.          DO I = 1,NT
  890.            DO J = 1,NT
  891.              RxxPHIR(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS1) 
  892.            ENDDO
  893.          ENDDO
  894. *        RxxPhiI
  895.          EPS2 = XEPS*(dnrm2(NT,Q1,1)/dnrm2(NT,PHII,1) + XEPS)         
  896. *        Perturbation de la Jacobienne: +EPS2         
  897.          DO JJ = 1,NT
  898.            Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ)
  899.          ENDDO
  900.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  901.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  902.      &            .false.)     
  903.          CALL COPYMAT(NT,JAC,JAC1)
  904. *        Perturbation de la Jacobienne: -EPS2
  905.          DO JJ = 1,NT
  906.            Q1(JJ) = Q1(JJ) - TWO*EPS2*PHII(JJ)
  907.          ENDDO
  908.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  909.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec,fref,
  910.      &            .false.)     
  911.          CALL COPYMAT(NT,JAC,JAC2)
  912.          DO JJ = 1,NT
  913.            Q1(JJ) = Q1(JJ) + EPS2*PHII(JJ)
  914.          ENDDO
  915. *        Differences centrees         
  916.          DO I = 1,NT
  917.            DO J = 1,NT
  918.              RxxPHII(I,J) = (JAC2(I,J)-JAC1(I,J))/(TWO*EPS2) 
  919.            ENDDO
  920.          ENDDO
  921. *        RxwPhiR-KAPPA*D1w*PhiI, RxwPhiI+KAPPA*D1w*PhiR 
  922.          CALL HBMZW(NT,NHBM,NDDL,MTQ,MTKAM,dZw)
  923.          DO I=1,NT
  924.           RxwPhiR(I) = ZERO
  925.           RxwPhiI(I) = ZERO
  926.           dRxDdk1(I) = ZERO
  927.           dRxDdk2(I) = ZERO
  928.           DO J = 1,NT
  929.            RxwPhiR(I)=RxwPhiR(I)+dZw(I,J)*PHIR(J)-KAPPA*D1w(I,J)*PHII(J)
  930.            RxwPhiI(I)=RxwPhiI(I)+dZw(I,J)*PHII(J)+KAPPA*D1w(I,J)*PHIR(J)
  931.            dRxDdk1(I)=dRxDdk1(I)-TWO*KAPPA*D2(I,J)*PHIR(J)
  932.      &                                                -DEL1(I,J)*PHII(J)
  933.            dRxDdk2(I)=dRxDdk2(I)-TWO*KAPPA*D2(I,J)*PHII(J)
  934.      &                                                +DEL1(I,J)*PHIR(J)
  935.           ENDDO 
  936.          ENDDO
  937. *        Rxmk2D = Rx-(KAPPA**2)*D2
  938.          CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D)
  939. *        Rw
  940.          CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)  
  941.          IF (BAL.EQ.1) THEN
  942.            DO I=1,NT
  943.              Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  944.            ENDDO
  945.          END IF   
  946. *        Construction de la matrice du probleme et du cote droit
  947.          DO I = 1,NT 
  948.            solNS(NT+I) = ZERO
  949.            solNS(2*NT+I) = ZERO
  950.            DO J = 1,NT
  951.              JNS(I,J) = RX(I,J) 
  952.              JNS(NT+I,J) = RxxPHIR(I,J)
  953.              JNS(2*NT+I,J) = RxxPHII(I,J)             
  954.              JNS(I,NT+J) = ZERO
  955.              JNS(I,2*NT+J) = ZERO             
  956.              JNS(NT+I,NT+J) = Rxmk2D(I,J)
  957.              JNS(NT+I,2*NT+J) = -KAPPA*DEL1(I,J)
  958.              JNS(2*NT+I,NT+J) = KAPPA*DEL1(I,J)
  959.              JNS(2*NT+I,2*NT+J) = Rxmk2D(I,J)    
  960.              solNS(NT+I)=solNS(NT+I)+Rxmk2D(I,J)*PHIR(J)
  961.      &                                          -KAPPA*DEL1(I,J)*PHII(J)             
  962.              solNS(2*NT+I)=solNS(2*NT+I)+Rxmk2D(I,J)*PHII(J)
  963.      &                                          +KAPPA*DEL1(I,J)*PHIR(J)
  964.            ENDDO
  965.            JNS(     I,3*NT+1) = ZERO 
  966.            JNS(  NT+I,3*NT+1) = dRxDdk1(I)
  967.            JNS(2*NT+I,3*NT+1) = dRxDdk2(I)
  968.            JNS(     I,3*NT+2) = Rw(I) 
  969.            JNS(  NT+I,3*NT+2) = RxwPhiR(I)
  970.            JNS(2*NT+I,3*NT+2) = RxwPhiI(I)
  971.            JNS(3*NT+1,I)      = ZERO
  972.            JNS(3*NT+1,NT+I)   = VQ(I)
  973.            JNS(3*NT+1,2*NT+I) = ZERO           
  974.            JNS(3*NT+2,I)      = ZERO
  975.            JNS(3*NT+2,NT+I)   = ZERO
  976.            JNS(3*NT+2,2*NT+I) = VQ(I)
  977.            solNS(I)           = -fvec0(I)
  978.            solNS(NT+I)        = -solNS(NT+I)
  979.            solNS(2*NT+I)      = -solNS(2*NT+I)
  980.          ENDDO
  981.  
  982.          solNS(3*NT+1) = -(DDOT(NT,VQ,1,PHIR,1)-ONE)
  983.          solNS(3*NT+2) = -(DDOT(NT,VQ,1,PHII,1))
  984.          JNS(3*NT+1,3*NT+1) = ZERO
  985.          JNS(3*NT+1,3*NT+2) = ZERO
  986.          JNS(3*NT+2,3*NT+1) = ZERO
  987.          JNS(3*NT+2,3*NT+2) = ZERO 
  988. *        Resolution
  989.          CALL DGESV(3*NT+2,1,JNS,3*NT+2,I4,solNS,3*NT+2,INFO) 
  990.          DO IOU=1,3*NT+2
  991.          ENDDO
  992. *        Corrections 
  993.          dWB = solNS(3*NT+2)
  994.          dK  = solNS(3*NT+1)
  995.          DO I = 1,NT
  996.           dXB(I) = solNS(I)
  997.           dPHIR(I)= solNS(NT+I)    
  998.           dPHII(I)= solNS(2*NT+I)    
  999.           Q1(I) = Q1(I) + dXB(I)
  1000.           PHIR(I) = PHIR(I) + dPHIR(I)
  1001.           PHII(I) = PHII(I) + dPHII(I)
  1002.          ENDDO
  1003.          OMEG = OMEG + dWB     
  1004.          KAPPA= KAPPA+ dK
  1005. *         write(*,*) 'w +',dWB,'=',OMEG
  1006. *         write(*,*) 'kappa +',dK,'=',KAPPA
  1007.          END IF
  1008. *=======================================================================     
  1009. *          Evaluation du residu       
  1010. *======================================================================= 
  1011.  
  1012. *          a) R
  1013.          CALL funcv(NT,NHBM,NDDL,NFFT,MTQ,MTKAM,MTPHI,MTLIAA,
  1014.      &              MTLIAB,MTPAS,MTFEX,LOCLFA,LOCLB1,CANAL,fvec0,fref,
  1015.      &            .false.)     
  1016.  
  1017. *-----------------------------------------------------------------------
  1018.          IF (TYPC.EQ.'L') THEN
  1019. *-----------------------------------------------------------------------
  1020. *          b) Rxs*Phi     
  1021.            CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)  
  1022.            DO I=1,NT
  1023.              res2(I) = ZERO
  1024.              DO J=1,NT
  1025.                RX(I,J) = ZZ(I,J) - JAC(I,J)
  1026.                res2(I) = res2(I) + RX(I,J)*PHIR(J)
  1027.              ENDDO         
  1028.            ENDDO
  1029. *          c) norm(Phi)= 1.
  1030.            res3 = DDOT(NT,PHIR,1,PHIR,1)-ONE
  1031. *          Y = ( R  ;  RQ * \phi  ;  \phi*\phi-1 )
  1032.            arrRes(1) = dnrm2(NT,fvec0,1)
  1033.            arrRes(2) = dnrm2(NT,res2,1)
  1034.            arrRes(3) = abs(res3)
  1035.            res = MAXVAL(arrRes,1)
  1036.            IF (IIMPI.GE.3)
  1037.      &     WRITE(IOIMP,991) its,(arrRes(iou),iou=1,3),res
  1038.  991       FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';',
  1039.      &     E10.3,';','} -> res=',E10.3)        
  1040.          ENDIF
  1041.  
  1042. *-----------------------------------------------------------------------
  1043.          IF (TYPC.EQ.'B') THEN
  1044. *-----------------------------------------------------------------------        
  1045. *          b) Rxs*Phi 
  1046.            CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)  
  1047.            CALL HBMDVEC(NT,NHBM,NDDL,Q1,ONE,V1)    
  1048.            znormV1 = dnrm2(NT,V1,1)
  1049.            DO I = 1,NT 
  1050.              V1(I) = V1(I)/znormV1
  1051.            ENDDO
  1052.            DO I = 1,NT
  1053.              res2(I) = ZERO
  1054.              DO J = 1,NT
  1055.                RX(I,J) = ZZ(I,J) - JAC(I,J) 
  1056.                RX2(I,J) = RX(I,J) + V1(I)*V1(J) 
  1057.                res2(I) = res2(I) + RX2(I,J)*PHIR(J)
  1058.              ENDDO
  1059.            ENDDO    
  1060. *          c) Rw'*Phi
  1061.            CALL HBMRW(NT,NDDL,1,Q1,OMEG,XM,XASM,XK,Rw,.true.)  
  1062.            IF (BAL.EQ.1) THEN
  1063.              DO I=1,NT
  1064.                Rw(I) = Rw(I) - TWO*OMEG*FEXA(I)
  1065.              ENDDO
  1066.            END IF 
  1067.            res3 = DDOT(NT,Rw,1,PHIR,1)
  1068. *          d) norm(Phi) = 1
  1069.            res4 = DDOT(NT,PHIR,1,PHIR,1)-ONE
  1070. *          Y = ( R  ;  RQ * \phi  ;  Rw\phi  ;  \phi*\phi-1 )
  1071.            arrRes2(1) = dnrm2(NT,fvec0,1)
  1072.            arrRes2(2) = dnrm2(NT,res2,1)
  1073.            arrRes2(3) = abs(res3)
  1074.            arrRes2(4) = abs(res4)
  1075.            res = MAXVAL(arrRes2,1)           
  1076.            IF (IIMPI.GE.3)
  1077.      &     WRITE(IOIMP,992) its,(arrRes2(iou),iou=1,4),res
  1078.  992       FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';',
  1079.      &     E10.3,';',E10.3,'} -> res=',E10.3)        
  1080.          ENDIF
  1081.  
  1082. *-----------------------------------------------------------------------
  1083.          IF (TYPC.EQ.'P') THEN
  1084. *-----------------------------------------------------------------------        
  1085.            CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)  
  1086.            DO I = 1,NT
  1087.              DO J = 1,NT
  1088.                RX(I,J) = ZZ(I,J) - JAC(I,J) 
  1089.              ENDDO
  1090.            ENDDO    
  1091.            CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D)
  1092.            CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1)
  1093. *          b) (Rx-((OMEG/2)**2)*D2)*PhiR - KAPPA*D1*PhiI     
  1094. *          c) (Rx-((OMEG/2)**2)*D2)*PhiI + KAPPA*D1*PhiR   
  1095.            DO I = 1,NT
  1096.              resV1(I) = ZERO 
  1097.              resV2(I) = ZERO
  1098.              DO J = 1,NT 
  1099.                resV1(I)=resV1(I)+Rxmk2D(I,J)*PHIR(J)
  1100.      &                          -KAPPA*DEL1(I,J)*PHII(J)
  1101.                resV2(I)=resV2(I)+Rxmk2D(I,J)*PHII(J)
  1102.      &                          +KAPPA*DEL1(I,J)*PHIR(J)
  1103.              ENDDO
  1104.            ENDDO
  1105. *          d) (VQ'*PHIR) - 1
  1106.            res4 = DDOT(NT,VQ,1,PHIR,1)-ONE
  1107. *          Y = ( R  ;  [D0-k²*D2]*\phi_R - k*D1*\phi_I  ; 
  1108. *                      [D0-k²*D2]*\phi_I + k*D1*\phi_R  ; Q*\phi_I-1 )
  1109.            arrRes2(1) = dnrm2(NT,fvec0,1)
  1110.            arrRes2(2) = dnrm2(NT,resV1,1)
  1111.            arrRes2(3) = dnrm2(NT,resV2,1)
  1112.            arrRes2(4) = abs(res4)
  1113.            res = MAXVAL(arrRes2,1)           
  1114.            IF (IIMPI.GE.3)
  1115.      &     WRITE(IOIMP,993) its,(arrRes2(iou),iou=1,4),res
  1116.  993       FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';',
  1117.      &     E10.3,';',E10.3,'} -> res=',E10.3)        
  1118.          ENDIF
  1119.  
  1120. *-----------------------------------------------------------------------
  1121.          IF (TYPC.EQ.'N') THEN
  1122. *-----------------------------------------------------------------------        
  1123.            CALL HBMZ(NT,NHBM,NDDL,MTQ,MTKAM,.true.)  
  1124.            DO J = 1,NT
  1125.              DO I = 1,NT
  1126.                RX(I,J) = ZZ(I,J) - JAC(I,J) 
  1127.              ENDDO
  1128.            ENDDO    
  1129.            CALL COLIMA(NT,NT,RX,D2,ONE,-1.D0*KAPPA**2,Rxmk2D)
  1130.            CALL COLIMA(NT,NT,D10,D1w,ONE,OMEG,DEL1)
  1131. *          b) (Rx-(KAPPA**2)*D2)*PhiR - KAPPA*D1*PhiI     
  1132. *          c) (Rx-(KAPPA**2)*D2)*PhiI + KAPPA*D1*PhiR     
  1133.            DO I = 1,NT
  1134.              resV1(I) = ZERO 
  1135.              resV2(I) = ZERO
  1136.              DO J = 1,NT 
  1137.                resV1(I)=resV1(I)+Rxmk2D(I,J)*PHIR(J)
  1138.      &                          -KAPPA*DEL1(I,J)*PHII(J)
  1139.                resV2(I)=resV2(I)+Rxmk2D(I,J)*PHII(J)
  1140.      &                          +KAPPA*DEL1(I,J)*PHIR(J)
  1141.              ENDDO
  1142.            ENDDO
  1143. *          d) (VQ'*PHIR) - 1
  1144.            res4 = DDOT(NT,VQ,1,PHIR,1)-ONE
  1145. *          e) VQ'*PHII
  1146.            res5 = DDOT(NT,VQ,1,PHII,1)
  1147. *          Y = ( R  ;  [D0-k²*D2]*\phi_R - k*D1*\phi_I  ; 
  1148. *                      [D0-k²*D2]*\phi_I + k*D1*\phi_R  ; 
  1149. *                                            vQ*\phi_R-1 ; vQ*\phi_I)
  1150.            arrRes3(1) = dnrm2(NT,fvec0,1)
  1151.            arrRes3(2) = dnrm2(NT,resV1,1)
  1152.            arrRes3(3) = dnrm2(NT,resV2,1)
  1153.            arrRes3(4) = abs(res4)
  1154.            arrRes3(5) = abs(res5)
  1155.            res = MAXVAL(arrRes3,1)            
  1156.            IF (IIMPI.GE.3)
  1157.      &     WRITE(IOIMP,994) its,(arrRes3(iou),iou=1,5),res
  1158. c      &     WRITE(IOIMP,994) its,KAPPA,(arrRes3(iou),iou=1,4),res
  1159.  994       FORMAT('+ hbmbif : iter ',I2,' Y={',E10.3,';',E10.3,';',
  1160.      &     E10.3,';',E10.3,';',E10.3,'} -> res=',E10.3)        
  1161.          ENDIF 
  1162. *-----------------------------------------------------------------------        
  1163. c            IF (IIMPI.GE.3) WRITE(IOIMP,*) 'Iter ',its,' -> res =',res
  1164.            IF (res.LE.TOLMIN) THEN             
  1165.              CHECK=.false.
  1166.              CBIF = CBIF + 1
  1167.              TYPBIF(CBIF) = TYPC
  1168.              DO I=1,NT
  1169.                QBIFU(I,CBIF) = Q1(I)
  1170.                QPSIR(I,CBIF) = PHIR(I)
  1171.                QPSII(I,CBIF) = PHII(I)
  1172.                Q1(NT) = QIN(I)
  1173.              ENDDO
  1174.              WBIFU(CBIF) = OMEG 
  1175.              WBIF2(CBIF) = KAPPA
  1176.              OMEG = OMEGIN
  1177. c              RETURN
  1178.              GOTO 666
  1179.            ELSE
  1180.              CHECK = .true.
  1181.            ENDIF           
  1182.        ENDDO 
  1183.        DO I = 1,NT 
  1184.          Q1(I) = QIN(I)
  1185.        ENDDO
  1186.        OMEG = OMEGIN
  1187. c      message
  1188.        IF (CHECK) THEN
  1189. c          Pas de convergence apres %i1 iterations. L'execution continue
  1190.            INTERR(1)=MAXITS
  1191.            CALL ERREUR(151) 
  1192. c            RETURN
  1193.              GOTO 666
  1194.        ELSEIF (IIMPI.GE.2) THEN
  1195.            WRITE(IOIMP,*) '--> convergence apres ',its,'iterations'
  1196.        ENDIF 
  1197. *      MENAGE       
  1198.  666   CONTINUE
  1199. c        SEGSUP,MWORK
  1200.        END       
  1201.  
  1202.  
  1203.  
  1204.  
  1205.  

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